ENH: dict features for Epochs objects

examples/plot_topo_compare_conditions.py

@@ -0,0 +1,65 @@
+"""
+===============================================
+Compare Evoked Reponses to Different Conditions
+===============================================
+
+"""
+
+# Authors: Denis Engemann <d.engemann@fz-juelich.de>
+#          Alexandre Gramfort <gramfort@nmr.mgh.harvard.edu>
+#
+# License: BSD (3-clause)
+
+print __doc__
+
+import pylab as pl
+import mne
+
+from mne.fiff import Raw, pick_types
+from mne.layouts import read_layout
+from mne.viz import plot_topo
+from mne.datasets import sample
+data_path = sample.data_path('.')
+
+###############################################################################
+# Set parameters
+raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'
+event_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'
+event_id = 1
+tmin = -0.2
+tmax = 0.5
+
+#   Setup for reading the raw data
+raw = Raw(raw_fname)
+events = mne.read_events(event_fname)
+
+#   Set up pick list: EEG + STI 014 - bad channels (modify to your needs)
+include = []  # or stim channels ['STI 014']
+exclude = raw.info['bads'] + ['EEG 053']  # bads + 1 more
+
+# pick EEG channels
+picks = pick_types(raw.info, meg=True, eeg=False, stim=False, eog=True,
+                                            include=include, exclude=exclude)
+# Read epochs
+epochs = mne.Epochs(raw, events, dict(audio_l=1, visual_r=3), tmin, tmax, picks=picks,
+                  baseline=(None, 0), reject=dict(eog=150e-6))
+
+# access sub-epochs by conditions labels
+epochs_au, epochs_vi = epochs['audio_l'], epochs['visual_r']
+
+print epochs_au, epochs_au
+
+layout = read_layout('Vectorview-all')
+
+evoked_au, evoked_vi = epochs_au.average(), epochs_vi.average()
+
+
+###############################################################################
+# Show topography for two different conditions
+
+pl.close('all')
+for evoked in [evoked_au, evoked_vi]:
+    pl.figure()
+    title = 'MNE sample data (condition : %s)' % evoked.comment
+    plot_topo(evoked, layout, title=title)
+    pl.show()

mne/epochs.py

@@ -30,6 +30,7 @@
 from .filter import resample
 from .event import _read_events_fif
 from . import verbose
+from .fixes import in1d
 
 
 class Epochs(object):
@@ -43,8 +44,13 @@ class Epochs(object):
     events : array, of shape [n_events, 3]
         Returned by the read_events function
 
-    event_id : int | None
-        The id of the event to consider. If None all events are used.
+    event_id : int | dict | None
+        The id of the event to consider. If dict,
+        the keys can later be used to acces associated events. Example:
+        dict(auditory=1, visual=3). If int, a dict will be created with
+        the id as string. If None, all events will be used with
+        and a dict is created with string integer names corresponding
+        to the event id integers.
 
     tmin : float
         Start time before event
@@ -105,6 +111,9 @@ class Epochs(object):
     info: dict
         Measurement info
 
+    event_id : dict
+        Names of  of conditions corresponding to event_ids.
+
     ch_names: list of string
         List of channels' names
 
@@ -121,7 +130,7 @@ class Epochs(object):
     get_data() : self
         Return all epochs as a 3D array [n_epochs x n_channels x n_times].
 
-    average(picks=None) : self
+    average(picks=None) : Evoked
         Return Evoked object containing averaged epochs as a
         2D array [n_channels x n_times].
 
@@ -139,6 +148,14 @@ class Epochs(object):
     resample() : self, int, int, int, string or list
         Resample preloaded data.
 
+    as_data_frame() : DataFrame
+        Export Epochs object as Pandas DataFrame for subsequent statistical
+        analyses.
+
+    to_nitime() : TimeSeries
+        Export Epochs object as nitime TimeSeries object for subsequent
+        analyses.
+
     Notes
     -----
     For indexing and slicing:
@@ -148,6 +165,13 @@ class Epochs(object):
     epochs[idx] : Epochs
         Return Epochs object with a subset of epochs (supports single
         index and python style slicing)
+
+    For subset selection using categorial labels:
+
+    epochs['name'] : Epochs
+        Return Epochs object with a subset of epochs corresponding to an
+        experimental condition as specified by event_id.
+
     """
     @verbose
     def __init__(self, raw, events, event_id, tmin, tmax, baseline=(None, 0),
@@ -159,10 +183,21 @@ def __init__(self, raw, events, event_id, tmin, tmax, baseline=(None, 0),
 
         self.raw = raw
         self.verbose = raw.verbose if verbose is None else verbose
-        self.event_id = event_id
+        self.name = name
+        if isinstance(event_id, dict):
+            if not all([isinstance(v, int) for v in event_id.values()]):
+                raise ValueError('Event IDs must be of type integer')
+            if not all([isinstance(k, str) for k in event_id]):
+                raise ValueError('Event names must be of type str')
+            self.event_id = event_id
+        elif isinstance(event_id, int):
+            self.event_id = {str(event_id): event_id}
+        elif event_id is None:
+            self.event_id = dict((str(e), e) for e in np.unique(events[:, 2]))
+        else:
+            raise ValueError('event_id must be dict or int.')
         self.tmin = tmin
         self.tmax = tmax
-        self.name = name
         self.keep_comp = keep_comp
         self.dest_comp = dest_comp
         self.baseline = baseline
@@ -209,10 +244,9 @@ def __init__(self, raw, events, event_id, tmin, tmax, baseline=(None, 0),
 
         #    Select the desired events
         self.events = events
-        if event_id is not None:
-            selected = np.logical_and(events[:, 1] == 0,
-                                      events[:, 2] == event_id)
-            self.events = self.events[selected]
+        overlap = in1d(events[:, 2], self.event_id.values())
+        selected = np.logical_and(events[:, 1] == 0, overlap)
+        self.events = events[selected]
 
         n_events = len(self.events)
 
@@ -463,26 +497,47 @@ def __repr__(self):
     def __getitem__(self, key):
         """Return an Epochs object with a subset of epochs
         """
-        if not self._bad_dropped:
-            warnings.warn("Bad epochs have not been dropped, indexing will be "
-                          "inaccurate. Use drop_bad_epochs() or preload=True")
 
-        epochs = cp.copy(self)  # XXX : should use deepcopy but breaks ...
-        epochs.events = np.atleast_2d(self.events[key])
+        if self.preload:
+            data = self._data
+            del self._data
+
+        epochs = self.copy()
 
         if self.preload:
-            if isinstance(key, slice):
-                epochs._data = self._data[key]
-            else:
-                key = np.atleast_1d(key)
-                epochs._data = self._data[key]
+            self._data = data
+
+        if not isinstance(key, str):
+            if not self._bad_dropped:
+                warnings.warn("Bad epochs have not been dropped, indexing will"
+                              " be inaccurate. Use drop_bad_epochs() or"
+                              " preload=True")
+
+            epochs.events = np.atleast_2d(self.events[key])
+
+            if self.preload:
+                if isinstance(key, slice):
+                    epochs._data = self._data[key]
+                else:
+                    key = np.atleast_1d(key)
+                    epochs._data = self._data[key]
+        elif key not in self.event_id:
+            raise KeyError('Event "%s" is not in Epochs.' % key)
+        else:
+            epochs.events = self._get_events(key)
+            if self.preload:
+                epochs._data = epochs._data[epochs.events[:, 0]]
+            epochs.name = (key if epochs.name == 'Unknown'
+                           else 'epochs_%s' % key)
+
         return epochs
 
     def average(self, picks=None):
         """Compute average of epochs
 
         Parameters
         ----------
+
         picks: None | array of int
             If None only MEG and EEG channels are kept
             otherwise the channels indices in picks are kept.
@@ -492,6 +547,7 @@ def average(self, picks=None):
         evoked : Evoked instance
             The averaged epochs
         """
+
         return self._compute_mean_or_stderr(picks, 'ave')
 
     def standard_error(self, picks=None):
@@ -510,12 +566,18 @@ def standard_error(self, picks=None):
         """
         return self._compute_mean_or_stderr(picks, 'stderr')
 
+    def _get_events(self, event_name):
+        """Aux function: get event ids"""
+        ids = None
+        if event_name in self.event_id:
+            ids = self.events[self.events[0:, 2] == self.event_id[event_name]]
+
+        return ids
+
     def _compute_mean_or_stderr(self, picks, mode='ave'):
         """Compute the mean or std over epochs and return Evoked"""
-        if mode == 'stderr':
-            _do_std = True
-        else:
-            _do_std = False
+
+        _do_std = True if mode == 'stderr' else False
         evoked = Evoked(None)
         evoked.info = cp.deepcopy(self.info)
         n_channels = len(self.ch_names)
@@ -610,7 +672,7 @@ def crop(self, tmin=None, tmax=None, copy=False):
         tmask = (self.times >= tmin) & (self.times <= tmax)
         tidx = np.where(tmask)[0]
 
-        this_epochs = self if not copy else cp.deepcopy(self)
+        this_epochs = self if not copy else self.copy()
         this_epochs.tmin = this_epochs.times[tidx[0]]
         this_epochs.tmax = this_epochs.times[tidx[-1]]
         this_epochs.times = this_epochs.times[tmask]
@@ -643,9 +705,12 @@ def resample(self, sfreq, npad=100, window='boxcar'):
     def copy(self):
         """ Return copy of Epochs instance
         """
-        raw = self.raw.copy()
+        raw = self.raw
+        del self.raw
         new = deepcopy(self)
+        self.raw = raw
         new.raw = raw
+
         return new
 
     def save(self, fname):
@@ -704,7 +769,6 @@ def save(self, fname):
         end_block(fid, FIFF.FIFFB_MEAS)
         end_file(fid)
 
-
     def as_data_frame(self, frame=True):
         """Get the epochs as Pandas panel of data frames
 
@@ -1011,7 +1075,7 @@ def read_epochs(fname, proj=True, verbose=None):
     epochs._projector, epochs.info = setup_proj(info)
     epochs.ch_names = info['ch_names']
     epochs.baseline = baseline
-    epochs.event_id = int(np.unique(events[:, 2]))
+    epochs.event_id = dict((str(e), e) for e in np.unique(events[:, 2]))
     fid.close()
 
     return epochs
@@ -1038,7 +1102,7 @@ def bootstrap(epochs, random_state=None):
                            'in the constructor.')
 
     rng = check_random_state(random_state)
-    epochs_bootstrap = cp.deepcopy(epochs)
+    epochs_bootstrap = epochs.copy()
     n_events = len(epochs_bootstrap.events)
     idx = rng.randint(0, n_events, n_events)
     epochs_bootstrap = epochs_bootstrap[idx]

mne/tests/test_epochs.py

@@ -4,7 +4,7 @@
 # License: BSD (3-clause)
 
 import os.path as op
-from nose.tools import assert_true, assert_equal
+from nose.tools import assert_true, assert_equal, assert_raises
 from numpy.testing import assert_array_equal, assert_array_almost_equal
 import numpy as np
 import copy as cp
@@ -381,3 +381,19 @@ def test_epoch_eq():
     assert_true(epochs_1.events.shape[0] != epochs_2.events.shape[0])
     equalize_epoch_counts(epochs_1, epochs_2)
     assert_true(epochs_1.events.shape[0] == epochs_2.events.shape[0])
+
+
+def test_access_by_name():
+    assert_raises(ValueError, Epochs, raw, events, {1: 42, 2: 42}, tmin,
+                 tmax, picks=picks)
+    assert_raises(ValueError, Epochs, raw, events, {'a': 'spam', 2: 'eggs'},
+                 tmin, tmax, picks=picks)
+    assert_raises(ValueError, Epochs, raw, events, {'a': 'spam', 2: 'eggs'},
+                 tmin, tmax, picks=picks)
+    assert_raises(ValueError, Epochs, raw, events, 'foo', tmin, tmax, picks=picks)
+    epochs = Epochs(raw, events, {'a': 1, 'b': 2}, tmin, tmax, picks=picks)
+    assert_raises(KeyError, epochs.__getitem__, 'bar')
+
+    data = epochs['a'].get_data()
+    event_a = events[events[:, 2] == 1]
+    assert_true(len(data) == len(event_a))