diff --git a/Makefile b/Makefile index 09cfd67ec9e..5f2492dcc77 100755 --- a/Makefile +++ b/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 diff --git a/examples/inverse/plot_make_inverse_operator.py b/examples/inverse/plot_make_inverse_operator.py new file mode 100644 index 00000000000..c5bab4d9cc4 --- /dev/null +++ b/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 +# +# 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() diff --git a/mne/cov.py b/mne/cov.py index 9fd121580e5..e5f7a48066a 100644 --- a/mne/cov.py +++ b/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,129 +54,6 @@ 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 @@ -209,25 +61,9 @@ def __repr__(self): 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 diff --git a/mne/forward.py b/mne/forward.py index 6b51ec523ba..61d017db125 100644 --- a/mne/forward.py +++ b/mne/forward.py @@ -341,7 +341,7 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, fwd['src'] = src # Handle the source locations and orientations - if fwd['source_ori'] == FIFF.FIFFV_MNE_FIXED_ORI or force_fixed == True: + if (fwd['source_ori'] == FIFF.FIFFV_MNE_FIXED_ORI) or force_fixed: nuse = 0 fwd['source_rr'] = np.zeros((fwd['nsource'], 3)) fwd['source_nn'] = np.zeros((fwd['nsource'], 3)) @@ -366,50 +366,49 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, fwd['sol_grad']['ncol'] = 3 * fwd['nsource'] print '[done]' + elif surf_ori: + # Rotate the local source coordinate systems + print '\tConverting to surface-based source orientations...' + nuse = 0 + pp = 0 + nuse_total = sum([s['nuse'] for s in src]) + fwd['source_rr'] = np.zeros((fwd['nsource'], 3)) + fwd['source_nn'] = np.empty((3 * nuse_total, 3), dtype=np.float) + for s in src: + rr = s['rr'][s['vertno'], :] + fwd['source_rr'][nuse:nuse + s['nuse'], :] = rr + for p in range(s['nuse']): + # Project out the surface normal and compute SVD + nn = s['nn'][s['vertno'][p], :][:, None] + U, S, _ = linalg.svd(np.eye(3, 3) - nn * nn.T) + # Make sure that ez is in the direction of nn + if np.sum(nn.ravel() * U[:, 2].ravel()) < 0: + U *= -1.0 + fwd['source_nn'][pp:pp + 3, :] = U.T + pp += 3 + + nuse += s['nuse'] + + surf_rot = _block_diag(fwd['source_nn'].T, 3) + fwd['sol']['data'] = fwd['sol']['data'] * surf_rot + if fwd['sol_grad'] is not None: + fwd['sol_grad']['data'] = np.dot(fwd['sol_grad']['data'] * \ + np.kron(surf_rot, np.eye(3))) + + print '[done]' else: - if surf_ori: - # Rotate the local source coordinate systems - print '\tConverting to surface-based source orientations...' - nuse = 0 - pp = 0 - nuse_total = sum([s['nuse'] for s in src]) - fwd['source_rr'] = np.zeros((fwd['nsource'], 3)) - fwd['source_nn'] = np.empty((3 * nuse_total, 3), dtype=np.float) - for s in src: - fwd['source_rr'][nuse:nuse + s['nuse'], :] = \ - s['rr'][s['vertno'], :] - for p in range(s['nuse']): - # Project out the surface normal and compute SVD - nn = s['nn'][s['vertno'][p], :].T - nn = nn[:, None] - U, S, _ = linalg.svd(np.eye(3, 3) - nn * nn.T) - # Make sure that ez is in the direction of nn - if np.sum(nn * U[:, 2]) < 0: - U *= -1 - - fwd['source_nn'][pp:pp + 3, :] = U.T - pp += 3 - - nuse += s['nuse'] - - surf_rot = _block_diag(fwd['source_nn'].T, 3) - fwd['sol']['data'] = fwd['sol']['data'] * surf_rot - if fwd['sol_grad'] is not None: - fwd['sol_grad']['data'] = np.dot(fwd['sol_grad']['data'] * \ - np.kron(surf_rot, np.eye(3))) + print '\tCartesian source orientations...' + nuse = 0 + fwd['source_rr'] = np.zeros((fwd['nsource'], 3)) + for s in src: + rr = s['rr'][s['vertno'], :] + fwd['source_rr'][nuse:nuse + s['nuse'], :] = rr + nuse += s['nuse'] - print '[done]' - else: - print '\tCartesian source orientations...' - nuse = 0 - fwd['source_rr'] = np.zeros((fwd['nsource'], 3)) - for s in src: - fwd['source_rr'][nuse:nuse + s['nuse'], :] = \ - s['rr'][s['vertno'], :] - nuse += s['nuse'] - - fwd['source_nn'] = np.kron(np.ones((fwd['nsource'], 1)), np.eye(3)) - print '[done]' + fwd['source_nn'] = np.kron(np.ones((fwd['nsource'], 1)), np.eye(3)) + print '[done]' + + fwd['surf_ori'] = surf_ori # Add channel information fwd['chs'] = chs @@ -417,3 +416,19 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, fwd = pick_channels_forward(fwd, include=include, exclude=exclude) return fwd + + +def compute_depth_prior(G, exp=0.8, limit=10.0): + """Compute weighting for depth prior + """ + n_pos = G.shape[1] // 3 + d = np.zeros(n_pos) + for k in xrange(n_pos): + Gk = G[:, 3 * k:3 * (k + 1)] + d[k] = linalg.svdvals(np.dot(Gk.T, Gk))[0] + w = 1.0 / d + wmax = np.min(w) * (limit ** 2) + wp = np.minimum(w, wmax) + wpp = (wp / wmax) ** exp + depth_prior = np.ravel(wpp[:, None] * np.ones((1, 3))) + return depth_prior diff --git a/mne/minimum_norm/__init__.py b/mne/minimum_norm/__init__.py index 899ec930866..fe3d453e511 100644 --- a/mne/minimum_norm/__init__.py +++ b/mne/minimum_norm/__init__.py @@ -1,3 +1,3 @@ -from .inverse import read_inverse_operator, apply_inverse, minimum_norm, \ - apply_inverse_raw, apply_inverse_epochs +from .inverse import read_inverse_operator, apply_inverse, \ + apply_inverse_raw, make_inverse_operator from .time_frequency import source_induced_power diff --git a/mne/minimum_norm/inverse.py b/mne/minimum_norm/inverse.py index 6930e9465e0..787155af6ab 100644 --- a/mne/minimum_norm/inverse.py +++ b/mne/minimum_norm/inverse.py @@ -1,9 +1,10 @@ # Authors: Alexandre Gramfort # Matti Hamalainen -# Rey Rene Ramirez, Ph.D. # # License: BSD (3-clause) +import warnings +from copy import deepcopy from math import sqrt import numpy as np from scipy import linalg @@ -16,9 +17,9 @@ from ..fiff.tree import dir_tree_find from ..fiff.pick import pick_channels -from ..cov import read_cov +from ..cov import read_cov, prepare_noise_cov +from ..forward import compute_depth_prior from ..source_space import read_source_spaces_from_tree, find_source_space_hemi -from ..forward import _block_diag from ..transforms import invert_transform, transform_source_space_to from ..source_estimate import SourceEstimate @@ -332,7 +333,7 @@ def prepare_inverse_operator(orig, nave, lambda2, dSPM): # Create the projection operator # inv['proj'], ncomp, _ = make_projector(inv['projs'], - inv['noise_cov']['names']) + inv['noise_cov']['names']) if ncomp > 0: print '\tCreated an SSP operator (subspace dimension = %d)' % ncomp @@ -353,7 +354,7 @@ def prepare_inverse_operator(orig, nave, lambda2, dSPM): # inv['whitener'] = np.dot(inv['whitener'], inv['noise_cov']['eigvec']) print ('\tCreated the whitener using a full noise covariance matrix ' - '(%d small eigenvalues omitted)' % (inv['noise_cov']['dim'] + '(%d small eigenvalues omitted)' % (inv['noise_cov']['dim'] - np.sum(nzero))) else: # @@ -732,225 +733,134 @@ def _xyz2lf(Lf_xyz, normals): return Lf_cortex -def minimum_norm(evoked, forward, whitener, method='dspm', - orientation='fixed', snr=3, loose=0.2, depth=True, - weight_exp=0.8, weight_limit=10, fmri=None, fmri_thresh=None, - fmri_off=0.1, pick_normal=False): - """Minimum norm estimate (MNE) +############################################################################### +# Assemble the inverse operator - Compute MNE, dSPM and sLORETA on evoked data starting from - a forward operator. +def make_inverse_operator(info, forward, noise_cov, loose=0.2, depth=0.8): + """Assemble inverse operator Parameters ---------- - evoked: Evoked - Evoked data to invert + info: dict + The measurement info to specify the channels to include. + Bad channels in info['bads'] are ignored. forward: dict Forward operator - whitener: Whitener - Whitening matrix derived from noise covariance matrix - method: 'wmne' | 'dspm' | 'sloreta' - The method to use - orientation: 'fixed' | 'free' | 'loose' - Type of orientation constraints 'fixed'. - snr: float - Signal-to noise ratio defined as in MNE (default: 3). + noise_cov: Covariance + The noise covariance matrix loose: float in [0, 1] Value that weights the source variances of the dipole components defining the tangent space of the cortical surfaces. - depth: bool - Flag to do depth weighting (default: True). - weight_exp: float - Order of the depth weighting. {0=no, 1=full normalization, default=0.8} - weight_limit: float - Maximal amount depth weighting (default: 10). - fmri: array of shape [n_sources] - Vector of fMRI values are the source points. - fmri_thresh: float - fMRI threshold. The source variances of source points with fmri smaller - than fmri_thresh will be multiplied by fmri_off. - fmri_off: float - Weight assigned to non-active source points according to fmri - and fmri_thresh. - pick_normal: bool - If True, rather than pooling the orientations by taking the norm, - only the radial component is kept. This is only implemented - when working with loose orientations. + depth: None | float in [0, 1] + Depth weighting coefficients. If None, no depth weighting is performed. + + # XXX : add support for megreg=0.0, eegreg=0.0 Returns ------- stc: dict Source time courses """ - assert method in ['wmne', 'dspm', 'sloreta'] - assert orientation in ['fixed', 'free', 'loose'] - - if not 0 <= loose <= 1: + is_fixed_ori = (forward['source_ori'] == FIFF.FIFFV_MNE_FIXED_ORI) + if is_fixed_ori and loose is not None: + warnings.warn('Ignoring loose parameter with forward operator with ' + 'fixed orientation.') + if not forward['surf_ori'] and loose is not None: + raise ValueError('Forward operator is not oriented in surface ' + 'coordinates. loose parameter should be None ' + 'not %s.' % loose) + + if loose is not None and not (0 <= loose <= 1): raise ValueError('loose value should be smaller than 1 and bigger than' - ' 0, or empty for no loose orientations.') - if not 0 <= weight_exp <= 1: - raise ValueError('weight_exp should be a scalar between 0 and 1') + ' 0, or None for not loose orientations.') + if depth is not None and not (0 < depth <= 1): + raise ValueError('depth should be a scalar between 0 and 1') - # Set regularization parameter based on SNR - lambda2 = 1.0 / snr ** 2 + fwd_ch_names = [c['ch_name'] for c in forward['chs']] + ch_names = [c['ch_name'] for c in info['chs'] + if (c['ch_name'] not in info['bads']) + and (c['ch_name'] in fwd_ch_names)] + n_chan = len(ch_names) - normals = [] - for s in forward['src']: - normals.append(s['nn'][s['inuse'] != 0]) - normals = np.concatenate(normals) + print "Computing inverse operator with %d channels." % n_chan - W, ch_names = whitener.W, whitener.ch_names + noise_cov = prepare_noise_cov(noise_cov, info, ch_names) + + W = np.zeros((n_chan, n_chan), dtype=np.float) + # + # Omit the zeroes due to projection + # + eig = noise_cov['eig'] + nzero = (eig > 0) + W[nzero, nzero] = 1.0 / np.sqrt(eig[nzero]) + n_nzero = sum(nzero) + # + # Rows of eigvec are the eigenvectors + # + W = np.dot(W, noise_cov['eigvec']) gain = forward['sol']['data'] - fwd_ch_names = [forward['chs'][k]['ch_name'] for k in range(gain.shape[0])] + + n_positions = gain.shape[1] / 3 + fwd_idx = [fwd_ch_names.index(name) for name in ch_names] gain = gain[fwd_idx] - print "Computing inverse solution with %d channels." % len(ch_names) + # Handle depth prior scaling + depth_prior = np.ones(gain.shape[1]) + if depth is not None: + depth_prior = compute_depth_prior(gain, exp=depth) - rank_noise = len(W) - print 'Total rank is %d' % rank_noise - - # processing lead field matrices, weights, areas & orientations - # Initializing. - n_positions = gain.shape[1] / 3 + print "Computing inverse operator with %d channels." % len(ch_names) - if orientation == 'fixed': + if is_fixed_ori: n_dip_per_pos = 1 - elif orientation in ['free', 'loose']: + else: n_dip_per_pos = 3 n_dipoles = n_positions * n_dip_per_pos - w = np.ones(n_dipoles) - - # compute power - if depth: - w = np.sum(gain ** 2, axis=0) - w = w.reshape(-1, 3).sum(axis=1) - w = w[:, None] * np.ones((1, n_dip_per_pos)) - w = w.ravel() - - if orientation == 'fixed': - print 'Appying fixed dipole orientations.' - gain = gain * _block_diag(normals.ravel()[None, :], 3).T - elif orientation == 'free': - print 'Using free dipole orientations. No constraints.' - elif orientation == 'loose': - print 'Transforming lead field matrix to cortical coordinate system.' - gain = _xyz2lf(gain, normals) - # Getting indices for tangential dipoles: [1, 2, 4, 5] - itangential = [k for k in range(n_dipoles) if n_dipoles % 3 != 0] - # Whiten lead field. print 'Whitening lead field matrix.' gain = np.dot(W, gain) - # Computing reciprocal of power. - w = 1.0 / w - - # apply areas - # if ~isempty(areas) - # display('wMNE> Applying areas to compute current source density.') - # areas = areas.^2; - # w = w .* areas; - # end - # clear areas - - # apply depth weighthing - if depth: - # apply weight limit - # Applying weight limit. - print 'Applying weight limit.' - weight_limit2 = weight_limit ** 2 - # limit = min(w(w>min(w) * weight_limit2)); % This is the Matti way. - # we do the Rey way (robust to possible weight discontinuity). - limit = np.min(w) * weight_limit2 - w[w > limit] = limit - - # apply weight exponent - # Applying weight exponent. - print 'Applying weight exponent.' - w = w ** weight_exp - # apply loose orientations - if orientation == 'loose': + orient_prior = np.ones(n_dipoles, dtype=np.float) + if loose is not None: print 'Applying loose dipole orientations. Loose value of %s.' % loose - w[itangential] *= loose + orient_prior[np.mod(np.arange(n_dipoles), 3) != 2] *= loose - # Apply fMRI Priors - if fmri is not None: - print 'Applying fMRI priors.' - w[fmri < fmri_thresh] *= fmri_off + source_cov = orient_prior * depth_prior - # Adjusting Source Covariance matrix to make trace of L*C_J*L' equal + # Adjusting Source Covariance matrix to make trace of G*R*G' equal # to number of sensors. print 'Adjusting source covariance matrix.' - source_std = np.sqrt(w) # sqrt(C_J) - trclcl = linalg.norm(gain * source_std[None, :], ord='fro') - source_std *= sqrt(rank_noise) / trclcl # correct C_J + source_std = np.sqrt(source_cov) gain *= source_std[None, :] + trace_GRGT = linalg.norm(gain, ord='fro') ** 2 + scaling_source_cov = n_nzero / trace_GRGT + source_cov *= scaling_source_cov + gain *= sqrt(scaling_source_cov) - # Compute SVD. - print 'Computing SVD of whitened and weighted lead field matrix.' - U, s, Vh = linalg.svd(gain, full_matrices=False) - ss = s / (s ** 2 + lambda2) - - # Compute whitened MNE operator. - Kernel = source_std[:, None] * np.dot(Vh.T, ss[:, None] * U.T) - - # Compute dSPM operator. - if method == 'dspm': - print 'Computing dSPM inverse operator.' - dspm_diag = np.sum(Kernel ** 2, axis=1) - if n_dip_per_pos == 1: - dspm_diag = np.sqrt(dspm_diag) - elif n_dip_per_pos in [2, 3]: - dspm_diag = dspm_diag.reshape(-1, n_dip_per_pos) - dspm_diag = np.sqrt(np.sum(dspm_diag, axis=1)) - dspm_diag = (np.ones((1, n_dip_per_pos)) * - dspm_diag[:, None]).ravel() - - Kernel /= dspm_diag[:, None] - - # whitened sLORETA imaging kernel - elif method == 'sloreta': - print 'Computing sLORETA inverse operator.' - if n_dip_per_pos == 1: - sloreta_diag = np.sqrt(np.sum(Kernel * gain.T, axis=1)) - Kernel /= sloreta_diag[:, None] - elif n_dip_per_pos in [2, 3]: - for k in n_positions: - start = k * n_dip_per_pos - stop = start + n_dip_per_pos - R = np.dot(Kernel[start:stop, :], gain[:, start:stop]) - SIR = linalg.matfuncs.sqrtm(R, linalg.pinv(R)) - Kernel[start:stop] = np.dot(SIR, Kernel[start:stop]) - - # Multiply inverse operator by whitening matrix, so no need to whiten data - Kernel = np.dot(Kernel, W) - sel = [evoked.ch_names.index(name) for name in ch_names] - sol = np.dot(Kernel, evoked.data[sel]) - - if n_dip_per_pos > 1: - if pick_normal: - print 'Picking only the normal components...', - if orientation != 'loose': - raise ValueError('The pick_normal parameter is only valid ' - 'when working with loose orientations.') - sol = sol[2::3] # take one every 3 sources ie. only the normal - else: - print 'Combining the current components...', - sol = combine_xyz(sol) - - src = forward['src'] - stc = SourceEstimate(None) - stc.data = sol - stc.tmin = float(evoked.first) / evoked.info['sfreq'] - stc.tstep = 1.0 / evoked.info['sfreq'] - stc.lh_vertno = src[0]['vertno'] - stc.rh_vertno = src[1]['vertno'] - stc._init_times() - print '[done]' + # now np.trace(np.dot(gain, gain.T)) == n_nzero + # print np.trace(np.dot(gain, gain.T)), n_nzero - return stc + print 'Computing SVD of whitened and weighted lead field matrix.' + eigen_fields, sing, eigen_leads = linalg.svd(gain, full_matrices=False) + + eigen_fields = dict(data=eigen_fields.T) + eigen_leads = dict(data=eigen_leads.T, nrow=eigen_leads.shape[1]) + depth_prior = dict(data=depth_prior) + orient_prior = dict(data=orient_prior) + source_cov = dict(data=source_cov) + nave = 1.0 + + inv_op = dict(eigen_fields=eigen_fields, eigen_leads=eigen_leads, + sing=sing, nave=nave, depth_prior=depth_prior, + source_cov=source_cov, noise_cov=noise_cov, + orient_prior=orient_prior, projs=deepcopy(info['projs']), + eigen_leads_weighted=False, source_ori=forward['source_ori'], + mri_head_t=deepcopy(forward['mri_head_t']), + src=deepcopy(forward['src'])) + + return inv_op diff --git a/mne/minimum_norm/tests/test_inverse.py b/mne/minimum_norm/tests/test_inverse.py index d40f10f1e84..07fdc6e6e43 100644 --- a/mne/minimum_norm/tests/test_inverse.py +++ b/mne/minimum_norm/tests/test_inverse.py @@ -1,31 +1,33 @@ import os.path as op import numpy as np -# from numpy.testing import assert_array_almost_equal, assert_equal +from numpy.testing import assert_array_almost_equal, assert_equal from ...datasets import sample from ...label import read_label from ...event import read_events from ...epochs import Epochs from ... import fiff, Covariance, read_forward_solution -from ..inverse import minimum_norm, apply_inverse, read_inverse_operator, \ - apply_inverse_raw, apply_inverse_epochs - +from ..inverse import apply_inverse, read_inverse_operator, \ + apply_inverse_raw, apply_inverse_epochs, \ + make_inverse_operator examples_folder = op.join(op.dirname(__file__), '..', '..', '..', 'examples') data_path = sample.data_path(examples_folder) fname_inv = op.join(data_path, 'MEG', 'sample', - 'sample_audvis-meg-oct-6-meg-inv.fif') + # 'sample_audvis-meg-eeg-oct-6-meg-eeg-inv.fif') + 'sample_audvis-meg-oct-6-meg-inv.fif') fname_data = op.join(data_path, 'MEG', 'sample', - 'sample_audvis-ave.fif') + 'sample_audvis-ave.fif') fname_cov = op.join(data_path, 'MEG', 'sample', - 'sample_audvis-cov.fif') + 'sample_audvis-cov.fif') fname_fwd = op.join(data_path, 'MEG', 'sample', - 'sample_audvis-meg-eeg-oct-6-fwd.fif') + 'sample_audvis-meg-oct-6-fwd.fif') + # 'sample_audvis-meg-eeg-oct-6-fwd.fif') fname_raw = op.join(data_path, 'MEG', 'sample', - 'sample_audvis_filt-0-40_raw.fif') + 'sample_audvis_filt-0-40_raw.fif') fname_event = op.join(data_path, 'MEG', 'sample', - 'sample_audvis_filt-0-40_raw-eve.fif') + 'sample_audvis_filt-0-40_raw-eve.fif') label = 'Aud-lh' fname_label = op.join(data_path, 'MEG', 'sample', 'labels', '%s.label' % label) @@ -36,17 +38,29 @@ lambda2 = 1.0 / snr ** 2 dSPM = True +def test_inverse_operator(): + """Test MNE inverse computation -def test_apply_mne_inverse(): - """Test MNE with precomputed inverse operator on Evoked + With and without precomputed inverse operator. """ - setno = 0 evoked = fiff.Evoked(fname_data, setno=0, baseline=(None, 0)) stc = apply_inverse(evoked, inverse_operator, lambda2, dSPM) assert np.all(stc.data > 0) assert np.all(stc.data < 35) + # Test MNE inverse computation starting from forward operator + noise_cov = Covariance(fname_cov) + 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) + + my_stc = apply_inverse(evoked, my_inv_op, lambda2, dSPM) + + assert_equal(stc.times, my_stc.times) + assert_array_almost_equal(stc.data, my_stc.data, 2) + def test_apply_mne_inverse_raw(): """Test MNE with precomputed inverse operator on Raw @@ -77,18 +91,3 @@ def test_apply_mne_inverse_epochs(): assert np.all(stcs[0].data > 0) assert np.all(stcs[0].data < 42) - -def test_compute_minimum_norm(): - """Test MNE inverse computation starting from forward operator - """ - setno = 0 - noise_cov = Covariance(fname_cov) - forward = read_forward_solution(fname_fwd) - evoked = fiff.Evoked(fname_data, setno=setno, baseline=(None, 0)) - whitener = noise_cov.get_whitener(evoked.info, mag_reg=0.1, - grad_reg=0.1, eeg_reg=0.1, pca=True) - stc = minimum_norm(evoked, forward, whitener, - orientation='loose', method='dspm', snr=3, loose=0.2) - - assert np.all(stc.data > 0) - # XXX : test something diff --git a/mne/stats/cluster_level.py b/mne/stats/cluster_level.py index 63dc41aa28d..a12e085fd50 100644 --- a/mne/stats/cluster_level.py +++ b/mne/stats/cluster_level.py @@ -15,7 +15,11 @@ def _get_components(x_in, connectivity): """get connected components from a mask and a connectivity matrix""" - from scikits.learn.utils._csgraph import cs_graph_components + try: + from sklearn.utils._csgraph import cs_graph_components + except: + from scikits.learn.utils._csgraph import cs_graph_components + mask = np.logical_and(x_in[connectivity.row], x_in[connectivity.col]) data = connectivity.data[mask] row = connectivity.row[mask] diff --git a/mne/tests/test_cov.py b/mne/tests/test_cov.py index e3638b84e94..3ec9531d6d4 100644 --- a/mne/tests/test_cov.py +++ b/mne/tests/test_cov.py @@ -6,8 +6,7 @@ from .. import Covariance, read_cov, Epochs, merge_events, \ find_events, write_cov_file, compute_raw_data_covariance, \ compute_covariance -from ..fiff import fiff_open, read_evoked, Raw -from ..datasets import sample +from ..fiff import fiff_open, Raw cov_fname = op.join(op.dirname(__file__), '..', 'fiff', 'tests', 'data', 'test-cov.fif') @@ -73,21 +72,3 @@ def test_cov_estimation_with_triggers(): assert cov_mne.ch_names == cov.ch_names assert (linalg.norm(cov.data - cov_mne.data, ord='fro') / linalg.norm(cov.data, ord='fro')) < 0.06 - - -def test_whitening_cov(): - """Whitening of evoked data and leadfields - """ - data_path = sample.data_path('.') - ave_fname = op.join(data_path, 'MEG', 'sample', - 'sample_audvis-ave.fif') - cov_fname = op.join(data_path, 'MEG', 'sample', - 'sample_audvis-cov.fif') - - # Reading - evoked = read_evoked(ave_fname, setno=0, baseline=(None, 0)) - - cov = Covariance(cov_fname) - cov.get_whitener(evoked.info) - - # XXX : test something