Skip to content

Conversation

@agramfort
Copy link
Member

closes #4736 and #4526

@jona-sassenhagen @ktavabi I took over both of your PRs and simplify to refactor / simplify

please have a look

ylim=dict(), invert_y=False, show_sensors=None,
show_legend=True, axes=None, title=None, show=True):
def _combine_grad(evoked, picks):
"""Creates a new instance of Evoked with combined gradiometers (RMSE)"""
Copy link
Contributor

Choose a reason for hiding this comment

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

Uuuuh, that's good to have.



def _check_loc_legal(loc, what='your choice'):
"""Check if a loc is a legal for MPL."""
Copy link
Contributor

Choose a reason for hiding this comment

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

"""Check if loc is a legal location for MPL subordinate axes."""
(or something like that, it's missing a noun)



def _format_evokeds_colors(evokeds, cmap, colors):
# set up labels and instances to have evokeds as a dict
Copy link
Contributor

Choose a reason for hiding this comment

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

Sorry, what? :) I don't get what this sentence means.

channels using global field power (GFP) computation, else it is taking
a plain mean.
This function is useful for comparing ER[P/F]s at a specific location.
Copy link
Contributor

Choose a reason for hiding this comment

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

This function is useful for comparing multiple ER[P/F]s - e.g., for multiple conditions - at a specific location.

It can plot a simple Evoked object or if supplied with a list or a dict
of lists of evoked instances, it will compute the average of the mean or
GFP across datasets. It can also compute some confidence intervals
using for example bootstrap estimates.
Copy link
Contributor

Choose a reason for hiding this comment

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

It plots 1. a simple evoked object, 2. a list or dict of evoked objects (e.g., for multiple conditions), 3. a list or dict of lists of evokeds (e.g., for multiple subjects in multiple conditions). In the last case, it can show a confidence interval (across e.g. subjects) using parametric or bootstrap estimation.

show_legend : bool | int
If not False, show a legend. If int, the position of the axes
show_legend : bool | str | int
If not False, show a legend. If int or str, the position of the axes
Copy link
Contributor

Choose a reason for hiding this comment

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

the position of the legend
or
the position of the legend axes

show_sensors: bool | int | str | None
If not False, channel locations are plotted on a small head circle.
If an int, the position of the axes (forwarded to
If int or str, the position of the axes (forwarded to
Copy link
Contributor

Choose a reason for hiding this comment

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

The position of the head circle
or
The position of the head circle axes

# dealing with continuous colors
the_colors = None
color_conds = None
color_order = None
Copy link
Contributor

Choose a reason for hiding this comment

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

the_colors, color_conds, color_order = None, None, None

mne/evoked.py Outdated
def _check_evokeds_ch_names_times(all_evoked):
evoked = all_evoked[0]
ch_names = evoked.ch_names
for e in all_evoked[1:]:
Copy link
Contributor

Choose a reason for hiding this comment

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

single-char var name ...

metadata = epochs.metadata
is_concrete = metadata["Concreteness"] > metadata["Concreteness"].median()
metadata["is_concrete"] = np.where(is_concrete, 'Concrete', 'Abstract')
is_concrete = metadata["NumberOfLetters"] > 5
Copy link
Contributor

Choose a reason for hiding this comment

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

is_long

@jona-sassenhagen
Copy link
Contributor

Thanks Alex, that looks very clean.

gradiometers are combined with RMSE. It means for example that for
a vectorview system with 204 gradiometers it will be transformed to
a 102 channels helmet.
Copy link
Contributor

Choose a reason for hiding this comment

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

For example/E.g., data from a Vectorview system with 204 gradiometers will be transformed to 102 channels.

gradiometer corresponding pairs with be combined.
If None, it defaults to all data channels, in which case the global
field power will be plotted for all channel type available.
gfp : bool
Copy link
Contributor

@ktavabi ktavabi Nov 15, 2017

Choose a reason for hiding this comment

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

If the selected channels are gradiometers, the signal from corresponding (gradiometer) pairs will be combined. If None, it defaults to all data channels, in which case the global field power will be plotted for all available channel types.

vlines : "auto" | list of float
A list in seconds at which to plot dashed vertical lines.
If "auto" and 0. is in the time interval of interest, it is set to [0.]
so a vertical bar is plotted at time 0.
Copy link
Contributor

@ktavabi ktavabi Nov 15, 2017

Choose a reason for hiding this comment

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

If "auto" and 0. ms is the time point of interest, it is set to [0.] and a vertical bar is plotted at time 0.

fig : Figure | list of Figures
The figure(s) in which the plot is drawn.
The figure(s) in which the plot is drawn. A list of figures
is returned only one multiple channel types are plotted.
Copy link
Contributor

Choose a reason for hiding this comment

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

When plotting multiple channel types, a list of figures, one for each channel type is returned.

@ktavabi
Copy link
Contributor

ktavabi commented Nov 15, 2017

Thanks @agramfort

When multiple channels are passed, this function combines them all, to
get one time course for each condition. If gfp is True it combines
channels using global field power (GFP) computation, else it is taking
a plain mean.
Copy link
Member

Choose a reason for hiding this comment

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

this preamble is 3 paragraphs, can we put this into Notes instead?

GFP across datasets. It can also compute some confidence intervals
using for example bootstrap estimates.
When passed with more than one planar gradiometers, the planar
Copy link
Member

Choose a reason for hiding this comment

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

When ``picks`` includes more than one planar gradiometer, ...

When passed with more than one planar gradiometers, the planar
gradiometers are combined with RMSE. It means for example that for
a vectorview system with 204 gradiometers it will be transformed to
a 102 channels helmet.
Copy link
Member

Choose a reason for hiding this comment

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

102-channel helmet

# The plot is styled with dictionary arguments, again using "/"-separated tags.
# We plot a MEG channel with a strong auditory response.
#
# For move advanced plotting using :func:`mne.viz.plot_compare_evokeds`
Copy link
Member

Choose a reason for hiding this comment

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

. See also

This function is useful for comparing ER[P/F]s at a specific location.
It can plot a simple Evoked object or if supplied with a list or a dict
Copy link
Member

Choose a reason for hiding this comment

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

:class:`mne.Evoked`

@larsoner
Copy link
Member

Seems to fix the problem for me. It would be good to add this test (I came up with this to help @ktavabi test his PR):

import numpy as np
from numpy.testing import assert_allclose
import matplotlib.pyplot as plt
import os.path as op
import mne

data_path = mne.datasets.testing.data_path()
fname_evoked = op.join(data_path, 'MEG', 'sample', 'sample_audvis-ave.fif')

rng = np.random.RandomState(0)
evoked = mne.read_evokeds(fname_evoked)[0]
evokeds = [evoked.copy() for _ in range(10)]
for evoked in evokeds:
    evoked.data += (rng.randn(*evoked.data.shape) *
                    np.std(evoked.data, axis=-1, keepdims=True))
for picks in ([0], [1], [2], [0, 2], [1, 2], [0, 1, 2],):
    figs = mne.viz.plot_compare_evokeds([evokeds], picks=picks, ci=0.95)
    if not isinstance(figs, list):
        figs = [figs]
    for fig in figs:
        ext = fig.axes[0].collections[0].get_paths()[0].get_extents()
        xs, ylim = ext.get_points().T
        assert_allclose(xs, evoked.times[[0, -1]])
        line = fig.axes[0].lines[0]
        xs = line.get_xdata()
        assert_allclose(xs, evoked.times)
        ys = line.get_ydata()
        assert (ys < ylim[1]).all()
        assert (ys > ylim[0]).all()
    plt.close('all')

@jona-sassenhagen
Copy link
Contributor

@agramfort want me to directly commit to the branch?

@codecov
Copy link

codecov bot commented Nov 17, 2017

Codecov Report

Merging #4755 into master will increase coverage by 0.01%.
The diff coverage is 92.83%.

@@            Coverage Diff             @@
##           master    #4755      +/-   ##
==========================================
+ Coverage   87.86%   87.87%   +0.01%     
==========================================
  Files         349      349              
  Lines       65819    66079     +260     
  Branches    11319    11383      +64     
==========================================
+ Hits        57829    58070     +241     
- Misses       5117     5130      +13     
- Partials     2873     2879       +6

@agramfort
Copy link
Member Author

agramfort commented Nov 17, 2017 via email

@jona-sassenhagen
Copy link
Contributor

Great, I'll try to get a bit of work done!

@jona-sassenhagen
Copy link
Contributor

Pushed a few lines directly to your branch, I will continue tomorrow.

@jona-sassenhagen
Copy link
Contributor

@larsoner I added your test.

@jona-sassenhagen
Copy link
Contributor

Green ...

raise ValueError('If evokeds is a dict and a cmap is passed, '
'you must specify the colors.')
# XXX : I am a bit concerned about the duplication of
# the colors and cmap parameters.
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 what do you think about this?
we don't use colors anywhere in doc and examples and it has a non-trivial interaction with cmap.

thoughts?

Copy link
Contributor

Choose a reason for hiding this comment

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

Hm? Maybe I misunderstand, but colors is all over the example.

colors specifies the mapping from conditions to line colors relatively to the color map, cmap specifies the colormap (and the name for the colorbar). I see two alternatives, both of which we ruled out:

  • use magic to read relative line color from the condition names (my original idea, which you ruled out)
  • demand colors values to be RGB or RGBA (possible, but requires more user work)

Essentially, option two is viable, but would require the user do something like

import seaborn as sns
from sklearn.preprocessing import QuantileTransformer

n_bins = 10
qt = QuantileTransformer(n_bins).fit(epochs.metadata[cond])
colorscale = sns.color_palette("viridis", n_bins)

colors = dict()
for ... in ...:
    colors[xxx] = colorscale[qt.transform(...)]

every time they run the function.

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 ok. forget it.

Copy link
Member

Choose a reason for hiding this comment

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

should the comment be removed, then?

@agramfort agramfort changed the title Fix compare evokeds [MRG] Fix compare evokeds Nov 20, 2017
@agramfort
Copy link
Member Author

good to go from my end.

Copy link
Member

@larsoner larsoner left a comment

Choose a reason for hiding this comment

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

Otherwise LGTM

mne/evoked.py Outdated
ch_names = evoked.ch_names
for ev in all_evoked[1:]:
assert ev.ch_names == ch_names, ValueError(
"%s and %s do not contain " "the same channels" % (evoked, ev))
Copy link
Member

Choose a reason for hiding this comment

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

assert statements are skipped with python -O, if this is meant to actually prevent some user error (rather than internal consistency check) this needs to be done the standard way

from numbers import Integral

import numpy as np
import matplotlib.lines as mlines
Copy link
Member

Choose a reason for hiding this comment

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

it probably doesn't matter for the backend choosing because you didn't import pyplot, but let's nest this one, too, to keep mne import times down

Copy link
Contributor

Choose a reason for hiding this comment

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

Ah I see, you told me to put it at the top of the function, not the file. So that's what I'll do.

raise ValueError('If evokeds is a dict and a cmap is passed, '
'you must specify the colors.')
# XXX : I am a bit concerned about the duplication of
# the colors and cmap parameters.
Copy link
Member

Choose a reason for hiding this comment

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

should the comment be removed, then?

- a list or dict of :class:`mne.Evoked` objects (e.g., for multiple
conditions),
- a list or dict of lists of :class:`mne.Evoked` (e.g., for multiple
subjects in multiple conditions).
Copy link
Member

Choose a reason for hiding this comment

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

no need to indent, doing so turns it from a standard bullet list to a quoted bullet list

if (tmin >= 0 or tmax <= 0) and vlines == [0.]:
vlines = list()

if vlines is "auto" and (tmin < 0 and tmax > 0):
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 use is with string comparisons


if vlines is "auto" and (tmin < 0 and tmax > 0):
vlines = [0.]
assert isinstance(vlines, list)
Copy link
Member

Choose a reason for hiding this comment

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

From the code above I don't think you've validated that it should be a list at this point, so this should be a full if not isinstance(vlines, (list, tuple)): raise TypeError(...)

if not isinstance(axes, list):
axes = [axes]
from .utils import _validate_if_list_of_axes
_validate_if_list_of_axes(axes, obligatory_len=len(ch_types))
Copy link
Member

Choose a reason for hiding this comment

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

why not axes = _validate_if_list_of_axes(...) and do the if not isinstance(...) in that function? It seems like part of the validation

Copy link
Contributor

Choose a reason for hiding this comment

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

I've decided against making that change here cause the function is somewhat widely used, mostly in contexts I don't know. I can do it in a separate PR if you still think it should be done (assign me if you want).

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 get it in a subsequent PR if we find time

# XXX: could possibly be refactored; plot_joint is doing a similar thing
if any([type_ not in _VALID_CHANNEL_TYPES for type_ in ch_types]):
raise ValueError("Non-data channel picked.")
axes = []
Copy link
Member

Choose a reason for hiding this comment

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

Could do instead:

axes = [plt.subplots(figsize=(8, 6))[1] for _ in range(len(ch_types))]

# calculate the CI
ci_dict = dict()
data_dict = dict()
for cond, this_evokeds in evokeds.items():
Copy link
Member

Choose a reason for hiding this comment

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

order here is not deterministic, probably better to do for cond in sorted(evokeds.keys()): this_evokeds = evokeds[key] ...


warnings.simplefilter('always') # enable b/c these tests throw warnings

rng = np.random.RandomState(0)
Copy link
Member

Choose a reason for hiding this comment

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

better to nest this in the function that uses it, lest the next contributor add a function and think it's okay to use the same rng (it could break your test with a heisenbug)

@jona-sassenhagen
Copy link
Contributor

jona-sassenhagen commented Nov 20, 2017 via email

@larsoner
Copy link
Member

I suppose that means we've asymptoted to a good result :)

Change it or not, then, up to you

@agramfort
Copy link
Member Author

@jona-sassenhagen are you on it?

@jona-sassenhagen
Copy link
Contributor

Yes, I will try to get it done today.

@jona-sassenhagen
Copy link
Contributor

Green.

@agramfort
Copy link
Member Author

@larsoner merge if you're happy

Copy link
Member

@larsoner larsoner left a comment

Choose a reason for hiding this comment

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

One small remaining problem then +1 for merge

mne/evoked.py Outdated
evoked = all_evoked[0]
ch_names = evoked.ch_names
for ev in all_evoked[1:]:
if not ev.ch_names == ch_names:
Copy link
Member

Choose a reason for hiding this comment

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

more easily expressed as ev.ch_names != ch_names but you can keep it this way

if vlines == "auto" and (tmin < 0 and tmax > 0):
vlines = [0.]
if not isinstance(vlines, (list, tuple)):
raise TypeError("vlines must be a list or tuple, not " + type(vlines))
Copy link
Member

Choose a reason for hiding this comment

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

This must not be covered by a test because it will error, you should use % as we usually do:

In [1]: 'a' + type('b')
Traceback (most recent call last):

  File "<ipython-input-1-03514d031c8f>", line 1, in <module>
    'a' + type('b')

TypeError: must be str, not type

if axes is not None:
if not isinstance(axes, list):
axes = [axes]
from .utils import _validate_if_list_of_axes
Copy link
Member

Choose a reason for hiding this comment

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

this shouldn't need to be nested, in fact there is already a from .utils import ... line at the top

@jona-sassenhagen
Copy link
Contributor

@ktavabi you are happy too?

@ktavabi
Copy link
Contributor

ktavabi commented Nov 21, 2017

I've been happy since the takeover. Thanks. Just waiting to get back on master :)

@larsoner larsoner merged commit 44bac44 into mne-tools:master Nov 22, 2017
@larsoner
Copy link
Member

Thanks @agramfort @ktavabi @jona-sassenhagen

@jona-sassenhagen
Copy link
Contributor

Thanks guys!

fraimondo pushed a commit to fraimondo/mne-python that referenced this pull request Nov 29, 2017
* [ENH] plot_compare_evoked improved line colouring

* FIX plot_compared_evokeds grad ci estimation

* Compute stats after computing RMS

* reorg + simplify grad handling of plot_compare_evokeds

* use auxilary functions

* add exception and test for corner case with both cmap and colors

* address comments

* cosmit, docstrings

* pep8

* add erics test for grad cis

* docstring

* comments

* pep8 :|

* comments
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.

4 participants