WIP epochs.plot_image(gfp=True)

[jona-sassenhagen]

I'm still annoyed by the fact that epochs.plot_image throws so many figures at you.

This allows one to plot the single-trial GFP.

[codecov-io]

Codecov Report

Merging #4227 into master will decrease coverage by 0.04%.
The diff coverage is 90.26%.

@@            Coverage Diff            @@
##           master   #4227      +/-   ##
=========================================
- Coverage   86.24%   86.2%   -0.05%     
=========================================
  Files         357     358       +1     
  Lines       64566   65037     +471     
  Branches     9831    9922      +91     
=========================================
+ Hits        55685   56064     +379     
- Misses       6180    6237      +57     
- Partials     2701    2736      +35

[jona-sassenhagen]

What do people think? Do we even want an epochs.plot_image(gfp=True) method?

I like it because I like "survey"-style plots that give you a super fast impression when called with defaults, and otherwise you always have to specify picks and you don't get the full answer either.

Also we should have epochs.plot_joint_image() (e.g. erpimage with joint topoplots), and for that we need this anyways.

[agramfort]

no objection with the gfp image

[jona-sassenhagen]

Alternative API suggestion, inspired by the tfr.plot_joint API:

epochs.plot_image(aggregate='gfp')
epochs.plot_image(aggregate=np.std)
epochs.plot_image(aggregate=np.mean)
epochs.plot_image(aggregate=np.median)

(vs. epochs.plot_image(gfp=True))

Thoughts?

[teonbrooks]

could you add an image of the plot when you have a chance? would be cool to check it out

[teonbrooks]

what does none yield?

[jona-sassenhagen]

Ok, need to doc this ...
The point is, for the GFP we don't need a diverging cmap (None -> Reds), so I've switched from the RdBu_r default to None.

[teonbrooks]

just curious, why on the average of epochs instead of the average of individual GFPs?

[jona-sassenhagen]

No reason I think. I wrote the code this way without thinking, then realized that there are two possible options, and decided to just doc the one I chose. If you think it should go the other way around, I'll do it like that.

I should probably reword it though. It's not the most obvious way of saying "data.std(ch_ax).mean(trial_ax) vs. data.mean(trial_ax).std(ch_ax)".

[teonbrooks]

"data.std(ch_ax).mean(trial_ax) vs. data.mean(trial_ax).std(ch_ax)".

that clarifies it much better

[mmagnuski]

I like aggregate option - would use it with picks often. I prefer reduction as the kwarg name though (which reminds me of my neglected welch reduction PR :) )

[agramfort]

YAGNI? will you really use more than GFP? Now on the positive side and think it's more explicit than GFP jargon :)

[jona-sassenhagen]

@agramfort I sometimes stuff model fits or similar stuff into mne insts for plotting convenience. If your values are e.g. r or R^2, the std doesn't make so much sense anymore, but maybe the median or some other metric would.

Or what if the data isn't averaged to a common ref anymore, but still on a common scale (e.g. the bipolar stuff)? Essentially a bunch of data transformations might require a different measure.

Not saying it's super important, but it's only like 3 lines of code, 2 lines of doc and 2 lines of test so ...

[jona-sassenhagen]

... any chance of having this one to actually default to gfp=True (or aggregate='gfp') in the future?

Reason: if called without specifying picks, the user gets 5 figures. That's a lot of plots, and the first 5 are usually far from the most important. The GFP is only 1-3 figures, and typically a lot more informative.

I doubt it would hurt any existing pipelines, because these probably don't rely on leaving picks unspecified and we changed the behavior for this case anyways a few releases ago.

[jona-sassenhagen]

@teonbrooks

(This is a very high SNR data set though.)

Actually, "inferno" or 'viridis' look much better than "Reds" for a lot of my tests

[jona-sassenhagen]

Votes on:

1. defaulting epochs.plot_image to plotting the gfp if no picks are given, yes/no?

2. having an aggregate param instead of a gfp param, yes/no?

3. name for the aggregate param? aggregate vs. reduce vs. ..? (In any case, if [MRG] Adding AverageTFR.plot_joint() method #4134 is merged first, I will follow its lead!)

I'm for the 1st option in all cases.

[agramfort]

you should not use std as you you don't want to subtract the mean for GFP

[agramfort]

giving a try to this PR gives me:

In [2]: epochs.plot_image(gfp=True)
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-2-a38d371c8036> in <module>()
----> 1 epochs.plot_image(gfp=True)

/Users/alex/work/src/mne-python/mne/epochs.py in plot_image(self, picks, sigma, vmin, vmax, colorbar, order, show, units, scalings, cmap, fig, axes, overlay_times, gfp)
   1009                                  show=show, units=units, scalings=scalings,
   1010                                  cmap=cmap, fig=fig, axes=axes,
-> 1011                                  overlay_times=overlay_times, gfp=gfp)
   1012
   1013     @verbose

/Users/alex/work/src/mne-python/mne/viz/epochs.py in plot_epochs_image(epochs, picks, sigma, vmin, vmax, colorbar, order, show, units, scalings, cmap, fig, axes, overlay_times, gfp)
    138             data_ = data[:, picks_, :].std(1)
    139             data_per_type.append(data_)
--> 140             evoked_data = scalings[ch_type] * evoked.data[picks_, :].std(0)
    141             evokeds.append(evoked_data)
    142         data = np.asarray(data_per_type)

IndexError: index 366 is out of bounds for axis 0 with size 366

I just tried

epochs.plot_image(gfp=True)

after running %run plot_visualize_epochs.py

[agramfort]

how about a group_by param as we have in raw.plot now? picks will be used to select what channels to consider and if group_by=None then show individual channels and if group_by='type' show GFP by channel type. But we could have group_by='selection' wdyt?

[jona-sassenhagen]

group_by='selection' would be for interactive mode?

[agramfort]

no it would give you the GFP image per selection of channels like raw.plot does it to make the butterfly rendering by selection

[schekroud]

Sorry if misunderstanding conv up to now but was directed here to provide input if possible (and doesn't seem to be in same train of thought? could be wrong)

is it possible to have a: epochs.plot_image(average_over_picks = bool) param where instead of plotting all chans in picks, you plot one image of average voltage/field at all chosen channels? sometimes valuable if looking at broader cortical (e.g. visual) response ?
would require iterating over epochs to calculate evoked for average of channels, then plotting in exact same way as current function?

something like:

e.g. use:
mne.epochs.plot_image(epochs, picks = picks, average_over_picks=True) and do something like this within the function....

n_epochs = np.shape(epochs)[0] # number of epochs in struct
selected = []
for i in range(len(n_epochs)):
    i_evoked = epochs[i].average() #has a shape of [nchans,timesamples]
    picked = i_evoked.data[picks,:] #get data from chosen channels
    chanaverage = np.nanmean(picked,0) #average across chosen channels
    selected.append(chanaverage) #add single epoch data, averaged across chosen chans -> new array

I could be on the wrong path here, but would create new array of single epoch data with only one 'channel' (average of picked channels) which could be stored back as an EpochArray to feed into the original epochs.plot_image code?

[jona-sassenhagen]

You're essentially asking for the "reduce" or "aggregate" option proposed upthread. Although usually the GFP is a better measure than the mean for eeg data ...

[schekroud]

ah ok, makes sense. mean would keep consistency with other non-python packages like FieldTrip, no? useful if people are working in both alongside each other as a sanity check (this is happening in my lab as a couple of us are trying to move people away from matlab...)

[jona-sassenhagen]

I'm sure if we can convince a few MATLAB users to switch with that, it'll be possible to implement this :)

[jona-sassenhagen]

@agramfort I finally tried the group_by thing in plot_raw and 1. we should have this for 10/20 data too (it's a predictable mapping from name to ROI) and I'd like to do that eventually, 2. i think that should be parallel to the aggregation.

partially because if you have a ROI-like parameter, grouping by ROI means aggregating via the mean actually makes sense and is very common. E.g.,

evoked.plot_image(group_by=my_ROI_dict, aggregate="mean")

Also you have to admit switching between grouping channels to aggregating over them via the GFP is really overloading the group_by keyword :)

[jona-sassenhagen]

Btw @schekroud don't use this style in Python:

n_epochs = np.shape(epochs)[0] # number of epochs in struct
selected = []
for i in range(len(n_epochs)):
    i_evoked = epochs[i].average() #has a shape of [nchans,timesamples]
    ...
    selected.append(some_function(i_evoked))

Almost always, you can do something like this:

for item in iterable_of_items:
    do_something_with(item)

And usually

collection = [do_something_with(item) for item in iterable_of_items]

That's much more readable!

In this case, I think you want something like:

[evoked.data.mean(0) for evoked in epochs.iter_evoked()]

or just straight-forward

epochs.get_data().mean(1)

This would give you the data averaged over channels.

My idea for this function is to have:

epochs.plot_image(aggregate=agg_kind)

where agg_kind can be False (for current behavior), 'gfp', or e.g. a function (e.g. np.mean) that will be applied to aggregate over channels.

[agramfort]

Also you have to admit switching between grouping channels to aggregating over them via the GFP is really overloading the group_by keyword :) what do you mean? I am not sure about the name aggregate but otherwise I

think the API is nice.

[jona-sassenhagen]

You mean you'd be okay with something like:

evoked.plot_image(group_by=my_ROI_dict, reduce="gfp")

evoked.plot_image(group_by='type', reduce=np.mean)

?

[agramfort]

yes but maybe aggregate is better? no strong feeling. I feel that aggregate is really geeky...

[schekroud]

thanks for the tips - only been python-ing for a (relatively) short period so brevity of code is an issue...

maybe grouping_measure is a better (read: less geeky) term than aggregate?
epochs.plot_image(group_by = ROI_dict, grouping_measure = 'gfp')
epochs.plot_image(group_by = ROI_dict, grouping_measure = np.mean)

[jona-sassenhagen]

@agramfort I think I'm the only one who prefers aggregate :)

We have mne.combine_evokeds, so combine is another option.

@schekroud grouping_measure is too long IMO and also the wrong term- the grouping separates channels from each other, the aggregation reduces them.

[agramfort]

+1 for combine it also matches the split-apply-combine approach of pandas

[teonbrooks]

My 2c rankings:

1. aggregate
2. combine
3. groupby

[agramfort]

we need both group_by and combine

[jona-sassenhagen]

Oh god, we need some consensus here :) So what do we call this param? @jaeilepp maybe you have an opinion?

[jaeilepp]

My vote goes for combine.

[jaeilepp]

What is going on here?

[jona-sassenhagen]

Don't review the code right now, I've been putting it on hold until we decide on the API.

[jona-sassenhagen]

Ok, then I'll do combine unless anybody has strong objections.

[agramfort]

yes but we agree that we need both group_by and combine right?

[jona-sassenhagen]

Should it default to None or to 'type' though? I guess 'type'?

[mmagnuski]

no strong feelings, combine is ok. I'll try to revive my welch-reduction PR in a few days and I'll change reduction to combine then.

[larsoner]

I do like "reduce" because it's standard to talk about map-reduce operations, but I think "combine" is more common in MNE so +1 for that from me

[jona-sassenhagen]

I'm working on this again ... next thing I'll do.

[jona-sassenhagen]

Did something very weird. I'll open a new PR.