pure python creation of inverse operator

Makefile

@@ -33,11 +33,10 @@ test-doc:
 	$(NOSETESTS) --with-doctest --doctest-tests --doctest-extension=rst doc/ doc/source/
 
 test-coverage:
-	$(NOSETESTS) --with-coverage --cover-package=mne
-
+	$(NOSETESTS) --with-coverage --cover-package=mne --cover-html --cover-html-dir=coverage
 
 trailing-spaces:
-	find -name "*.py" |xargs sed -i 's/[ \t]*$$//'
+	find . -name "*.py" | xargs perl -pi -e 's/[ \t]*$$//'
 
 ctags:
 	# make tags for symbol based navigation in emacs and vim

examples/inverse/plot_make_inverse_operator.py

@@ -0,0 +1,53 @@
+"""
+===============================================================
+Assemble inverse operator and compute MNE-dSPM inverse solution
+===============================================================
+
+Assemble an EEG inverse operator and compute dSPM inverse solution
+on MNE evoked dataset and stores the solution in stc files for
+visualisation.
+
+"""
+
+# Author: Alexandre Gramfort <gramfort@nmr.mgh.harvard.edu>
+#
+# License: BSD (3-clause)
+
+print __doc__
+
+import pylab as pl
+import mne
+from mne.datasets import sample
+from mne.fiff import Evoked
+from mne.minimum_norm import make_inverse_operator, apply_inverse
+
+data_path = sample.data_path('..')
+fname_fwd = data_path + '/MEG/sample/sample_audvis-eeg-oct-6-fwd.fif'
+fname_cov = data_path + '/MEG/sample/sample_audvis-cov.fif'
+fname_evoked = data_path + '/MEG/sample/sample_audvis-ave.fif'
+
+setno = 0
+snr = 3.0
+lambda2 = 1.0 / snr ** 2
+dSPM = True
+
+# Load data
+evoked = Evoked(fname_evoked, setno=setno, baseline=(None, 0))
+forward = mne.read_forward_solution(fname_fwd, surf_ori=True)
+noise_cov = mne.Covariance(fname_cov)
+inverse_operator = make_inverse_operator(evoked.info, forward, noise_cov,
+                                         loose=0.2, depth=0.8)
+
+# Compute inverse solution
+stc = apply_inverse(evoked, inverse_operator, lambda2, dSPM, pick_normal=False)
+
+# Save result in stc files
+stc.save('mne_dSPM_inverse')
+
+###############################################################################
+# View activation time-series
+pl.close('all')
+pl.plot(1e3 * stc.times, stc.data[::150, :].T)
+pl.xlabel('time (ms)')
+pl.ylabel('dSPM value')
+pl.show()

mne/cov.py

@@ -23,31 +23,6 @@
 from .epochs import _is_good
 
 
-def rank(A, tol=1e-8):
-    s = linalg.svd(A, compute_uv=0)
-    return np.sum(np.where(s > s[0] * tol, 1, 0))
-
-
-def _get_whitener(A, rnk, pca, ch_type):
-    # whitening operator
-    D, V = linalg.eigh(A, overwrite_a=True)
-    I = np.argsort(D)[::-1]
-    D = D[I]
-    V = V[:, I]
-    D = 1.0 / D
-    if not pca:  # No PCA case.
-        print 'Not doing PCA for %s.' % ch_type
-        W = np.sqrt(D)[:, None] * V.T
-    else:  # Rey's approach. MNE has been changed to implement this.
-        print 'Setting small %s eigenvalues to zero.' % ch_type
-        D[rnk:] = 0.0
-        W = np.sqrt(D)[:, None] * V.T
-        # This line will reduce the actual number of variables in data
-        # and leadfield to the true rank.
-        W = W[:rnk]
-    return W
-
-
 class Covariance(object):
     """Noise covariance matrix"""
 
@@ -79,155 +54,16 @@ def save(self, fname):
         """save covariance matrix in a FIF file"""
         write_cov_file(fname, self._cov)
 
-    def get_whitener(self, info, mag_reg=0.1, grad_reg=0.1, eeg_reg=0.1,
-                     pca=True):
-        """Compute whitener based on a list of channels
-
-        Parameters
-        ----------
-        info : dict
-            Measurement info of data to apply the whitener.
-            Defines data channels and which are the bad channels
-            to be ignored.
-        mag_reg : float
-            Regularization of the magnetometers.
-            Recommended between 0.05 and 0.2
-        grad_reg : float
-            Regularization of the gradiometers.
-            Recommended between 0.05 and 0.2
-        eeg_reg : float
-            Regularization of the EGG channels.
-            Recommended between 0.05 and 0.2
-        pca : bool
-            If True, whitening is restricted to the space of
-            the data. It makes sense when data have a low rank
-            due to SSP or maxfilter.
-
-        Returns
-        -------
-        W : instance of Whitener
-        """
-
-        if not 0 <= grad_reg <= 1:
-            raise ValueError('grad_reg should be a scalar between 0 and 1')
-        if not 0 <= mag_reg <= 1:
-            raise ValueError('mag_reg should be a scalar between 0 and 1')
-        if not 0 <= eeg_reg <= 1:
-            raise ValueError('eeg_reg should be a scalar between 0 and 1')
-
-        if pca and self.kind == 'diagonal':
-            print "Setting pca to False with a diagonal covariance matrix."
-            pca = False
-
-        bads = info['bads']
-        C_idx = [k for k, name in enumerate(self.ch_names)
-                 if name in info['ch_names'] and name not in bads]
-        ch_names = [self.ch_names[k] for k in C_idx]
-        C_noise = self.data[np.ix_(C_idx, C_idx)]  # take covariance submatrix
-
-        # Create the projection operator
-        proj, ncomp, _ = make_projector(info['projs'], ch_names)
-        if ncomp > 0:
-            print '\tCreated an SSP operator (subspace dimension = %d)' % ncomp
-            C_noise = np.dot(proj, np.dot(C_noise, proj.T))
-
-        # Regularize Noise Covariance Matrix.
-        variances = np.diag(C_noise)
-        ind_meg = pick_types(info, meg=True, eeg=False, exclude=bads)
-        names_meg = [info['ch_names'][k] for k in ind_meg]
-        C_ind_meg = [ch_names.index(name) for name in names_meg]
-
-        ind_grad = pick_types(info, meg='grad', eeg=False, exclude=bads)
-        names_grad = [info['ch_names'][k] for k in ind_grad]
-        C_ind_grad = [ch_names.index(name) for name in names_grad]
-
-        ind_mag = pick_types(info, meg='mag', eeg=False, exclude=bads)
-        names_mag = [info['ch_names'][k] for k in ind_mag]
-        C_ind_mag = [ch_names.index(name) for name in names_mag]
-
-        ind_eeg = pick_types(info, meg=False, eeg=True, exclude=bads)
-        names_eeg = [info['ch_names'][k] for k in ind_eeg]
-        C_ind_eeg = [ch_names.index(name) for name in names_eeg]
-
-        has_meg = len(ind_meg) > 0
-        has_eeg = len(ind_eeg) > 0
-
-        if self.kind == 'diagonal':
-            C_noise = np.diag(variances)
-            rnkC_noise = len(variances)
-            print 'Rank of noise covariance is %d' % rnkC_noise
-        else:
-            # estimate noise covariance matrix rank
-            # Loop on all the required data types (MEG MAG, MEG GRAD, EEG)
-
-            if has_meg:  # Separate rank of MEG
-                rank_meg = rank(C_noise[C_ind_meg][:, C_ind_meg])
-                print 'Rank of MEG part of noise covariance is %d' % rank_meg
-            if has_eeg:  # Separate rank of EEG
-                rank_eeg = rank(C_noise[C_ind_eeg][:, C_ind_eeg])
-                print 'Rank of EEG part of noise covariance is %d' % rank_eeg
-
-            for ind, reg in zip([C_ind_grad, C_ind_mag, C_ind_eeg],
-                                [grad_reg, mag_reg, eeg_reg]):
-                if len(ind) > 0:
-                    # add constant on diagonal
-                    C_noise[ind, ind] += reg * np.mean(variances[ind])
-
-            if has_meg and has_eeg:  # Sets cross terms to zero
-                C_noise[np.ix_(C_ind_meg, C_ind_eeg)] = 0.0
-                C_noise[np.ix_(C_ind_eeg, C_ind_meg)] = 0.0
-
-        # whitening operator
-        if has_meg:
-            W_meg = _get_whitener(C_noise[C_ind_meg][:, C_ind_meg], rank_meg,
-                                  pca, 'MEG')
-
-        if has_eeg:
-            W_eeg = _get_whitener(C_noise[C_ind_eeg][:, C_ind_eeg], rank_eeg,
-                                  pca, 'EEG')
-
-        if has_meg and not has_eeg:  # Only MEG case.
-            W = W_meg
-        elif has_eeg and not has_meg:  # Only EEG case.
-            W = W_eeg
-        elif has_eeg and has_meg:  # Bimodal MEG and EEG case.
-            # Whitening of MEG and EEG separately, which assumes zero
-            # covariance between MEG and EEG (i.e., a block diagonal noise
-            # covariance). This was recommended by Matti as EEG does not
-            # measure all the signals from the same environmental noise sources
-            # as MEG.
-            W = np.r_[np.c_[W_meg, np.zeros((W_meg.shape[0], W_eeg.shape[1]))],
-                      np.c_[np.zeros((W_eeg.shape[0], W_meg.shape[1])), W_eeg]]
-
-        whitener = Whitener(W, names_meg + names_eeg)
-        return whitener
-
     def __repr__(self):
         s = "kind : %s" % self.kind
         s += ", size : %s x %s" % self.data.shape
         s += ", data : %s" % self.data
         return "Covariance (%s)" % s
 
 
-class Whitener(object):
-    """Whitener
-
-    Attributes
-    ----------
-    W : array
-        Whiten matrix
-    ch_names : list of strings
-        Channel names (columns of W)
-    """
-
-    def __init__(self, W, ch_names):
-        self.W = W
-        self.ch_names = ch_names
-
 ###############################################################################
 # IO
 
-
 def read_cov(fid, node, cov_kind):
     """Read a noise covariance matrix
 
@@ -336,6 +172,7 @@ def read_cov(fid, node, cov_kind):
 
     return None
 
+
 ###############################################################################
 # Estimate from data
 
@@ -554,3 +391,75 @@ def write_cov_file(fname, cov):
         raise '%s', inst
 
     end_file(fid)
+
+
+###############################################################################
+# Prepare for inverse modeling
+
+def rank(A, tol=1e-8):
+    s = linalg.svd(A, compute_uv=0)
+    return np.sum(np.where(s > s[0] * tol, 1, 0))
+
+
+def _get_whitener(A, pca, ch_type):
+    # whitening operator
+    rnk = rank(A)
+    eig, eigvec = linalg.eigh(A, overwrite_a=True)
+    eigvec = eigvec.T
+    eig[:-rnk] = 0.0
+    print 'Setting small %s eigenvalues to zero.' % ch_type
+    if not pca:  # No PCA case.
+        print 'Not doing PCA for %s.' % ch_type
+    else:
+        print 'Doing PCA for %s.' % ch_type
+        # This line will reduce the actual number of variables in data
+        # and leadfield to the true rank.
+        eigvec = eigvec[:-rnk].copy()
+    return eig, eigvec
+
+
+def prepare_noise_cov(noise_cov, info, ch_names):
+    C_ch_idx = [noise_cov.ch_names.index(c) for c in ch_names]
+    C = noise_cov.data[C_ch_idx][:, C_ch_idx]
+
+    # Create the projection operator
+    proj, ncomp, _ = make_projector(info['projs'], ch_names)
+    if ncomp > 0:
+        print '\tCreated an SSP operator (subspace dimension = %d)' % ncomp
+        C = np.dot(proj, np.dot(C, proj.T))
+
+    pick_meg = pick_types(info, meg=True, eeg=False, exclude=info['bads'])
+    pick_eeg = pick_types(info, meg=False, eeg=True, exclude=info['bads'])
+    meg_names = [info['chs'][k]['ch_name'] for k in pick_meg]
+    C_meg_idx = [k for k in range(len(C)) if ch_names[k] in meg_names]
+    eeg_names = [info['chs'][k]['ch_name'] for k in pick_eeg]
+    C_eeg_idx = [k for k in range(len(C)) if ch_names[k] in eeg_names]
+
+    has_meg = len(C_meg_idx) > 0
+    has_eeg = len(C_eeg_idx) > 0
+
+    if has_meg:
+        C_meg = C[C_meg_idx][:, C_meg_idx]
+        C_meg_eig, C_meg_eigvec = _get_whitener(C_meg, False, 'MEG')
+
+    if has_eeg:
+        C_eeg = C[C_eeg_idx][:, C_eeg_idx]
+        C_eeg_eig, C_eeg_eigvec = _get_whitener(C_eeg, False, 'EEG')
+
+    n_chan = len(ch_names)
+    eigvec = np.zeros((n_chan, n_chan), dtype=np.float)
+    eig = np.zeros(n_chan, dtype=np.float)
+
+    if has_meg:
+        eigvec[np.ix_(C_meg_idx, C_meg_idx)] = C_meg_eigvec
+        eig[C_meg_idx] = C_meg_eig
+    if has_eeg:
+        eigvec[np.ix_(C_eeg_idx, C_eeg_idx)] = C_eeg_eigvec
+        eig[C_eeg_idx] = C_eeg_eig
+
+    assert(len(C_meg_idx) + len(C_eeg_idx) == n_chan)
+
+    noise_cov = dict(data=C, eig=eig, eigvec=eigvec, dim=len(ch_names),
+                     diag=False, names=ch_names)
+
+    return noise_cov