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 9c3e90b61bc84de19b9f9c125444ac4640f703c4 Author: Alexandre Gramfort <[email protected]> Date: Tue Jan 18 18:34:17 2011 -0500 examples with new sample data + some docstrings --- examples/compute_mne_inverse.py | 9 ++- examples/read_bem_surfaces.py | 22 ++++-- examples/read_forward.py | 3 +- examples/read_inverse.py | 26 +++---- examples/read_stc.py | 2 +- mne/inverse.py | 151 ++++++++++++++++++++++------------------ mne/source_space.py | 2 +- 7 files changed, 119 insertions(+), 96 deletions(-) diff --git a/examples/compute_mne_inverse.py b/examples/compute_mne_inverse.py index 94efae0..ed97060 100644 --- a/examples/compute_mne_inverse.py +++ b/examples/compute_mne_inverse.py @@ -5,20 +5,19 @@ Compute MNE inverse solution """ print __doc__ -from mne import fiff +import mne -fname_inv = 'MNE-sample-data/MEG/sample/sample_audvis-ave-7-meg-inv.fif' +fname_inv = 'MNE-sample-data/MEG/sample/sample_audvis-meg-oct-6-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) +res = mne.compute_inverse(fname_data, setno, fname_inv, lambda2, dSPM) import pylab as pl pl.plot(res['sol'][::100,:].T) -pl.xlabel('time (s)') +pl.xlabel('time (ms)') pl.ylabel('Source amplitude') pl.show() diff --git a/examples/read_bem_surfaces.py b/examples/read_bem_surfaces.py index 127f9c4..c4faa5d 100644 --- a/examples/read_bem_surfaces.py +++ b/examples/read_bem_surfaces.py @@ -5,19 +5,27 @@ Reading BEM surfaces from a forward solution """ print __doc__ -from mne import fiff +import mne -fname = 'MNE-sample-data/subjects/sample/bem/sample-5120-bem-sol.fif' +fname = 'MNE-sample-data/subjects/sample/bem/sample-5120-5120-5120-bem-sol.fif' -surf = fiff.read_bem_surfaces(fname) +surfaces = mne.read_bem_surfaces(fname) -print "Number of surfaces : %d" % len(surf) +print "Number of surfaces : %d" % len(surfaces) ############################################################################### # Show result +head_col = (0.95, 0.83, 0.83) # light pink +skull_col = (0.91, 0.89, 0.67) +brain_col = (0.67, 0.89, 0.91) # light blue +colors = [head_col, skull_col, brain_col] + # 3D source space -points = surf[0]['rr'] -faces = surf[0]['tris'] from enthought.mayavi import mlab -mlab.triangular_mesh(points[:,0], points[:,1], points[:,2], faces) +mlab.clf() +for c, surf in zip(colors, surfaces): + points = surf['rr'] + faces = surf['tris'] + mlab.triangular_mesh(points[:,0], points[:,1], points[:,2], faces, + color=c, opacity=0.3) diff --git a/examples/read_forward.py b/examples/read_forward.py index 75a9bad..c6377f0 100644 --- a/examples/read_forward.py +++ b/examples/read_forward.py @@ -7,8 +7,7 @@ print __doc__ import mne -# fname = 'MNE-sample-data/MEG/sample/sample_audvis-ave-7-fwd.fif' -fname = 'sm01a5-ave-oct-6-fwd.fif' +fname = 'MNE-sample-data/MEG/sample/sample_audvis-meg-oct-6-fwd.fif' data = mne.read_forward_solution(fname) leadfield = data['sol']['data'] diff --git a/examples/read_inverse.py b/examples/read_inverse.py index 96e8c47..00025ed 100644 --- a/examples/read_inverse.py +++ b/examples/read_inverse.py @@ -7,7 +7,7 @@ print __doc__ import mne -fname = 'MNE-sample-data/MEG/sample/sample_audvis-ave-7-meg-inv.fif' +fname = 'MNE-sample-data/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif' inv = mne.read_inverse_operator(fname) @@ -16,14 +16,16 @@ print "fMRI prior: %s" % inv['fmri_prior'] print "Number of sources: %s" % inv['nsource'] print "Number of channels: %s" % inv['nchan'] -# ############################################################################### -# # Show result -# -# # 3D source space -# lh_points = inv['src'][0]['rr'] -# lh_faces = inv['src'][0]['use_tris'] -# rh_points = inv['src'][1]['rr'] -# rh_faces = inv['src'][1]['use_tris'] -# from enthought.mayavi import mlab -# mlab.triangular_mesh(lh_points[:,0], lh_points[:,1], lh_points[:,2], lh_faces) -# mlab.triangular_mesh(rh_points[:,0], rh_points[:,1], rh_points[:,2], rh_faces) +############################################################################### +# Show result + +# 3D source space +import numpy as np +lh_points = inv['src'][0]['rr'] +lh_faces = inv['src'][0]['use_tris'] +rh_points = inv['src'][1]['rr'] +rh_faces = inv['src'][1]['use_tris'] +from enthought.mayavi import mlab +mlab.triangular_mesh(lh_points[:,0], lh_points[:,1], lh_points[:,2], lh_faces) +mlab.triangular_mesh(rh_points[:,0], rh_points[:,1], rh_points[:,2], rh_faces) +mlab.show() diff --git a/examples/read_stc.py b/examples/read_stc.py index 17cd649..d5acd0e 100644 --- a/examples/read_stc.py +++ b/examples/read_stc.py @@ -10,7 +10,7 @@ print __doc__ import mne -fname = 'MNE-sample-data/MEG/sample/sample_audvis-ave-7-meg-lh.stc' +fname = 'MNE-sample-data/MEG/sample/sample_audvis-meg-lh.stc' stc = mne.read_stc(fname) diff --git a/mne/inverse.py b/mne/inverse.py index 8e67bae..53d92b7 100644 --- a/mne/inverse.py +++ b/mne/inverse.py @@ -66,14 +66,14 @@ def read_inverse_operator(fname): raise ValueError, 'Modalities not found' inv = dict() - inv['methods'] = tag.data; + inv['methods'] = tag.data tag = find_tag(fid, invs, FIFF.FIFF_MNE_SOURCE_ORIENTATION) if tag is None: fid.close() raise ValueError, 'Source orientation constraints not found' - inv['source_ori'] = tag.data; + inv['source_ori'] = tag.data tag = find_tag(fid, invs, FIFF.FIFF_MNE_SOURCE_SPACE_NPOINTS) if tag is None: @@ -110,26 +110,29 @@ def read_inverse_operator(fname): fid.close() raise ValueError, 'Singular values not found' - inv['sing'] = tag.data + inv['sing'] = tag.data inv['nchan'] = len(inv['sing']) # # The eigenleads and eigenfields # inv['eigen_leads_weighted'] = False try: - inv['eigen_leads'] = _read_named_matrix(fid, invs, FIFF.FIFF_MNE_INVERSE_LEADS) + inv['eigen_leads'] = _read_named_matrix(fid, invs, + FIFF.FIFF_MNE_INVERSE_LEADS) except: - inv['eigen_leads_weighted'] = True - try: - inv.eigen_leads = _read_named_matrix(fid,invs,FIFF.FIFF_MNE_INVERSE_LEADS_WEIGHTED); - except Exception as inst: - raise ValueError, '%s' % inst + inv['eigen_leads_weighted'] = True + try: + inv.eigen_leads = _read_named_matrix(fid, invs, + FIFF.FIFF_MNE_INVERSE_LEADS_WEIGHTED) + except Exception as inst: + raise ValueError, '%s' % inst # # Having the eigenleads as columns is better for the inverse calculations # inv['eigen_leads'] = _transpose_named_matrix(inv['eigen_leads']) try: - inv['eigen_fields'] = _read_named_matrix(fid, invs, FIFF.FIFF_MNE_INVERSE_FIELDS) + inv['eigen_fields'] = _read_named_matrix(fid, invs, + FIFF.FIFF_MNE_INVERSE_FIELDS) except Exception as inst: raise ValueError, '%s' % inst @@ -164,13 +167,13 @@ def read_inverse_operator(fname): FIFF.FIFFV_MNE_DEPTH_PRIOR_COV) print '\tDepth priors read.' except: - inv['depth_prior'] = []; + inv['depth_prior'] = [] try: inv['fmri_prior'] = read_cov(fid, invs, FIFF.FIFFV_MNE_FMRI_PRIOR_COV) print '\tfMRI priors read.' except: - inv['fmri_prior'] = []; + inv['fmri_prior'] = [] # # Read the source spaces @@ -182,7 +185,7 @@ def read_inverse_operator(fname): raise ValueError, 'Could not read the source spaces (%s)' % inst for s in inv['src']: - s['id'] = find_source_space_hemi(s) + s['id'] = find_source_space_hemi(s) # # Get the MRI <-> head coordinate transformation @@ -192,14 +195,15 @@ def read_inverse_operator(fname): fid.close() raise ValueError, 'MRI/head coordinate transformation not found' else: - mri_head_t = tag.data; + mri_head_t = tag.data if mri_head_t['from_'] != FIFF.FIFFV_COORD_MRI or \ mri_head_t['to'] != FIFF.FIFFV_COORD_HEAD: mri_head_t = _invert_transform(mri_head_t) if mri_head_t['from_'] != FIFF.FIFFV_COORD_MRI or \ mri_head_t['to'] != FIFF.FIFFV_COORD_HEAD: fid.close() - raise ValueError, 'MRI/head coordinate transformation not found' + raise ValueError, ('MRI/head coordinate transformation ' + 'not found') inv['mri_head_t'] = mri_head_t # @@ -215,7 +219,7 @@ def read_inverse_operator(fname): # # Number of averages is initially one # - inv['nave'] = 1; + inv['nave'] = 1 # # We also need the SSP operator # @@ -223,24 +227,25 @@ def read_inverse_operator(fname): # # Some empty fields to be filled in later # - inv['proj'] = [] # This is the projector to apply to the data - inv['whitener'] = [] # This whitens the data - inv['reginv'] = [] # This the diagonal matrix implementing + inv['proj'] = [] # This is the projector to apply to the data + inv['whitener'] = [] # This whitens the data + inv['reginv'] = [] # This the diagonal matrix implementing # regularization and the inverse inv['noisenorm'] = [] # These are the noise-normalization factors # nuse = 0 for k in range(len(inv['src'])): - try: - inv['src'][k] = _transform_source_space_to(inv['src'][k], + try: + inv['src'][k] = _transform_source_space_to(inv['src'][k], inv['coord_frame'], mri_head_t) - except Exception as inst: - fid.close() - raise ValueError, 'Could not transform source space (%s)', inst + except Exception as inst: + fid.close() + raise ValueError, 'Could not transform source space (%s)', inst - nuse += inv['src'][k]['nuse'] + nuse += inv['src'][k]['nuse'] - print '\tSource spaces transformed to the inverse solution coordinate frame' + print ('\tSource spaces transformed to the inverse solution ' + 'coordinate frame') # # Done! # @@ -267,7 +272,7 @@ def combine_xyz(vec): raise ValueError, ('Input must be a 1D vector with ' '3N entries') - s = _block_diag(vec[None,:], 3) + s = _block_diag(vec[None, :], 3) comb = (s * s.T).diagonal() return comb @@ -295,12 +300,12 @@ def prepare_inverse_operator(orig, nave, lambda2, dSPM): # 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['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'] + inv['eigen_leads']['data'] = sqrt(scale) * inv['eigen_leads']['data'] print ('\tScaled noise and source covariance from nave = %d to ' @@ -309,7 +314,7 @@ def prepare_inverse_operator(orig, nave, lambda2, dSPM): # # Create the diagonal matrix for computing the regularized inverse # - inv['reginv'] = inv['sing'] / (inv['sing'] * inv['sing'] + lambda2); + inv['reginv'] = inv['sing'] / (inv['sing'] * inv['sing'] + lambda2) print '\tCreated the regularized inverter' # # Create the projection operator @@ -330,22 +335,23 @@ def prepare_inverse_operator(orig, nave, lambda2, dSPM): # 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 + 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)) + '(%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]) + 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) @@ -358,12 +364,13 @@ def prepare_inverse_operator(orig, nave, lambda2, dSPM): 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'] + 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']) + np.sum(inv['eigen_leads']['data'][k, :] + * inv['reginv']) noise_norm[k] = np.sum(one**2) # @@ -388,31 +395,39 @@ def prepare_inverse_operator(orig, nave, lambda2, dSPM): inv['noisenorm'] = np.diag(1.0 / np.abs(noise_norm)) # XXX print '[done]' else: - inv['noisenorm'] = []; + 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? - % + """Compute inverse solution + + Computes a L2-norm inverse solution + Actual code using these principles might be different because + the inverse operator is often reused across data sets. + + Parameters + ---------- + fname: string + File name of the data file + setno: int + Data set number + fname_inv: string + File name of the inverse operator + nave: int + Number of averages (scales the noise covariance) + If negative, the number of averages in the data will be used XXX + lambda2: float + The regularization parameter + dSPM: bool + do dSPM ? + + Returns + ------- + res: dict + Inverse solution """ # @@ -453,17 +468,17 @@ def compute_inverse(fname_data, setno, fname_inv, lambda2, dSPM=True, # 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) + # + # 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] * \ + # + # 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) @@ -471,7 +486,7 @@ def compute_inverse(fname_data, setno, fname_inv, lambda2, dSPM=True, 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])) + sol1[:, k] = np.sqrt(combine_xyz(sol[:, k])) sol = sol1 diff --git a/mne/source_space.py b/mne/source_space.py index 1b10c8d..074c82a 100644 --- a/mne/source_space.py +++ b/mne/source_space.py @@ -144,7 +144,7 @@ def _read_one_source_space(fid, this, open_here): fid.close() raise ValueError, 'Vertex data not found' - res['rr'] = tag.data + res['rr'] = tag.data.astype(np.float) # make it double precision for mayavi if res['rr'].shape[0] != res['np']: 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
