search

choldgraf / ecogtools

A collection of tools for analyzing electrocorticography data
BSD 2-Clause "Simplified" License
17 stars 4 forks source link

Q?: pipeline gamma from raw #1

Open kingjr opened 8 years ago

kingjr commented 8 years ago

Hey :)

Do you have a pipeline for extract gamma activity from the raw data, instead of doing it at the epoch level?

choldgraf commented 8 years ago

Hey! Yep, the function I wrote for this is here:

https://github.com/choldgraf/ecogtools/blob/master/tfr.py#L244

It lets you specify a few sub-bands if you want. It'll do a bandpass for each, and will average them together (scaling them first, optionally).

kingjr commented 8 years ago

Thanks Chris. I wrote this to make it memory/CPU efficient; it works with multitapers and morlet, I can add the hilbert too; it's probably going to be fast too.

The main design change is that

Also, in my Matlab code I transformed the power in log base before normalizing. But I didn't write why, nor a ref... So I'm not sure it was a good idea in the first place

def extract_amplitude(inst, freqs, normalize=True, method='morlet',
                      picks=None, log=False, copy=False, split_size=1):
    from tqdm import tqdm  # progress bar
    from mne.io.pick import _pick_data_channels
    from mne.time_frequency import cwt_morlet
    from mne.time_frequency.tfr import _induced_power_mtm
    from jr.stats import median_abs_deviation

    if picks is None:
        picks = _pick_data_channels(inst.info, with_ref_meg=False)
    if copy:
        inst = inst.copy()

    data = inst._data[picks, :]
    sfreq = inst.info['sfreq']

    # Loop across chunks of channels to optimize memory/CPU
    chans_splits = np.array_split(np.arange(len(picks)),
                                  np.ceil(float(len(picks)) / split_size))
    for chans in tqdm(chans_splits):
        if method == 'morlet':
            x = cwt_morlet(data[chans, :], sfreq, freqs,
                           use_fft=True, n_cycles=7.0, zero_mean=False)
            power = (x * x.conj()).real
            del x
        elif method == 'mtm':
            power, itc = _induced_power_mtm(data[None, chans, :], sfreq=sfreq,
                                            frequencies=freqs, n_jobs=-1,
                                            zero_mean=True, verbose='INFO')
            del itc
        else:
            raise NotImplementedError('Unknown method %s' % method)
        if log:
            np.log(1 + power, out=power)  # in place
        if normalize:
            mean = np.median(power, axis=2)[:, :, np.newaxis]
            std = median_abs_deviation(power, axis=2)[:, :, np.newaxis]
            std[std == 0] = 1.
            power -= mean
            power /= std
        mean_power = np.mean(power, axis=1)
        inst._data[picks[chans], :] = mean_power
        del power, mean_power
    return inst

e.g.

freqs = np.logspace(np.log10(60), np.log10(150), 10)
raw = extract_amplitude(raw, freqs, method='morlet', split_size=10)
kingjr commented 8 years ago

for from jr.stats import median_abs_deviation see https://github.com/kingjr/jr-tools/blob/master/jr/stats/base.py#L484

choldgraf commented 8 years ago

Nice! Thanks for the updated function. I've seen papers do log-transforming and others that don't, I don't think there's a rule one way or another. In my case, I usually z-score data at some point anyway (within frequency bands) so I haven't worried too much about it.

Do you have any interest in a co-maintained set of ecog tools? We could create a new repo for it and I could port over the relevant code from this repository into that one. I'd be happy to come up with functions that are crowdsourced for usefulness and correctness (well, even if the crowd is an N of 2). I'm imagining functions that are useful, but maybe too specific to contribute to MNE. Or maybe it'd be too specific? I just imagine that there are lots of ecog pythonistas out there implementing their own functions, no reason we should all be reinventing the wheel.

kingjr commented 8 years ago

Absolutely. I'm happy to contribute here directly, no need to port it elsewhere (except if you want to keep a dirty repo for yourself.). Eventually, I can transfer ecoggui for the elec location here too. (I'm waiting to get my next patient before looking at it again, hence the waiting)

I think the good move would be to use this as a buffer before MNE. Once we have a bit of experience with the functions we wrote we can propose to make them available to MNE with an example for each. WDYT?

choldgraf commented 8 years ago

I think that's a good idea - it's pretty much how I've tried to structure code here already (using MNE docstring style as much as possible and structuring objects similarly etc). Want me to just give you admin rights? If there are functions in here that are too specific to my stuff then I can put that in another repository.

kingjr commented 8 years ago

Sure. Shall we set up Travis ?

On Friday, 11 March 2016, Chris Holdgraf notifications@github.com wrote:

I think that's a good idea - it's pretty much how I've tried to structure code here already (using MNE docstring style as much as possible and structuring objects similarly etc). Want me to just give you admin rights? If there are functions in here that are too specific to my stuff then I can put that in another repository.

— Reply to this email directly or view it on GitHub https://github.com/choldgraf/ecogtools/issues/1#issuecomment-195631624.

choldgraf commented 8 years ago

Sounds like a good idea - I don't have much experience setting Travis up but it's something I've been wanting to figure out for a while now. I assume it is basically just setting up a random number generator that has a 99% probability of telling you that your code is wrong on each query, is that correct? ;)

On Fri, Mar 11, 2016 at 8:41 PM, Jean-Rémi KING notifications@github.com wrote:

Sure. Shall we set up Travis ?

On Friday, 11 March 2016, Chris Holdgraf notifications@github.com wrote:

I think that's a good idea - it's pretty much how I've tried to structure code here already (using MNE docstring style as much as possible and structuring objects similarly etc). Want me to just give you admin rights? If there are functions in here that are too specific to my stuff then I can put that in another repository.

— Reply to this email directly or view it on GitHub <https://github.com/choldgraf/ecogtools/issues/1#issuecomment-195631624 .

— Reply to this email directly or view it on GitHub https://github.com/choldgraf/ecogtools/issues/1#issuecomment-195659534.

kingjr commented 8 years ago

I actually don't know either. Good place to learn.

In case you haven't seen, they just created an mne-sandbox rep. I'm not sure whether this wouldn't be a good place too. On the plus side we'd have their feedbacks right away. On the minus one, it'll be heavier (more testing, more integration with all types of MEG specificies, e.g. projections etc. WDYT?

choldgraf commented 8 years ago

Hey - jumping back on this now that graduation craziness has ended. I am a bit torn between having an ecog-specific repository and pushing things to MNE-sandbox. I guess the question is where there is enough ecog-specific stuff to warrant its own repository. If we think that's the case, then we can create an MNE-ecog repo or something like this. However, IMO if there's significant overlap with EEG/MEG then I'm happy to just contribute this to MNE-sandbox/python. It'd be cool if we could improve functionality in general for eeg / ieeg now that both types are supported natively.

So I guess my intuition is that we should make an mne_sandbox.ecog submodule and then build things in there.

We should also figure out a good ECoG dataset to use for testing. Maybe we can take an existing dataset, anonymize it, remove metadata about the task except for when "events" happen, and choose a bunch of semi-random points on the MNI brain that roughly reflect the grid coverage.