[med-svn] [python-mne] 153/353: ENH : add new function to regularize noise covariance
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:24:48 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 39dd68eed8551797f6f27babe22f705068a0b6b8
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
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
More information about the debian-med-commit
mailing list