This is an automated email from the git hooks/post-receive script. yoh pushed a commit to annotated tag v0.1 in repository python-mne.
commit b475cbc0ea0bf927df8c48dbc0cd9585f48b7377 Author: Alexandre Gramfort <[email protected]> Date: Fri Jan 14 15:56:09 2011 -0500 ENH : now computing MNE solution --- examples/compute_mne_inverse.py | 21 ++++ examples/read_evoked.py | 2 +- examples/read_raw.py | 2 +- fiff/__init__.py | 2 +- fiff/bem_surfaces.py | 2 +- fiff/cov.py | 6 +- fiff/ctf.py | 2 +- fiff/epochs.py | 8 +- fiff/evoked.py | 6 +- fiff/forward.py | 2 +- fiff/inverse.py | 257 +++++++++++++++++++++++++++++++++++++++- fiff/open.py | 2 +- fiff/pick.py | 77 +++++++++++- fiff/raw.py | 2 +- fiff/source_space.py | 2 +- 15 files changed, 365 insertions(+), 28 deletions(-) diff --git a/examples/compute_mne_inverse.py b/examples/compute_mne_inverse.py new file mode 100644 index 0000000..9673d69 --- /dev/null +++ b/examples/compute_mne_inverse.py @@ -0,0 +1,21 @@ +"""Compute MNE inverse solution +""" +print __doc__ + +import fiff + +fname_inv = 'MNE-sample-data/MEG/sample/sample_audvis-ave-7-meg-inv.fif' +fname_data = 'MNE-sample-data/MEG/sample/sample_audvis-ave.fif' + +# inv = fiff.read_inverse_operator(fname) +setno = 0 +lambda2 = 10 +dSPM = True + +res = fiff.compute_inverse(fname_data, setno, fname_inv, lambda2, dSPM) + +import pylab as pl +pl.plot(res['sol'][::100,:].T) +pl.xlabel('time (s)') +pl.ylabel('Source amplitude') +pl.show() diff --git a/examples/read_evoked.py b/examples/read_evoked.py index 7aec99e..1aae8f2 100644 --- a/examples/read_evoked.py +++ b/examples/read_evoked.py @@ -17,6 +17,6 @@ fiff.write_evoked('evoked.fif', data) import pylab as pl pl.plot(data['evoked']['times'], data['evoked']['epochs'][:306,:].T) -pl.xlabel('time (ms)') +pl.xlabel('time (s)') pl.ylabel('MEG data (T)') pl.show() diff --git a/examples/read_raw.py b/examples/read_raw.py index c078e23..82314dc 100644 --- a/examples/read_raw.py +++ b/examples/read_raw.py @@ -21,6 +21,6 @@ raw.close() # Show MEG data pl.close('all') pl.plot(times, data.T) -pl.xlabel('time (ms)') +pl.xlabel('time (s)') pl.ylabel('MEG data (T)') pl.show() diff --git a/fiff/__init__.py b/fiff/__init__.py index 72ca143..25ddfe2 100644 --- a/fiff/__init__.py +++ b/fiff/__init__.py @@ -10,7 +10,7 @@ from .event import read_events, write_events from .forward import read_forward_solution from .stc import read_stc, write_stc from .bem_surfaces import read_bem_surfaces -from .inverse import read_inverse_operator +from .inverse import read_inverse_operator, compute_inverse from .pick import pick_types from .meas_info import get_current_comp from .epochs import read_epochs diff --git a/fiff/bem_surfaces.py b/fiff/bem_surfaces.py index 1ec7546..a1fdbe8 100644 --- a/fiff/bem_surfaces.py +++ b/fiff/bem_surfaces.py @@ -207,5 +207,5 @@ def _complete_surface_info(this): if size > 0: this['nn'][p,:] = this['nn'][p,:] / size - print '[done]\n' + print '[done]' return this diff --git a/fiff/cov.py b/fiff/cov.py index 36ae8d1..5cb19f6 100644 --- a/fiff/cov.py +++ b/fiff/cov.py @@ -67,7 +67,7 @@ def read_cov(fid, node, cov_kind): # Diagonal is stored data = tag.data diagmat = True - print '\t%d x %d diagonal covariance (kind = %d) found.\n' % (dim, dim, cov_kind) + print '\t%d x %d diagonal covariance (kind = %d) found.' % (dim, dim, cov_kind) else: from scipy import sparse @@ -84,11 +84,11 @@ def read_cov(fid, node, cov_kind): data.flat[::dim+1] /= 2.0 diagmat = False; - print '\t%d x %d full covariance (kind = %d) found.\n' % (dim, dim, cov_kind) + print '\t%d x %d full covariance (kind = %d) found.' % (dim, dim, cov_kind) else: diagmat = False data = tag.data - print '\t%d x %d sparse covariance (kind = %d) found.\n' % (dim, dim, cov_kind) + print '\t%d x %d sparse covariance (kind = %d) found.' % (dim, dim, cov_kind) # Read the possibly precomputed decomposition tag1 = find_tag(fid, this, FIFF.FIFF_MNE_COV_EIGENVALUES) diff --git a/fiff/ctf.py b/fiff/ctf.py index 354d650..fe1c59a 100644 --- a/fiff/ctf.py +++ b/fiff/ctf.py @@ -188,7 +188,7 @@ def read_ctf_comp(fid, node, chs): del col_cals if len(compdata) > 0: - print '\tRead %d compensation matrices\n' % len(compdata) + print '\tRead %d compensation matrices' % len(compdata) return compdata diff --git a/fiff/epochs.py b/fiff/epochs.py index 4928ca9..38bd8f0 100644 --- a/fiff/epochs.py +++ b/fiff/epochs.py @@ -65,7 +65,7 @@ def read_epochs(raw, events, event_id, tmin, tmax, picks=None, for proj in raw['info']['projs']: proj['active'] = True - print '%d projection items activated\n' % len(raw['info']['projs']) + print '%d projection items activated' % len(raw['info']['projs']) # Create the projector proj, nproj = fiff.proj.make_projector_info(raw['info']) @@ -79,7 +79,7 @@ def read_epochs(raw, events, event_id, tmin, tmax, picks=None, # Set up the CTF compensator current_comp = fiff.get_current_comp(raw['info']) if current_comp > 0: - print 'Current compensation grade : %d\n' % current_comp + print 'Current compensation grade : %d' % current_comp if keep_comp: dest_comp = current_comp @@ -98,7 +98,7 @@ def read_epochs(raw, events, event_id, tmin, tmax, picks=None, # raise ValueError, 'Raw file name does not end properly' # # events = fiff.read_events(event_name) - # print 'Events read from %s\n' % event_name + # print 'Events read from %s' % event_name # else: # # Binary file # if event_name.endswith('-eve.fif'): @@ -121,7 +121,7 @@ def read_epochs(raw, events, event_id, tmin, tmax, picks=None, # raise ValueError, ('This new format event file is not compatible' # ' with the raw data') # else: - # print 'The text event file %s is in the old format\n' % event_name + # print 'The text event file %s is in the old format' % event_name # # Offset with first sample # events[:,0] += raw['first_samp'] diff --git a/fiff/evoked.py b/fiff/evoked.py index 5038f2d..eaf7ab4 100644 --- a/fiff/evoked.py +++ b/fiff/evoked.py @@ -28,7 +28,7 @@ def read_evoked(fname, setno=0): if setno < 0: raise ValueError, 'Data set selector must be positive' - print 'Reading %s ...\n' % fname + print 'Reading %s ...' % fname fid, tree, _ = fiff_open(fname) # Read the measurement info @@ -73,7 +73,7 @@ def read_evoked(fname, setno=0): naspect += nsaspects print '\t%d evoked data sets containing a total of %d data aspects' \ - ' in %s\n' % (len(evoked), naspect, fname) + ' in %s' % (len(evoked), naspect, fname) if setno >= naspect or setno < 0: fid.close() @@ -174,7 +174,7 @@ def read_evoked(fname, setno=0): tag = read_tag(fid, pos) epoch.append(tag) - print '\t\tnave = %d aspect type = %d\n' % (nave, aspect_kind) + print '\t\tnave = %d aspect type = %d' % (nave, aspect_kind) nepoch = len(epoch) if nepoch != 1 and nepoch != info.nchan: diff --git a/fiff/forward.py b/fiff/forward.py index 87a014c..6a7b2cc 100644 --- a/fiff/forward.py +++ b/fiff/forward.py @@ -447,7 +447,7 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, if nuse == 0: raise ValueError, 'Nothing remains after picking' - print '\t%d out of %d channels remain after picking\n' % (nuse, + print '\t%d out of %d channels remain after picking' % (nuse, fwd['nchan']) # Pick the correct rows of the forward operator diff --git a/fiff/inverse.py b/fiff/inverse.py index ed4ab45..b18c67d 100644 --- a/fiff/inverse.py +++ b/fiff/inverse.py @@ -1,13 +1,18 @@ +from math import sqrt +import numpy as np + from .constants import FIFF from .open import fiff_open from .tag import find_tag from .matrix import _read_named_matrix, _transpose_named_matrix from .cov import read_cov -from .proj import read_proj +from .proj import read_proj, make_projector from .tree import dir_tree_find from .source_space import read_source_spaces, find_source_space_hemi from .forward import _invert_transform, _transform_source_space_to - +from .evoked import read_evoked +from .pick import pick_channels_evoked +from .forward import _block_diag def read_inverse_operator(fname): """Read the inverse operator decomposition from a FIF file @@ -89,7 +94,7 @@ def read_inverse_operator(fname): raise ValueError, 'Source orientation information not found' inv['source_nn'] = tag.data - print '[done]\n' + print '[done]' # # The SVD decomposition... # @@ -151,13 +156,13 @@ def read_inverse_operator(fname): try: inv['depth_prior'] = read_cov(fid, invs, FIFF.FIFFV_MNE_DEPTH_PRIOR_COV) - print '\tDepth priors read.\n' + print '\tDepth priors read.' except: inv['depth_prior'] = []; try: inv['fmri_prior'] = read_cov(fid, invs, FIFF.FIFFV_MNE_FMRI_PRIOR_COV) - print '\tfMRI priors read.\n' + print '\tfMRI priors read.' except: inv['fmri_prior'] = []; @@ -235,4 +240,244 @@ def read_inverse_operator(fname): # fid.close() - return inv \ No newline at end of file + return inv + +############################################################################### +# Compute inverse solution + +def combine_xyz(vec): + """ + % + % function [comb] = mne_combine_xyz(vec) + % + % Compute the three Cartesian components of a vector together + % + % + % vec - Input row or column vector [ x1 y1 z1 ... x_n y_n z_n ] + % comb - Output vector [x1^2+y1^2+z1^2 ... x_n^2+y_n^2+z_n^2 ] + % + """ + if vec.ndim != 1 or (vec.size % 3) != 0: + raise ValueError, ('Input must be a 1D vector with ' + '3N entries') + + s = _block_diag(vec[None,:], 3) + comb = (s * s.T).diagonal() + return comb + + +def prepare_inverse_operator(orig, nave, lambda2, dSPM): + """ + % + % [inv] = mne_prepare_inverse_operator(orig,nave,lambda2,dSPM) + % + % Prepare for actually computing the inverse + % + % orig - The inverse operator structure read from a file + % nave - Number of averages (scales the noise covariance) + % lambda2 - The regularization factor + % dSPM - Compute the noise-normalization factors for dSPM? + % + """ + + if nave <= 0: + raise ValueError, 'The number of averages should be positive' + + print 'Preparing the inverse operator for use...' + inv = orig.copy() + # + # Scale some of the stuff + # + scale = float(inv['nave']) / nave + inv['noise_cov']['data'] = scale * inv['noise_cov']['data'] + inv['noise_cov']['eig'] = scale * inv['noise_cov']['eig'] + inv['source_cov']['data'] = scale * inv['source_cov']['data'] + # + if inv['eigen_leads_weighted']: + inv['eigen_leads']['data'] = sqrt(scale) * inv['eigen_leads']['data'] + + + print ('\tScaled noise and source covariance from nave = %d to ' + 'nave = %d' % (inv['nave'], nave)) + inv['nave'] = nave + # + # Create the diagonal matrix for computing the regularized inverse + # + inv['reginv'] = inv['sing'] / (inv['sing'] * inv['sing'] + lambda2); + print '\tCreated the regularized inverter' + # + # Create the projection operator + # + inv['proj'], ncomp, _ = make_projector(inv['projs'], + inv['noise_cov']['names']) + if ncomp > 0: + print '\tCreated an SSP operator (subspace dimension = %d)' % ncomp + + # + # Create the whitener + # + inv['whitener'] = np.zeros((inv['noise_cov']['dim'], + inv['noise_cov']['dim'])) + if inv['noise_cov']['diag'] == 0: + # + # Omit the zeroes due to projection + # + nnzero = 0 + for k in range(ncomp, inv['noise_cov']['dim']): + if inv['noise_cov']['eig'][k] > 0: + inv['whitener'][k,k] = 1.0 / sqrt(inv['noise_cov']['eig'][k]) + nnzero += 1 + + # + # Rows of eigvec are the eigenvectors + # + 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'] - nnzero)) + else: + # + # No need to omit the zeroes due to projection + # + for k in range(inv['noise_cov']['dim']): + inv['whitener'][k,k] = 1.0 / sqrt(inv['noise_cov']['data'][k]) + + print ('\tCreated the whitener using a diagonal noise covariance ' + 'matrix (%d small eigenvalues discarded)' % ncomp) + + # + # Finally, compute the noise-normalization factors + # + if dSPM: + print '\tComputing noise-normalization factors...' + noise_norm = np.zeros(inv['eigen_leads']['nrow']) + if inv['eigen_leads_weighted']: + for k in range(inv['eigen_leads']['nrow']): + one = inv['eigen_leads']['data'][k,:] * inv['reginv'] + noise_norm[k] = np.sum(one**2) + else: + for k in range(inv['eigen_leads']['nrow']): + one = sqrt(inv['source_cov']['data'][k]) * \ + np.sum(inv['eigen_leads']['data'][k,:] * inv['reginv']) + noise_norm[k] = np.sum(one**2) + + # + # Compute the final result + # + if inv['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI: + # + # The three-component case is a little bit more involved + # The variances at three consequtive entries must be squeared and + # added together + # + # Even in this case return only one noise-normalization factor + # per source location + # + noise_norm = np.sqrt(combine_xyz(noise_norm)) + # + # This would replicate the same value on three consequtive + # entries + # + # noise_norm = kron(sqrt(mne_combine_xyz(noise_norm)),ones(3,1)); + + inv['noisenorm'] = np.diag(1.0 / np.abs(noise_norm)) # XXX + print '[done]' + else: + inv['noisenorm'] = []; + + return inv + + +def compute_inverse(fname_data, setno, fname_inv, lambda2, dSPM=True, + nave=None): + """ + % + % [res] = mne_ex_compute_inverse(fname_data,setno,fname_inv,nave,lambda2,dSPM) + % + % An example on how to compute a L2-norm inverse solution + % Actual code using these principles might be different because + % the inverse operator is often reused across data sets. + % + % + % fname_data - Name of the data file + % setno - Data set number + % fname_inv - Inverse operator file name + % nave - Number of averages (scales the noise covariance) + % If negative, the number of averages in the data will be + % used + % lambda2 - The regularization factor + % dSPM - do dSPM? + % + """ + + # + # Read the data first + # + data = read_evoked(fname_data, setno) + # + # Then the inverse operator + # + inv = read_inverse_operator(fname_inv) + # + # Set up the inverse according to the parameters + # + if nave is None: + nave = data['evoked']['nave'] + + inv = prepare_inverse_operator(inv, nave, lambda2, dSPM) + # + # Pick the correct channels from the data + # + data = pick_channels_evoked(data, inv['noise_cov']['names']) + print 'Picked %d channels from the data' % data['info']['nchan'] + print 'Computing inverse...', + # + # Simple matrix multiplication followed by combination of the + # three current components + # + # This does all the data transformations to compute the weights for the + # eigenleads + # + trans = reduce(np.dot, [np.diag(inv['reginv']), + inv['eigen_fields']['data'], + inv['whitener'], + inv['proj'], + data['evoked']['epochs']]) + # + # Transformation into current distributions by weighting the eigenleads + # with the weights computed above + # + if inv['eigen_leads_weighted']: + # + # R^0.5 has been already factored in + # + print '(eigenleads already weighted)...', + sol = np.dot(inv['eigen_leads']['data'], trans) + else: + # + # R^0.5 has to factored in + # + print '(eigenleads need to be weighted)...', + sol = np.sqrt(inv['source_cov']['data'])[:,None] * \ + np.dot(inv['eigen_leads']['data'], trans) + + + if inv['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI: + print 'combining the current components...', + sol1 = np.zeros((sol.shape[0]/3, sol.shape[1])) + for k in range(sol.shape[1]): + sol1[:,k] = np.sqrt(combine_xyz(sol[:,k])) + + sol = sol1 + + if dSPM: + print '(dSPM)...', + sol = np.dot(inv['noisenorm'], sol) + + res = dict() + res['inv'] = inv + res['sol'] = sol + res['tmin'] = float(data['evoked']['first']) / data['info']['sfreq'] + res['tstep'] = 1.0 / data['info']['sfreq'] + print '[done]' + + return res diff --git a/fiff/open.py b/fiff/open.py index 93c2db0..7cb57d7 100644 --- a/fiff/open.py +++ b/fiff/open.py @@ -67,7 +67,7 @@ def fiff_open(fname, verbose=False): tree, _ = make_dir_tree(fid, directory, verbose=verbose) if verbose: - print '[done]\n' + print '[done]' # Back to the beginning fid.seek(0) diff --git a/fiff/pick.py b/fiff/pick.py index c6bf36e..31dd241 100644 --- a/fiff/pick.py +++ b/fiff/pick.py @@ -1,3 +1,5 @@ +from copy import copy + import numpy as np from .constants import FIFF @@ -13,10 +15,10 @@ def pick_channels(ch_names, include, exclude): List of channels include : list of string - List of channels to include. If None include all available. + List of channels to include. If empty include all available. exclude : list of string - List of channels to exclude. If None do not exclude any channel. + List of channels to exclude. If empty do not exclude any channel. Returns ------- @@ -25,7 +27,7 @@ def pick_channels(ch_names, include, exclude): """ sel = [] for k, name in enumerate(ch_names): - if (include is None or name in include) and name not in exclude: + if (include is [] or name in include) and name not in exclude: sel.append(k) sel = np.unique(sel) np.sort(sel) @@ -82,3 +84,72 @@ def pick_types(info, meg=True, eeg=False, stim=False, include=[], exclude=[]): sel = pick_channels(info['ch_names'], myinclude, exclude) return sel + + +def pick_info(info, sel=[]): + """Restrict an info structure to a selection of channels + + Parameters + ---------- + info : dict + Info structure from evoked or raw data + sel : list of int + Indices of channels to include + + Returns + ------- + res : dict + Info structure restricted to a selection of channels + """ + + res = copy(info) + if len(sel) == 0: + raise ValueError, 'Warning : No channels match the selection.' + + res['chs'] = [res['chs'][k] for k in sel] + res['ch_names'] = [res['ch_names'][k] for k in sel] + res['nchan'] = len(sel) + return res + + +def pick_channels_evoked(orig, include=[], exclude=[]): + """Pick channels from evoked data + + Parameters + ---------- + orig : dict + One evoked data + + include : list of string, (optional) + List of channels to include. (if None, include all available) + + exclude : list of string, (optional) + Channels to exclude (if None, do not exclude any) + + Returns + ------- + res : dict + Evoked data restricted to selected channels. If include and + exclude are None it returns orig without copy. + """ + + if include is None and exclude is None: + return orig + + sel = pick_channels(orig['info']['ch_names'], include=include, + exclude=exclude) + + if len(sel) == 0: + raise ValueError, 'Warning : No channels match the selection.' + + res = orig.copy() + # + # Modify the measurement info + # + res['info'] = pick_info(res['info'], sel) + # + # Create the reduced data set + # + res['evoked']['epochs'] = res['evoked']['epochs'][sel,:] + + return res diff --git a/fiff/raw.py b/fiff/raw.py index a529f90..d9eb0c9 100644 --- a/fiff/raw.py +++ b/fiff/raw.py @@ -164,7 +164,7 @@ def setup_read_raw(fname, allow_maxshield=False): data['first_samp'], data['last_samp'], float(data['first_samp']) / data['info']['sfreq'], float(data['last_samp']) / data['info']['sfreq']) - print 'Ready.\n' + print 'Ready.' return data diff --git a/fiff/source_space.py b/fiff/source_space.py index 9381250..fc8daaa 100644 --- a/fiff/source_space.py +++ b/fiff/source_space.py @@ -85,7 +85,7 @@ def read_source_spaces(source, add_geom=False, tree=None): src.append(this) - print '\t%d source spaces read\n' % len(spaces) + print '\t%d source spaces read' % len(spaces) if open_here: fid.close() -- 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
