Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[MRG] plot_ica_properties #3275

Merged
merged 62 commits into from
Jun 27, 2016
Merged

Conversation

mmagnuski
Copy link
Member

@mmagnuski mmagnuski commented Jun 1, 2016

Adressess #3269, #2747 by adding plot_ica_properties (and probably a method for ICA class etc.).

Example

An example of plot_ica_properties (name is open for discussion of course) in inline mode (ep is an instance of mne.Epochs here):

It is definitelly more cluttered than when plotted in a separate window. I will work on this (probably just change figure size). ## Discussion

@jona-sassenhagen @agramfort @espressofiend @Eric89GXL @jaeilepp
I will need your input on several issues:

  • API:
    • name of the base function / ica method
    • input arguments (for example: topo_kwargs to allow the user to control how the topoplot is made?). Some useful candidates (with rather temporary names):
      • dB - to plot spectrum in dB
      • cmap - to control topoplot and epochs image colormap
      • interactive - to run in interactive mode or not (if we generate non-interactive image, we do not need buttons etc., but see also point 2)
  • PropertiesExplorer (or some other name) class for future PR adding interactivity. I know you prefer not to add too many classes to mne-python, but I think it would be useful in this case. If you are ok with adding such a class I could start working on it in this PR. Some reasons for why I think it would be useful:
    • it is easier to write/maintain code for interactive visualizations if all the data and plot refresh methods are kept in one object and not in partial functions

    • it is also easier to extend the interactive plotter (from a user and developer perspective) either by subclassing it or attaching callbacks to the figure/axis that take the plotter object as one of the arguments. If all relevant data is kept in partial functions it is much harder to use it without rewriting the original mne code.
      Example: if I wanted to extend the properties plotter to show the unmixing matrix weights as a bar plot when I click on the topo map I could just do:

      prop = ica.plot_properties(epochs, interactive=True, show=False)
      def my_callback(event, obj=prop):
        if event.inaxes.get_label() == 'topo':
            fig, ax = plt.subplots()
            wgt = obj.ic_weights[obj.ic_idx]
            ax.bar(range(len(wgt)), wgt)
      prop.fig.canvas.mpl_connect('button_press_event', my_callback)
      prop.show()

      The discussion on the class can wait until my next PR if you prefer.

  • and just BTW - I would be really glad for some tips on how you reload mne when developing. I was testing the code in Jupyter notebook and restarting the kernel every time (and recalculating everything including ica) can be a bit tiresome. I tired using importlib.reload but it didn't seem to work... How do you reload mne to test incremental changes?
  • I will also need some help with making this PR work with channel types other than "eeg" - I have no experience in how mne handles MEG sensors

TODOs

  • fix spectrum std (power distribution is skewed, symmetric std doesn't work well)
  • make this work with continuous signal (segment raw data into 1-second (or something like that) "epochs")
  • decrease clutter so that it looks better with inline mode
  • plot_properties arguments
  • make topoplot work also with channel types other than eeg (I will need help here)
  • add tests
  • segment method for _BaseRaw
  • psd_kws or freq_kws
  • should picks be allowed a list of values - creating multiple figures ?
  • check a problem with two figures poping up when only one should be present
  • sigma param?
  • check the warning with non-epoched data - RuntimeWarning: Mean of empty slice.
  • update ICA example
  • tests for _segment_raw?
  • speed up ica.plot_properties tests:
  • default value for fmax
  • update what's new

@mmagnuski
Copy link
Member Author

Changing figure size helps with clutter:

I tried detecting whether "time" and "epoch variance" text objects are too close, but it is harder than I expected - for example textobj.get_window_extent() returns bounding boxes with zero width and zero height 😞

@jona-sassenhagen
Copy link
Contributor

Looks pretty good.

@@ -3106,3 +3106,32 @@ def average_movements(epochs, head_pos=None, orig_sfreq=None, picks=None,
_remove_meg_projs(evoked) # remove MEG projectors, they won't apply now
logger.info('Created Evoked dataset from %s epochs' % (count,))
return (evoked, mapping) if return_mapping else evoked


def segment_raw(raw, segment_length=1., verbose=True):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be exposed or private?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't have a strong opinion here, I quite often use segmentation but maybe it is not common enough to expose it. I can add the _ prefix if you prefer.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

let's make it private

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I use something like this ever so often too, actually. But I'm okay with either.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes please make it private and don't set verbose to True by default. See how we use the verbose decorator elsewhere

@jona-sassenhagen
Copy link
Contributor

Consider using this for the topoplots: #3165
It should make it easy to do MEG.

axes[2].plot(1e3 * evoked.times, evoked.data[comp_idx].ravel())
# spectrum
axes[3].plot(freqs, psds_mean, color='k')
axes[3].fill_between(freqs, psds_mean - neg_sd, psds_mean + pos_sd,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Another option would be to use descriptive names (erp_axes, topomap_axes, ...).

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am definitelly going to add labels to the axes (will be useful for interactive stuff), but I can change the variable names here too.

@jona-sassenhagen
Copy link
Contributor

Oh, and please also make it work for regular Epochs :) (like EEGLAB)

@mmagnuski
Copy link
Member Author

I thought it does... or maybe you mean something else by "regular Epochs"?
And thanks for the plot_topomap suggestion - much more convenient.

raise ValueError('inst should be an instance of Raw or Epochs,'
' got %s instead.' % type(inst))
if isinstance(picks, int):
picks = [picks]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're doing this twice.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

😄

@jona-sassenhagen
Copy link
Contributor

I mean channel data.

@mmagnuski
Copy link
Member Author

Ah, right, I was thinking about another PR, but I can do it here if you insist :)
So how would you like it to work? plot_properties method for Raw, Epochs and ICA where component properties are only plotted with ica.plot_properties(inst)?

evoked = src.average(picks)
# smooth_data is kept in a separate variable for future
# interactive changing of sigma
smooth_data = ndimage.gaussian_filter1d(ica_data[comp_idx],
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Have you thought about hacking plot_epochs_image so that its individual axes are somehow accessible instead?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I thought it might be too much work but it might be useful if we want to have topo_kwargs and epochs_image_kwargs for example...
plot_epochs_image creates a separate figure so I would need to make it accept axes keyword for example.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Most of the plot functions have axes params so I think this might be good.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also we're big into DRY so I think it should be done somehow.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, I'll try to look into it today/tomorrow.

@jona-sassenhagen
Copy link
Contributor

jona-sassenhagen commented Jun 3, 2016

Ah, right, I was thinking about another PR, but I can do it here if you insist :)

No, it's okay, just keep it in mind - don't make it impossible to add that functionality later.

seg_len = np.arange(0.5, 10., 0.5)
nice_len = inst._data.shape[1] / (250. * inst.info['sfreq'])
seg_len = seg_len[np.argmin(np.abs(seg_len - nice_len))]
inst = segment_raw(inst, segment_length=seg_len, verbose=False)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

don't like this very much. 1s is arbitrary + all the various Epochs params. I suggest you do not show any raster plot for Raw data. Just when looking at properties of Epochs instances. I means though that figure layout will be different for Raw and Epochs.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is what EEGLAB does and it's actually very useful. It's a very good way of showing what and where particular artefacts are, what kind of activity a channel captures.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I also like what eeglab does for continuous files in pop_prop. It is not as informative as for epoched data but you get some "birds eye view" instead of nothing.
BTW - I do not set segmentation to 1s - the code above looks for a value from range np.arange(0.5, 10., 0.5) that is closest to segmet length that divides that data into 250 chunks. 250 is arbitrary though :)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

IMO least surprise would be using EEGLAB's default.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok you win :)
Make the segment function private though please

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jona-sassenhagen by EEGLAB's default you mean 1 second length segments?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jona-sassenhagen by EEGLAB's default you mean 1 second length segments?

I think EEGLAB uses 2 seconds?

Predictability is probably better than a smart/'smart' auto magic here.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not perfectly sure here because it can mean that something looks predictably worse :), but I didn't check it thoroughly anyway so the simple 2s behavior should be ok.

@mmagnuski
Copy link
Member Author

mmagnuski commented Jun 4, 2016

axes argument in plot_epochs_image

After discussing this with @jona-sassenhagen I went ahead and added axes argument to plot_epochs_image so that this function can be reused here. So now one can create abominations like this:

A couple of issues:

  • regardless of points below adding axes to plot_epochs_image increases its reusability. The docstring suggested so far that passing a figure to plot_epochs_image causes it to fill two axes of the figure but actually the axes were created. Now if axes list is passed the axes are being used - even if one of the axes is in one figure and the other in another one.
  • this will work ok in noninteractive mode but for interactive mode (moving between components with arrows etc.) will not be that useful because it is generally faster/better to update properties of an existing axis than to clear it and fill up with content again. This is not a major problem - it just means some code repetition will be necessary for interaction
  • there is one thing I would prefer to change in plot_epochs_image - I think setting vmin and vmax should be performed after using ndimage.gaussian_filter1d because smoothing changes min, max limits
  • are you are ok with image_kwargs dict argumen to plot_ica_properties that allows to pass arguments to plot_epochs_image (plot_ica_properties(ica, epochs, image_kwargs={'sigma': 2.5}))? It is not very convenient but it is easy to handle. However - I think that at least sigma argumet should be made accessible directly without the need to use image_kwargs (plot_ica_properties(ica, epochs, sigma=2.5))

@mmagnuski
Copy link
Member Author

Now plot_properties uses plot_epochs_image, below are the changes:
before_using_plot_epochs_image after_using_plot_epochs_image
The image looks less happy, but more elegant :) (I hate how despite knowing all the dangers of jet colormap I still like it...), there is also no smoothing along y now (previously sigma was 2. by default).
I am thinking about making erp line more consistent with spectrum line and variance scatter - changing its color to black and adding gray fill_between (the fill_between could be added to plot_epochs_image if you prefer - plot_std by default set to False for example).
The last aesthetic change I think I will add is adding more "air" at the top and bottom of spectrum plot (but I will have to see how it interacts with dB kwarg).


def _segment_raw(raw, segment_length=1., verbose=True):
"""Divide continuous raw data into equal-sized
consecutive chunks.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

chunks -> epochs

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

don't forget this one

@agramfort
Copy link
Member

@mmagnuski let us know when travis is happy and you're done with docstrings etc.

@jona-sassenhagen
Copy link
Contributor

there is one thing I would prefer to change in plot_epochs_image - I think setting vmin and vmax should be performed after using ndimage.gaussian_filter1d because smoothing changes min, max limits

Probably a different PR?

are you are ok with image_kwargs dict argumen to plot_ica_properties that allows to pass arguments to plot_epochs_image (plot_ica_properties(ica, epochs, image_kwargs={'sigma': 2.5}))? It is not very convenient but it is easy to handle. However - I think that at least sigma argumet should be made accessible directly without the need to use image_kwargs (plot_ica_properties(ica, epochs, sigma=2.5))

I'd prefer to keep it all in the image_kwargs dict actually. I think the choice of sigma in plot_epochs_image is often unfortunate, but that should be handled there, not here.

Haven't gone over the code, but the results look good.

@jona-sassenhagen
Copy link
Contributor

Is 120 Hz a good default for the max freq? Wouldn't something like 80 or so be better to more clearly see alpha vs. beta vs. theta vs. no low freq peaks?

@mmagnuski
Copy link
Member Author

mmagnuski commented Jun 4, 2016

@jona-sassenhagen @Eric89GXL @agramfort
My thoughts on the API:
plot_properties (instead of plot_ica_properties) that takes inst (Epochs or Raw, just like now) and following kwargs:

  • ica argument to pass ica (if ica is not passed, plot channel properties - next PR)
  • picks - if more than one, plot multiple figures
  • interactive, when set to True:
    • throws error when matplotlib backend is set to inline?
    • adds buttons and interactivity (instead of plotting multiple figures (when picks is a list longer than one element) allows to move left/right across picks with buttons and keyboard)
    • does not return figure and axes but PropertiesExplorer object (I will leave this discussion for another PR)
  • dB controlling whether spectrum is plotted in dB
  • cmap to control topo and image color maps with one param
  • colorbar? to add colorbar (with plot_epochs_image, will require additional axis)
  • topo_kwargs, image_kwargs to pass additional arguments to plot_topo and plot_epochs_image
  • any other? plot_std?

Then Epochs and Raw would have plot_properties(...) method that calls plot_properties(inst, ...) and ICA plot_properties(inst, ...) method would call plot_properties(inst, ica=self, ...).
I am thinking about a case when somebody does ica.get_sources(inst).plot_properties() - is there info about component weights in resulting object?

Currently plot_ica_properties returns fig and axes (list of axes) but given that multiple figures should be drawn if picks is a list longer than one it doesn't make much sense to return axes - so I'm going to return only the fig which would be a figure or a list of figures (or a plotter object if interactive=True in a future PR...)

@jona-sassenhagen
Copy link
Contributor

picks - if more than one, plot multiple figures

I'm not 100% sure if this works without conflict, API wise - picks usually does something different. For example, maybe you want to only look at the MEG channels for your ICA, so you pick the MEG channels. Then you need to further supply the component to be plotted. Do you get what I mean?
I don't have a good idea for what to call it, though.

Rest I either agree, or am indifferent on.

@mmagnuski
Copy link
Member Author

@jona-sassenhagen Yeah, I understand, I was following the way plot_epochs_image handles picks.

@jona-sassenhagen
Copy link
Contributor

Should stay consistent with that somehow, too. I'm not sure.

I guess just following the plot_epochs_image approach is best.

@agramfort
Copy link
Member

Stick to ICA objects only for now otherwise PR will take much longer to be merged and don't be too creative ;)

@mmagnuski
Copy link
Member Author

there is one thing I would prefer to change in plot_epochs_image - I think setting vmin and vmax should be performed after using ndimage.gaussian_filter1d because smoothing changes min, max limits

Probably a different PR?

Yes, good idea.

I'd prefer to keep it all in the image_kwargs dict actually. I think the choice of sigma in plot_epochs_image is often unfortunate, but that should be handled there, not here.

I didn't think of setting the parameter in plot_properties - just exposing it for convenience and passing it to plot_epochs_image like:

if sigma is not None:
    image_kwargs.update(sigma=sigma)

I just expect sigma to be used quite frequently and image_kwargs={sigma: 1.5} isn't very handy.

Is 120 Hz a good default for the max freq? Wouldn't something like 80 or so be better to more clearly see alpha vs. beta vs. theta vs. no low freq peaks?

There is no default currently, all freqs are plotted. Maybe an additional argument?

Stick to ICA objects only for now otherwise PR will take much longer to be merged and don't be too creative ;)

Yes, I was planning to stick to ICA objs now - I was just trying to plan ahead and sketch my thoughts to see if you agree. I will try tro restrict my creativity ;)

@mmagnuski
Copy link
Member Author

An example of current functionality:
plot_properties_v0 2

is_correct_type = np.array([isinstance(x, mpl.axes._axes.Axes)
for x in axes])
if not np.all(is_correct_type):
first_bad = np.where(is_correct_type == False)[0][0]
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

flake doesn't like this :(
it doesn't notice it is numpy array and tells me:
comparison to False should be if cond is False: or if not cond:
Pretty weird, but I'll change to numpy.logical_not and remove the comparison

Copy link
Member

@larsoner larsoner Jun 4, 2016 via email

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mmagnuski
Copy link
Member Author

@agramfort
Rebased.
Sorry for the number of commits - what rules should I follow in the future so that you don't have to squash?

@agramfort agramfort merged commit 6f69b79 into mne-tools:master Jun 27, 2016
@agramfort
Copy link
Member

thanks @mmagnuski for the hard work and the patience

@espressofiend feel free to complain if it's not good enough for your needs.

@jona-sassenhagen
Copy link
Contributor

I love this, great work @mmagnuski.

jona-sassenhagen added a commit to jona-sassenhagen/mne-python that referenced this pull request Jun 27, 2016
* 'master' of git://github.com/mne-tools/mne-python:
  [MRG] plot_ica_properties (mne-tools#3275)
  [MRG] Example link (mne-tools#3346)
  GUI (coregistration):  scale step +/- instead of multiplicative (mne-tools#3345)
  Fix bug in set_bipolar_reference (mne-tools#3343)
  fixed spelling mistake (mne-tools#3341)
  Fixed spelling mistake (mne-tools#3339)
@mmagnuski
Copy link
Member Author

Thanks! And thanks for your helpful comments!
If you have any tips on how commits should be structured - let me know!

@mmagnuski mmagnuski deleted the plot_properties branch June 27, 2016 22:26
@dengemann
Copy link
Member

Thanks @mmagnuski for doing this! I did not have time to look during the last weeks unfortunately. This is really great work!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

7 participants