[med-svn] [python-mne] 176/376: API : changing epoch rejection API

Yaroslav Halchenko debian at onerussian.com
Fri Nov 27 17:22:31 UTC 2015


This is an automated email from the git hooks/post-receive script.

yoh pushed a commit to annotated tag v0.1
in repository python-mne.

commit 0fd67ce84b3a2cf2fc256786957fe0104ea01eed
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date:   Wed Mar 30 23:08:39 2011 -0400

    API : changing epoch rejection API
---
 examples/plot_from_raw_to_epochs_to_evoked.py      |   5 +-
 examples/plot_read_epochs.py                       |   4 +-
 examples/plot_whitened_evoked_data.py              |   5 +-
 .../plot_cluster_1samp_test_time_frequency.py      |   5 +-
 examples/stats/plot_cluster_stats_evoked.py        |  11 +-
 .../stats/plot_cluster_stats_time_frequency.py     |  13 +-
 examples/stats/plot_sensor_permutation_test.py     |   5 +-
 examples/time_frequency/plot_time_frequency.py     |   5 +-
 mne/epochs.py                                      | 172 ++++++++++-----------
 mne/misc.py                                        |  31 +++-
 mne/tests/test_epochs.py                           |  12 +-
 11 files changed, 143 insertions(+), 125 deletions(-)

diff --git a/examples/plot_from_raw_to_epochs_to_evoked.py b/examples/plot_from_raw_to_epochs_to_evoked.py
index ac77f86..5914cfa 100644
--- a/examples/plot_from_raw_to_epochs_to_evoked.py
+++ b/examples/plot_from_raw_to_epochs_to_evoked.py
@@ -39,9 +39,8 @@ exclude = raw.info['bads'] + ['EEG 053'] # bads + 1 more
 picks = fiff.pick_types(raw.info, meg=False, eeg=True, stim=False, eog=True,
                                             include=include, exclude=exclude)
 # Read epochs
-epochs = mne.Epochs(raw, events, event_id,
-                            tmin, tmax, picks=picks, baseline=(None, 0))
-epochs.reject(eeg=40e-6, eog=150e-6)
+epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
+                    baseline=(None, 0), reject=dict(eeg=40e-6, eog=150e-6))
 evoked = epochs.average() # average epochs and get an Evoked dataset.
 
 evoked.save('sample_audvis_eeg-ave.fif') # save evoked data to disk
diff --git a/examples/plot_read_epochs.py b/examples/plot_read_epochs.py
index 6f30c51..feeed3b 100644
--- a/examples/plot_read_epochs.py
+++ b/examples/plot_read_epochs.py
@@ -38,8 +38,8 @@ picks = fiff.pick_types(raw.info, meg=True, eeg=False, stim=True, eog=True,
 
 # Read epochs
 epochs = mne.Epochs(raw, events, event_id, tmin, tmax,
-                    picks=picks, baseline=(None, 0), preload=True)
-epochs.reject(grad=4000e-13, mag=4e-12, eog=150e-6)
+                    picks=picks, baseline=(None, 0), preload=True,
+                    reject=dict(grad=4000e-13, mag=4e-12, eog=150e-6))
 evoked = epochs.average() # average epochs to get the evoked response
 
 ###############################################################################
diff --git a/examples/plot_whitened_evoked_data.py b/examples/plot_whitened_evoked_data.py
index 890df5c..c1e5db4 100644
--- a/examples/plot_whitened_evoked_data.py
+++ b/examples/plot_whitened_evoked_data.py
@@ -37,9 +37,8 @@ events = mne.find_events(raw)
 exclude = raw.info['bads'] + ['EEG 053'] # bads + 1 more
 picks = fiff.pick_types(raw.info, meg=False, eeg=True, stim=False, eog=True,
                         exclude=exclude)
-epochs = mne.Epochs(raw, events, event_id,
-                    tmin, tmax, picks=picks, baseline=(None, 0))
-epochs.reject(eeg=40e-6, eog=150e-6)
+epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
+                    baseline=(None, 0), reject=dict(eeg=40e-6, eog=150e-6))
 evoked = epochs.average() # average epochs and get an Evoked dataset.
 
 cov = mne.Covariance(cov_fname)
diff --git a/examples/stats/plot_cluster_1samp_test_time_frequency.py b/examples/stats/plot_cluster_1samp_test_time_frequency.py
index 29312fd..8290ffa 100644
--- a/examples/stats/plot_cluster_1samp_test_time_frequency.py
+++ b/examples/stats/plot_cluster_1samp_test_time_frequency.py
@@ -51,9 +51,8 @@ picks = fiff.pick_types(raw.info, meg='grad', eeg=False, eog=True,
 
 # Load condition 1
 event_id = 1
-epochs = mne.Epochs(raw, events, event_id,
-                    tmin, tmax, picks=picks, baseline=(None, 0))
-epochs.reject(grad=4000e-13, eog=150e-6)
+epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
+                    baseline=(None, 0), reject(grad=4000e-13, eog=150e-6))
 data = epochs.get_data() # as 3D matrix
 data *= 1e13 # change unit to fT / cm
 # Time vector
diff --git a/examples/stats/plot_cluster_stats_evoked.py b/examples/stats/plot_cluster_stats_evoked.py
index 05bf10d..5ab6d02 100644
--- a/examples/stats/plot_cluster_stats_evoked.py
+++ b/examples/stats/plot_cluster_stats_evoked.py
@@ -40,15 +40,14 @@ include = [channel]
 # Read epochs for the channel of interest
 picks = fiff.pick_types(raw.info, meg=False, eog=True, include=include)
 event_id = 1
-epochs1 = mne.Epochs(raw, events, event_id,
-                            tmin, tmax, picks=picks, baseline=(None, 0))
-epochs1.reject(grad=4000e-13, eog=150e-6)
+reject = dict(grad=4000e-13, eog=150e-6)
+epochs1 = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
+                     baseline=(None, 0), reject=reject)
 condition1 = epochs1.get_data() # as 3D matrix
 
 event_id = 2
-epochs2 = mne.Epochs(raw, events, event_id,
-                            tmin, tmax, picks=picks, baseline=(None, 0))
-epochs2.reject(grad=4000e-13, eog=150e-6)
+epochs2 = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
+                     baseline=(None, 0), reject=reject)
 condition2 = epochs2.get_data() # as 3D matrix
 
 condition1 = condition1[:,0,:] # take only one channel to get a 2D array
diff --git a/examples/stats/plot_cluster_stats_time_frequency.py b/examples/stats/plot_cluster_stats_time_frequency.py
index ac3eb05..bf0392d 100644
--- a/examples/stats/plot_cluster_stats_time_frequency.py
+++ b/examples/stats/plot_cluster_stats_time_frequency.py
@@ -53,18 +53,19 @@ picks = fiff.pick_types(raw.info, meg='grad', eeg=False, eog=True,
 ch_name = raw.info['ch_names'][picks[0]]
 
 # Load condition 1
+reject = dict(grad=4000e-13, eog=150e-6)
 event_id = 1
-epochs_condition_1 = mne.Epochs(raw, events, event_id,
-                    tmin, tmax, picks=picks, baseline=(None, 0))
-epochs_condition_1.reject(grad=4000e-13, eog=150e-6)
+epochs_condition_1 = mne.Epochs(raw, events, event_id, tmin, tmax,
+                                picks=picks, baseline=(None, 0),
+                                reject=reject)
 data_condition_1 = epochs_condition_1.get_data() # as 3D matrix
 data_condition_1 *= 1e13 # change unit to fT / cm
 
 # Load condition 2
 event_id = 2
-epochs_condition_2 = mne.Epochs(raw, events, event_id,
-                    tmin, tmax, picks=picks, baseline=(None, 0))
-epochs_condition_2.reject(grad=4000e-13, eog=150e-6)
+epochs_condition_2 = mne.Epochs(raw, events, event_id, tmin, tmax,
+                                picks=picks, baseline=(None, 0),
+                                reject=reject)
 data_condition_2 = epochs_condition_2.get_data() # as 3D matrix
 data_condition_2 *= 1e13 # change unit to fT / cm
 
diff --git a/examples/stats/plot_sensor_permutation_test.py b/examples/stats/plot_sensor_permutation_test.py
index 7dd8e9c..c219469 100644
--- a/examples/stats/plot_sensor_permutation_test.py
+++ b/examples/stats/plot_sensor_permutation_test.py
@@ -42,9 +42,8 @@ exclude = raw.info['bads'] + ['MEG 2443', 'EEG 053'] # bads + 2 more
 # pick MEG Magnetometers
 picks = fiff.pick_types(raw.info, meg='grad', eeg=False, stim=False, eog=True,
                                             include=include, exclude=exclude)
-epochs = mne.Epochs(raw, events, event_id,
-                            tmin, tmax, picks=picks, baseline=(None, 0))
-epochs.reject(grad=4000e-13, eog=150e-6)
+epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
+                    baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6))
 data = epochs.get_data()
 times = epochs.times
 
diff --git a/examples/time_frequency/plot_time_frequency.py b/examples/time_frequency/plot_time_frequency.py
index 7133d52..df35773 100644
--- a/examples/time_frequency/plot_time_frequency.py
+++ b/examples/time_frequency/plot_time_frequency.py
@@ -41,9 +41,8 @@ exclude = raw.info['bads'] + ['MEG 2443', 'EEG 053'] # bads + 2 more
 picks = fiff.pick_types(raw.info, meg='grad', eeg=False, eog=True,
                                 stim=False, include=include, exclude=exclude)
 
-epochs = mne.Epochs(raw, events, event_id,
-                    tmin, tmax, picks=picks, baseline=(None, 0))
-epochs.reject(grad=4000e-13, eog=150e-6)
+epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
+                    baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6))
 data = epochs.get_data() # as 3D matrix
 evoked = epochs.average() # compute evoked fields
 
diff --git a/mne/epochs.py b/mne/epochs.py
index ea65192..b149396 100644
--- a/mne/epochs.py
+++ b/mne/epochs.py
@@ -36,9 +36,6 @@ class Epochs(object):
     keep_comp : boolean
         Apply CTF gradient compensation
 
-    info : dict
-        Measurement info
-
     baseline: None (default) or tuple of length 2
         The time interval to apply baseline correction.
         If None do not apply it. If baseline is (a, b)
@@ -53,6 +50,22 @@ class Epochs(object):
         or wait before accessing each epoch (more memory
         efficient but can be slower).
 
+    reject : dict
+        Epoch rejection parameters based on peak to peak amplitude.
+        Valid keys are 'grad' | 'mag' | 'eeg' | 'eog' | 'ecg'.
+        If reject is None then no rejection is done.
+        Values are float. Example:
+        reject = dict(grad=4000e-13, # T / m (gradiometers)
+                      mag=4e-12, # T (magnetometers)
+                      eeg=40e-6, # uV (EEG channels)
+                      eog=250e-6 # uV (EOG channels)
+                      )
+
+    flat : dict
+        Epoch rejection parameters based on flatness of signal
+        Valid keys are 'grad' | 'mag' | 'eeg' | 'eog' | 'ecg'
+        If flat is None then no rejection is done.
+
     Methods
     -------
     get_epoch(i) : self
@@ -69,7 +82,7 @@ class Epochs(object):
 
     def __init__(self, raw, events, event_id, tmin, tmax, baseline=(None, 0),
                 picks=None, name='Unknown', keep_comp=False, dest_comp=0,
-                preload=False):
+                preload=False, reject=None, flat=None):
         self.raw = raw
         self.event_id = event_id
         self.tmin = tmin
@@ -80,6 +93,8 @@ class Epochs(object):
         self.dest_comp = dest_comp
         self.baseline = baseline
         self.preload = preload
+        self.reject = reject
+        self.flat = flat
 
         # Handle measurement info
         self.info = copy.copy(raw.info)
@@ -147,12 +162,12 @@ class Epochs(object):
         self.times = np.arange(int(round(tmin*sfreq)), int(round(tmax*sfreq)),
                           dtype=np.float) / sfreq
 
+        # setup epoch rejection
+        self._reject_setup()
+
         if self.preload:
             self._data = self._get_data_from_disk()
 
-    def __len__(self):
-        return len(self.events)
-
     def get_epoch(self, idx):
         """Load one epoch
 
@@ -204,9 +219,52 @@ class Epochs(object):
         n_times = len(self.times)
         n_events = len(self.events)
         data = np.empty((n_events, n_channels, n_times))
+        cnt = 0
+        n_reject = 0
         for k in range(n_events):
-            data[k] = self._get_epoch_from_disk(k)
-        return data
+            e = self._get_epoch_from_disk(k)
+            if self._is_good_epoch(e):
+                data[cnt] = self._get_epoch_from_disk(k)
+                cnt += 1
+            else:
+                n_reject += 1
+        print "Rejecting %d epochs." % n_reject
+        return data[:cnt]
+
+    def _is_good_epoch(self, e):
+        """Test is epoch e is good
+        """
+        if self.reject is not None:
+            for key, thresh in self.reject.iteritems():
+                idx = self._channel_type_idx[key]
+                name = key.upper()
+                if len(idx) > 0:
+                    e_idx = e[idx]
+                    deltas = np.max(e_idx, axis=1) - np.min(e_idx, axis=1)
+                    idx_max_delta = np.argmax(deltas)
+                    delta = deltas[idx_max_delta]
+                    if delta > thresh:
+                        ch_name = self.ch_names[idx[idx_max_delta]]
+                        print '\tRejecting epoch based on %s : %s (%s > %s).' \
+                                    % (name, ch_name, delta, thresh)
+                        return False
+        if self.flat is not None:
+            for key, thresh in self.flat.iteritems():
+                idx = self._channel_type_idx[key]
+                name = key.upper()
+                if len(idx) > 0:
+                    e_idx = e[idx]
+                    deltas = np.max(e_idx, axis=1) - np.min(e_idx, axis=1)
+                    idx_min_delta = np.argmin(deltas)
+                    delta = deltas[idx_min_delta]
+                    if delta < thresh:
+                        ch_name = self.ch_names[idx[idx_min_delta]]
+                        print '\tRejecting flat epoch based on %s : %s (%s < %s).' \
+                                    % (name, ch_name, delta, thresh)
+                        return False
+
+        return True
+
 
     def get_data(self):
         """Get all epochs as a 3D array
@@ -221,84 +279,26 @@ class Epochs(object):
         else:
             return self._get_data_from_disk()
 
-    def reject(self, grad=None, mag=None, eeg=None, eog=None):
-        """Reject some epochs based on threshold values
-
-        Parameters
-        ----------
-        grad : float
-            Max value for gradiometers. (about 4000 fT / cm = 4000e-13 T / m).
-            If None do not reject based on gradiometers.
-        mag : float
-            Max value for magnetometers. (about  4000 fT = 4e-12 T)
-            If None do not reject based on magnetometers.
-        eeg : float
-            Max value for EEG. (about 40e-6 uV)
-            If None do not reject based on EEG.
-        eog : float
-            Max value for EEG. (about 250e-6 uV)
-            If None do not reject based on EOG.
-
-        Returns
-        -------
-        data : array of shape [n_epochs, n_channels, n_times]
-            The epochs data
+    def _reject_setup(self):
+        """Setup reject process
         """
-        grad_idx = []
-        mag_idx = []
-        eeg_idx = []
-        eog_idx = []
-        for idx, ch in enumerate(self.ch_names):
-            if grad is not None and  channel_type(self.info, idx) == 'grad':
-                grad_idx.append(idx)
-            if mag is not None and channel_type(self.info, idx) == 'mag':
-                mag_idx.append(idx)
-            if eeg is not None and channel_type(self.info, idx) == 'eeg':
-                eeg_idx.append(idx)
-            if eog is not None and channel_type(self.info, idx) == 'eog':
-                eog_idx.append(idx)
-
-        if grad is not None and len(grad_idx) == 0:
-            raise ValueError("No GRAD channel found. Cannot reject based on GRAD.")
-        elif grad is not None:
-            print "grad reject : %s fT/cm" % (1e13*grad)
-        if mag is not None and len(mag_idx) == 0:
-            raise ValueError("No MAG channel found. Cannot reject based on MAG.")
-        elif mag is not None:
-            print "mag  reject : %s fT" % (1e15*mag)
-        if eeg is not None and len(eeg_idx) == 0:
-            raise ValueError("No EEG channel found. Cannot reject based on EEG.")
-        elif eeg is not None:
-            print "EEG  reject : %s uV" % (1e6*eeg)
-        if eog is not None and len(eog_idx) == 0:
-            raise ValueError("No EOG channel found. Cannot reject based on EOG.")
-        elif eog is not None:
-            print "EOG  reject : %s uV" % (1e6*eog)
-
-        good_epochs = []
-        for k, e in enumerate(self):
-            for thresh, idx, name in zip([grad, mag, eeg, eog],
-                                         [grad_idx, mag_idx, eeg_idx, eog_idx],
-                                         ['GRAD', 'MAG', 'EEG', 'EOG']):
-                if len(idx) > 0:
-                    e_idx = e[idx]
-                    deltas = np.max(e_idx, axis=1) - np.min(e_idx, axis=1)
-                    idx_max_delta = np.argmax(deltas)
-                    delta = deltas[idx_max_delta]
-                    if delta > thresh:
-                        ch_name = self.ch_names[idx[idx_max_delta]]
-                        print '\tRejecting epoch based on %s : %s (%s > %s).' \
-                                    % (name, ch_name, delta, thresh)
-                        break
-            else:
-                good_epochs.append(k)
-
-        n_good_epochs = len(good_epochs)
-        print "Keeping %d epochs (%d bad)" % (n_good_epochs,
-                                              len(self.events) - n_good_epochs)
-        self.events = self.events[good_epochs]
-        if self.preload:
-            self._data = self._data[good_epochs]
+        if self.reject is None and self.flat is None:
+            return
+
+        idx = dict(grad=[], mag=[], eeg=[], eog=[], ecg=[])
+        for k, ch in enumerate(self.ch_names):
+            for key in idx.keys():
+                if channel_type(self.info, k) == key:
+                    idx[key].append(k)
+
+        for key in idx.keys():
+            if (self.reject is not None and key in self.reject) \
+                    or (self.flat is not None and key in self.flat):
+                if len(idx[key]) == 0:
+                    raise ValueError("No %s channel found. Cannot reject based"
+                                 " on %s." % (key.upper(), key.upper()))
+
+        self._channel_type_idx = idx
 
     def __iter__(self):
         """To iteration over epochs easy.
diff --git a/mne/misc.py b/mne/misc.py
index 271776c..1ca6a2f 100644
--- a/mne/misc.py
+++ b/mne/misc.py
@@ -64,14 +64,35 @@ def read_reject_parameters(fname):
     except:
         print("Error while reading %s" % fname)
 
-    reject_names = ['gradReject', 'magReject', 'eegReject', 'eogReject']
-    reject_pynames = ['grad', 'mag', 'eeg', 'eog']
-    reject_params = dict()
+    reject_names = ['gradReject', 'magReject', 'eegReject', 'eogReject', 'ecgReject']
+    reject_pynames = ['grad', 'mag', 'eeg', 'eog', 'ecg']
+    reject = dict()
     for line in lines:
         words = line.split()
         print words
         if words[0] in reject_names:
-            reject_params[reject_pynames[reject_names.index(words[0])]] = \
+            reject[reject_pynames[reject_names.index(words[0])]] = \
                                                                 float(words[1])
 
-    return reject_params
\ No newline at end of file
+    return reject
+
+def read_flat_parameters(fname):
+    """Read flat channel rejection parameters from .cov or .ave config file"""
+
+    try:
+        with open(fname, 'r') as f:
+            lines = f.readlines()
+    except:
+        print("Error while reading %s" % fname)
+
+    reject_names = ['gradFlat', 'magFlat', 'eegFlat', 'eogFlat', 'ecgFlat']
+    reject_pynames = ['grad', 'mag', 'eeg', 'eog', 'ecg']
+    flat = dict()
+    for line in lines:
+        words = line.split()
+        print words
+        if words[0] in reject_names:
+            flat[reject_pynames[reject_names.index(words[0])]] = \
+                                                                float(words[1])
+
+    return flat
diff --git a/mne/tests/test_epochs.py b/mne/tests/test_epochs.py
index df1170f..6c78989 100644
--- a/mne/tests/test_epochs.py
+++ b/mne/tests/test_epochs.py
@@ -44,12 +44,14 @@ def test_reject_epochs():
     picks = fiff.pick_types(raw.info, meg=True, eeg=True, stim=True,
                             eog=True, include=['STI 014'])
     epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
-                        baseline=(None, 0))
-    n_epochs = len(epochs)
-    epochs.reject(grad=1000e-12, mag=4e-12, eeg=80e-6, eog=150e-6)
-    n_clean_epochs = len(epochs)
+                        baseline=(None, 0),
+                        reject=dict(grad=1000e-12, mag=4e-12, eeg=80e-6,
+                                    eog=150e-6))
+    data = epochs.get_data()
+    n_events = len(epochs.events)
+    n_clean_epochs = len(data)
     # Should match
     # mne_process_raw --raw test_raw.fif --projoff \
     #   --saveavetag -ave --ave test.ave --filteroff
-    assert n_epochs > n_clean_epochs
+    assert n_events > n_clean_epochs
     assert n_clean_epochs == 3

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git



More information about the debian-med-commit mailing list