[med-svn] [python-mne] 118/353: FIX: cov comp. for multiple event types if keep_sample_mean=False
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:24:40 UTC 2015
This is an automated email from the git hooks/post-receive script.
yoh pushed a commit to tag 0.4
in repository python-mne.
commit 2e287f47470798c1221f81ca9758989b124d3c1e
Author: Martin Luessi <mluessi at nmr.mgh.harvard.edu>
Date: Tue Mar 13 14:44:22 2012 -0400
FIX: cov comp. for multiple event types if keep_sample_mean=False
---
mne/cov.py | 102 ++++++++++++++++++++++++++++++++++++--------------
mne/fiff/__init__.py | 2 +-
mne/fiff/proj.py | 17 ++++++++-
mne/tests/test_cov.py | 17 +++++++--
4 files changed, 105 insertions(+), 33 deletions(-)
diff --git a/mne/cov.py b/mne/cov.py
index 053a261..c5689a4 100644
--- a/mne/cov.py
+++ b/mne/cov.py
@@ -6,12 +6,14 @@
import copy
import os
from math import floor, ceil
+import warnings
+
import numpy as np
from scipy import linalg
from . import fiff
from .fiff.write import start_file, end_file
-from .fiff.proj import make_projector
+from .fiff.proj import make_projector, proj_equal
from .fiff import fiff_open
from .fiff.pick import pick_types, channel_indices_by_type
from .fiff.constants import FIFF
@@ -231,48 +233,92 @@ def compute_covariance(epochs, keep_sample_mean=True):
The noise covariance is typically estimated on pre-stim periods
when the stim onset if defined from events.
+ If the covariance is computed for multiple event types (events
+ with different IDs), an Epochs object for each event type has to
+ be created and a list of Epochs has to be passed to this function.
+
+ Note: Baseline correction should be used when creating the Epochs.
+ Otherwise the computed covariance matrix will be inaccurate.
+
+ Note: For multiple event types, it is also possible to create a
+ single Epochs object with events obtained using
+ merge_events(). However, the resulting covariance matrix
+ will only be correct if keep_sample_mean is True.
+
Parameters
----------
- epochs : instance of Epochs
+ epochs : instance of Epochs, or a list of Epochs objects
The epochs
keep_sample_mean : bool
- If False data are centered at each instant before computing
- the covariance.
+ If False, the average response over epochs is computed for
+ each event type and subtracted during the covariance
+ computation. This is useful if the evoked response from a
+ previous stimulus extends into the baseline period of the next.
+
Returns
-------
cov : instance of Covariance
The computed covariance.
"""
+ if not isinstance(epochs, list):
+ epochs = [epochs]
+
+ # check for baseline correction
+ for epochs_t in epochs:
+ if epochs_t.baseline is None:
+ warnings.warn('Epochs are not baseline corrected, covariance '
+ 'matrix may be inaccurate')
+
+ bads = epochs[0].info['bads']
+ projs = epochs[0].info['projs']
+ ch_names = epochs[0].ch_names
+
+ # make sure Epochs are compatible
+ for epochs_t in epochs[1:]:
+ if epochs_t.info['bads'] != bads:
+ raise ValueError('Epochs must have same bad channels')
+ if epochs_t.ch_names != ch_names:
+ raise ValueError('Epochs must have same channel names')
+ for proj_a, proj_b in zip(epochs_t.info['projs'], projs):
+ if not proj_equal(proj_a, proj_b):
+ raise ValueError('Epochs must have same projectors')
+
+ n_epoch_types = len(epochs)
data = 0.0
- data_mean = 0.0
- n_samples = 0
- n_epochs = 0
- picks_meeg = pick_types(epochs.info, meg=True, eeg=True, eog=False)
- ch_names = [epochs.ch_names[k] for k in picks_meeg]
- for e in epochs:
- e = e[picks_meeg]
- if not keep_sample_mean:
- data_mean += e
- data += np.dot(e, e.T)
- n_samples += e.shape[1]
- n_epochs += 1
-
- if n_samples == 0:
+ data_mean = list(np.zeros(n_epoch_types))
+ n_samples = np.zeros(n_epoch_types, dtype=np.int)
+ n_epochs = np.zeros(n_epoch_types, dtype=np.int)
+
+ picks_meeg = pick_types(epochs[0].info, meg=True, eeg=True, eog=False)
+ ch_names = [epochs[0].ch_names[k] for k in picks_meeg]
+ for i, epochs_t in enumerate(epochs):
+ for e in epochs_t:
+ e = e[picks_meeg]
+ if not keep_sample_mean:
+ data_mean[i] += e
+ data += np.dot(e, e.T)
+ n_samples[i] += e.shape[1]
+ n_epochs[i] += 1
+
+ n_samples_tot = int(np.sum(n_samples))
+
+ if n_samples_tot == 0:
raise ValueError('Not enough samples to compute the noise covariance'
- ' matrix : %d samples' % n_samples)
+ ' matrix : %d samples' % n_samples_tot)
if keep_sample_mean:
- nfree = n_samples
- data /= nfree
+ data /= n_samples_tot
else:
n_samples_epoch = n_samples / n_epochs
- nfree = n_samples_epoch * (n_epochs - 1)
- data /= nfree
- data -= 1.0 / nfree * np.dot(data_mean, data_mean.T)
+ norm_const = np.sum(n_samples_epoch * (n_epochs - 1))
+ for i, mean in enumerate(data_mean):
+ data -= 1.0 / n_epochs[i] * np.dot(mean, mean.T)
+ data /= norm_const
+
cov = Covariance(None)
cov.data = data
cov.ch_names = ch_names
- cov.nfree = nfree
+ cov.nfree = n_samples_tot
# XXX : do not compute eig and eigvec now (think it's better...)
eig = None
@@ -280,11 +326,11 @@ def compute_covariance(epochs, keep_sample_mean=True):
# Store structure for fif
cov._cov = dict(kind=1, diag=False, dim=len(data), names=ch_names,
- data=data, projs=copy.deepcopy(epochs.info['projs']),
- bads=epochs.info['bads'], nfree=n_samples, eig=eig,
+ data=data, projs=copy.deepcopy(epochs[0].info['projs']),
+ bads=epochs[0].info['bads'], nfree=n_samples_tot, eig=eig,
eigvec=eigvec)
- print "Number of samples used : %d" % n_samples
+ print "Number of samples used : %d" % n_samples_tot
print '[done]'
return cov
diff --git a/mne/fiff/__init__.py b/mne/fiff/__init__.py
index 1e0fc1c..3738047 100644
--- a/mne/fiff/__init__.py
+++ b/mne/fiff/__init__.py
@@ -15,5 +15,5 @@ from .pick import pick_types, pick_channels, pick_types_evoked, \
pick_types_forward
from .compensator import get_current_comp
-from .proj import compute_spatial_vectors
+from .proj import compute_spatial_vectors, proj_equal
from .cov import read_cov, write_cov
diff --git a/mne/fiff/proj.py b/mne/fiff/proj.py
index e667d68..81032e7 100644
--- a/mne/fiff/proj.py
+++ b/mne/fiff/proj.py
@@ -25,6 +25,21 @@ class Projection(dict):
return "Projection (%s)" % s
+def proj_equal(a, b):
+ """ Test if two projectors are equal """
+
+ equal = a['active'] == b['active']\
+ and a['kind'] == b['kind']\
+ and a['desc'] == b['desc']\
+ and a['data']['col_names'] == b['data']['col_names']\
+ and a['data']['row_names'] == b['data']['row_names']\
+ and a['data']['ncol'] == b['data']['ncol']\
+ and a['data']['nrow'] == b['data']['nrow']\
+ and np.all(a['data']['data'] == b['data']['data'])
+
+ return equal
+
+
def read_proj(fid, node):
"""Read spatial projections from a FIF file.
@@ -340,7 +355,7 @@ def compute_spatial_vectors(epochs, n_grad=2, n_mag=2, n_eeg=2):
['planar', 'axial', 'eeg']):
if n == 0:
continue
- data_ind = data[ind][:,ind]
+ data_ind = data[ind][:, ind]
U = linalg.svd(data_ind, full_matrices=False,
overwrite_a=True)[0][:, :n]
for k, u in enumerate(U.T):
diff --git a/mne/tests/test_cov.py b/mne/tests/test_cov.py
index af405a5..901f7ed 100644
--- a/mne/tests/test_cov.py
+++ b/mne/tests/test_cov.py
@@ -57,8 +57,9 @@ def test_cov_estimation_with_triggers():
event_ids = [1, 2, 3, 4]
reject = dict(grad=10000e-13, mag=4e-12, eeg=80e-6, eog=150e-6)
- events = merge_events(events, event_ids, 1234)
- epochs = Epochs(raw, events, 1234, tmin=-0.2, tmax=0,
+ # cov with merged events and keep_sample_mean=True
+ events_merged = merge_events(events, event_ids, 1234)
+ epochs = Epochs(raw, events_merged, 1234, tmin=-0.2, tmax=0,
baseline=(-0.2, -0.1), proj=True,
reject=reject)
@@ -68,11 +69,21 @@ def test_cov_estimation_with_triggers():
assert_true((linalg.norm(cov.data - cov_mne.data, ord='fro')
/ linalg.norm(cov.data, ord='fro')) < 0.005)
+ # cov using a list of epochs and keep_sample_mean=True
+ epochs = [Epochs(raw, events, ev_id, tmin=-0.2, tmax=0,
+ baseline=(-0.2, -0.1), proj=True, reject=reject)
+ for ev_id in event_ids]
+
+ cov2 = compute_covariance(epochs, keep_sample_mean=True)
+ assert_array_almost_equal(cov.data, cov2.data)
+ assert_true(cov.ch_names == cov2.ch_names)
+
+ # cov with keep_sample_mean=False using a list of epochs
cov = compute_covariance(epochs, keep_sample_mean=False)
cov_mne = Covariance(cov_fname)
assert_true(cov_mne.ch_names == cov.ch_names)
assert_true((linalg.norm(cov.data - cov_mne.data, ord='fro')
- / linalg.norm(cov.data, ord='fro')) < 0.06)
+ / linalg.norm(cov.data, ord='fro')) < 0.005)
# test IO when computation done in Python
cov.save('test-cov.fif') # test saving
--
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