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 39dd68eed8551797f6f27babe22f705068a0b6b8 Author: Alexandre Gramfort <[email protected]> Date: Tue Apr 17 14:42:01 2012 +0200 ENH : add new function to regularize noise covariance --- mne/cov.py | 87 +++++++++++++++++++++++++++++++++- mne/minimum_norm/tests/test_inverse.py | 1 + mne/tests/test_cov.py | 15 ++++++ mne/viz.py | 15 ++++-- 4 files changed, 112 insertions(+), 6 deletions(-) diff --git a/mne/cov.py b/mne/cov.py index aad0d8a..389599e 100644 --- a/mne/cov.py +++ b/mne/cov.py @@ -13,9 +13,9 @@ from scipy import linalg from . import fiff from .fiff.write import start_file, end_file -from .fiff.proj import make_projector, proj_equal +from .fiff.proj import make_projector, proj_equal, activate_proj from .fiff import fiff_open -from .fiff.pick import pick_types, channel_indices_by_type +from .fiff.pick import pick_types, channel_indices_by_type, pick_channels_cov from .fiff.constants import FIFF from .epochs import _is_good @@ -423,3 +423,86 @@ def prepare_noise_cov(noise_cov, info, ch_names): diag=False, names=ch_names) return noise_cov + + +def regularize(cov, info, mag=0.1, grad=0.1, eeg=0.1, exclude=None, + proj=True): + """Regularize noise covariance matrix + + This method works by adding a constant to the diagonal for each + channel type separatly. Special care is taken to keep the + rank of the data constant. + + Parameters + ---------- + cov: Covariance + The noise covariance matrix. + info: dict + The measurement info (used to get channel types and bad channels) + mag: float + Regularization factor for MEG magnetometers + grad: float + Regularization factor for MEG gradiometers + eeg: float + Regularization factor for EEG + exclude: list + List of channels to mark as bad. If None, bads channels + are extracted from info and cov['bads']. + proj: bool + Apply or not projections to keep rank of data. + + Return + ------ + reg_cov : Covariance + The regularized covariance matrix. + """ + if exclude is None: + exclude = info['bads'] + cov['bads'] + + sel_eeg = pick_types(info, meg=False, eeg=True, exclude=exclude) + sel_mag = pick_types(info, meg='mag', eeg=False, exclude=exclude) + sel_grad = pick_types(info, meg='grad', eeg=False, exclude=exclude) + + info_ch_names = info['ch_names'] + cov = pick_channels_cov(cov, include=info_ch_names, exclude=exclude) + ch_names = cov.ch_names + idx_eeg = [ch_names.index(info_ch_names[c]) for c in sel_eeg] + idx_mag = [ch_names.index(info_ch_names[c]) for c in sel_mag] + idx_grad = [ch_names.index(info_ch_names[c]) for c in sel_grad] + + C = cov['data'] + + assert len(C) == (len(idx_eeg) + len(idx_mag) + len(idx_grad)) + + if proj: + projs = info['projs'] + cov['projs'] + projs = activate_proj(projs) + + for desc, idx, reg in [('EEG', idx_eeg, eeg), ('MAG', idx_mag, mag), + ('GRAD', idx_grad, grad)]: + if len(idx) == 0 or reg == 0.0: + print " %s regularization : None" % desc + continue + + print " %s regularization : %s" % (desc, reg) + + this_C = C[idx][:, idx] + if proj: + this_ch_names = [ch_names[k] for k in idx] + P, ncomp, _ = make_projector(projs, this_ch_names) + U = linalg.svd(P)[0][:, :-ncomp] + if ncomp > 0: + print ' Created an SSP operator for %s (dimension = %d)' % \ + (desc, ncomp) + this_C = np.dot(U.T, np.dot(this_C, U)) + + sigma = np.mean(np.diag(this_C)) + this_C.flat[::len(this_C) + 1] += reg * sigma # modify diag inplace + if proj and ncomp > 0: + this_C = np.dot(U, np.dot(this_C, U.T)) + + C[np.ix_(idx, idx)] = this_C + + cov['data'] = C + + return cov diff --git a/mne/minimum_norm/tests/test_inverse.py b/mne/minimum_norm/tests/test_inverse.py index 8b896df..793f59e 100644 --- a/mne/minimum_norm/tests/test_inverse.py +++ b/mne/minimum_norm/tests/test_inverse.py @@ -108,6 +108,7 @@ def test_apply_inverse_operator(): # Test MNE inverse computation starting from forward operator evoked = fiff.Evoked(fname_data, setno=0, baseline=(None, 0)) fwd_op = read_forward_solution(fname_fwd, surf_ori=True) + my_inv_op = make_inverse_operator(evoked.info, fwd_op, noise_cov, loose=0.2, depth=0.8) diff --git a/mne/tests/test_cov.py b/mne/tests/test_cov.py index f1fe6aa..95c1ac7 100644 --- a/mne/tests/test_cov.py +++ b/mne/tests/test_cov.py @@ -2,8 +2,10 @@ import os.path as op from nose.tools import assert_true from numpy.testing import assert_array_almost_equal +import numpy as np from scipy import linalg +from ..cov import regularize from .. import Covariance, Epochs, merge_events, \ find_events, compute_raw_data_covariance, \ compute_covariance @@ -113,3 +115,16 @@ def test_arithmetic_cov(): assert_array_almost_equal(cov_sum.nfree, cov.nfree) assert_array_almost_equal(cov_sum.data, cov.data) assert_true(cov_sum.ch_names == cov.ch_names) + + +def test_regularize_cov(): + """Test cov regularization + """ + noise_cov = Covariance(cov_fname) + raw = Raw(raw_fname) + # Regularize noise cov + reg_noise_cov = regularize(noise_cov, raw.info, + mag=0.1, grad=0.1, eeg=0.1, proj=True) + assert_true(noise_cov['dim'] == reg_noise_cov['dim']) + assert_true(noise_cov['data'].shape == reg_noise_cov['data'].shape) + assert_true(np.mean(noise_cov['data'] < reg_noise_cov['data']) < 0.08) diff --git a/mne/viz.py b/mne/viz.py index aa651b3..2933a87 100644 --- a/mne/viz.py +++ b/mne/viz.py @@ -85,7 +85,7 @@ def plot_evoked(evoked, picks=None, unit=True, show=True): pl.show() -def plot_cov(cov, info, exclude=None, colorbar=True, proj=False, show_svd=True, +def plot_cov(cov, info, exclude=[], colorbar=True, proj=False, show_svd=True, show=True): """Plot Covariance data @@ -113,9 +113,12 @@ def plot_cov(cov, info, exclude=None, colorbar=True, proj=False, show_svd=True, sel_eeg = pick_types(info, meg=False, eeg=True, exclude=exclude) sel_mag = pick_types(info, meg='mag', eeg=False, exclude=exclude) sel_grad = pick_types(info, meg='grad', eeg=False, exclude=exclude) - idx_eeg = [ch_names.index(info_ch_names[c]) for c in sel_eeg] - idx_mag = [ch_names.index(info_ch_names[c]) for c in sel_mag] - idx_grad = [ch_names.index(info_ch_names[c]) for c in sel_grad] + idx_eeg = [ch_names.index(info_ch_names[c]) + for c in sel_eeg if info_ch_names[c] in ch_names] + idx_mag = [ch_names.index(info_ch_names[c]) + for c in sel_mag if info_ch_names[c] in ch_names] + idx_grad = [ch_names.index(info_ch_names[c]) + for c in sel_grad if info_ch_names[c] in ch_names] idx_names = [(idx_eeg, 'EEG covariance', 'uV', 1e6), (idx_grad, 'Gradiometers', 'fT/cm', 1e13), @@ -162,6 +165,10 @@ def plot_cov(cov, info, exclude=None, colorbar=True, proj=False, show_svd=True, pl.imshow(C[idx][:, idx], interpolation="nearest") pl.title(name) pl.subplots_adjust(0.04, 0.0, 0.98, 0.94, 0.2, 0.26) + try: + pl.tight_layout() # XXX : recent pylab feature + except: + pass if show: pl.show() -- 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
