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 <[email protected]> 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 _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
