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 7786a0fd31b695fcce3900a5c0178cdadfa6efb3 Author: Alexandre Gramfort <[email protected]> Date: Mon Jan 24 14:14:35 2011 -0500 ENH: make MNE-dSPM computation work --- examples/plot_compute_mne_inverse.py | 26 +++++++++++++++----------- mne/cov.py | 8 +------- mne/fiff/matrix.py | 6 ++---- mne/inverse.py | 26 ++++++++++++-------------- mne/stc.py | 29 +++++++++++++++-------------- 5 files changed, 45 insertions(+), 50 deletions(-) diff --git a/examples/plot_compute_mne_inverse.py b/examples/plot_compute_mne_inverse.py index 99beefc..8cf7f7f 100644 --- a/examples/plot_compute_mne_inverse.py +++ b/examples/plot_compute_mne_inverse.py @@ -1,7 +1,11 @@ """ -============================ -Compute MNE inverse solution -============================ +================================= +Compute MNE-dSPM inverse solution +================================= + +Compute dSPM inverse solution on MNE sample data set +and stores the solution in stc files for visualisation. + """ # Author: Alexandre Gramfort <[email protected]> @@ -26,16 +30,16 @@ dSPM = True res = mne.compute_inverse(fname_data, setno, fname_inv, lambda2, dSPM, baseline=(None, 0)) -# XXX : kind of ugly -import numpy as np -res['vertices'] = np.r_[res['inv']['src'][0]['vertno']] -# res['vertices'] = np.r_[res['inv']['src'][0]['vertno'], -# res['inv']['src'][1]['vertno']] -# res['data'] = res['sol'] -res['data'] = res['sol'][:len(res['vertices'])] +lh_vertices = res['inv']['src'][0]['vertno'] +rh_vertices = res['inv']['src'][1]['vertno'] +lh_data = res['sol'][:len(lh_vertices)] +rh_data = res['sol'][len(rh_vertices):] # Save result in stc file -mne.write_stc('mne_dSPM_inverse-lh.stc', res) +mne.write_stc('mne_dSPM_inverse-lh.stc', tmin=res['tmin'], tstep=res['tstep'], + vertices=lh_vertices, data=lh_data) +mne.write_stc('mne_dSPM_inverse-rh.stc', tmin=res['tmin'], tstep=res['tstep'], + vertices=rh_vertices, data=rh_data) import pylab as pl pl.plot(res['sol'][::100,:].T) diff --git a/mne/cov.py b/mne/cov.py index a35b427..16882ef 100644 --- a/mne/cov.py +++ b/mne/cov.py @@ -84,14 +84,9 @@ def read_cov(fid, node, cov_kind): # Lower diagonal is stored vals = tag.data data = np.zeros((dim, dim)) - q = 0 - for j in range(dim): - for k in range(j): - data[j, k] = vals[q]; - q += 1 + data[np.tril(np.ones((dim, dim))) > 0] = vals data = data + data.T data.flat[::dim+1] /= 2.0 - diagmat = False; print '\t%d x %d full covariance (kind = %d) found.' % (dim, dim, cov_kind) else: @@ -109,7 +104,6 @@ def read_cov(fid, node, cov_kind): eig = None eigvec = None - # Read the projection operator projs = read_proj(fid, this) diff --git a/mne/fiff/matrix.py b/mne/fiff/matrix.py index 199f7a3..595a432 100644 --- a/mne/fiff/matrix.py +++ b/mne/fiff/matrix.py @@ -10,10 +10,8 @@ from .tag import find_tag, has_tag def _transpose_named_matrix(mat): """Transpose mat inplace (no copy) """ - mat['nrow'] = mat['ncol'] - mat['ncol'] = mat['nrow'] - mat['row_names'] = mat['col_names'] - mat['col_names'] = mat['row_names'] + mat['nrow'], mat['ncol'] = mat['ncol'], mat['nrow'] + mat['row_names'], mat['col_names'] = mat['col_names'], mat['row_names'] mat['data'] = mat['data'].T return mat diff --git a/mne/inverse.py b/mne/inverse.py index 41331d8..197e6fc 100644 --- a/mne/inverse.py +++ b/mne/inverse.py @@ -123,7 +123,7 @@ def read_inverse_operator(fname): except: inv['eigen_leads_weighted'] = True try: - inv.eigen_leads = _read_named_matrix(fid, invs, + inv['eigen_leads'] = _read_named_matrix(fid, invs, FIFF.FIFF_MNE_INVERSE_LEADS_WEIGHTED) except Exception as inst: raise ValueError, '%s' % inst @@ -158,7 +158,7 @@ def read_inverse_operator(fname): # Read the various priors # try: - inv.orient_prior = read_cov(fid, invs, FIFF.FIFFV_MNE_ORIENT_PRIOR_COV) + inv['orient_prior'] = read_cov(fid, invs, FIFF.FIFFV_MNE_ORIENT_PRIOR_COV) print '\tOrientation priors read.' except Exception as inst: inv['orient_prior'] = [] @@ -366,13 +366,12 @@ def prepare_inverse_operator(orig, nave, lambda2, dSPM): 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) + noise_norm[k] = sqrt(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) + inv['eigen_leads']['data'][k, :] * inv['reginv'] + noise_norm[k] = sqrt(np.sum(one**2)) # # Compute the final result @@ -393,7 +392,7 @@ def prepare_inverse_operator(orig, nave, lambda2, dSPM): # # noise_norm = kron(sqrt(mne_combine_xyz(noise_norm)),ones(3,1)); - inv['noisenorm'] = np.diag(1.0 / np.abs(noise_norm)) # XXX + inv['noisenorm'] = 1.0 / np.abs(noise_norm) print '[done]' else: inv['noisenorm'] = [] @@ -467,11 +466,11 @@ def compute_inverse(fname_data, setno, fname_inv, lambda2, dSPM=True, # 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']]) + trans = inv['reginv'][:,None] * reduce(np.dot, [inv['eigen_fields']['data'], + inv['whitener'], + inv['proj'], + data['evoked']['epochs']]) + # # Transformation into current distributions by weighting the eigenleads # with the weights computed above @@ -490,7 +489,6 @@ def compute_inverse(fname_data, setno, fname_inv, lambda2, dSPM=True, 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])) @@ -501,7 +499,7 @@ def compute_inverse(fname_data, setno, fname_inv, lambda2, dSPM=True, if dSPM: print '(dSPM)...', - sol = np.dot(inv['noisenorm'], sol) + sol *= inv['noisenorm'][:, None] res = dict() res['inv'] = inv diff --git a/mne/stc.py b/mne/stc.py index fdf9be8..9184ac2 100644 --- a/mne/stc.py +++ b/mne/stc.py @@ -62,38 +62,39 @@ def read_stc(filename): return stc -def write_stc(filename, stc): +def write_stc(filename, tmin, tstep, vertices, data): """Write an STC file Parameters ---------- filename: string The name of the STC file - - stc: dict - The STC structure. It has the following keys: - tmin The first time point of the data in seconds - tstep Time between frames in seconds - vertices vertex indices (0 based) - data The data matrix (nvert * ntime) + tmin: float + The first time point of the data in seconds + tstep: float + Time between frames in seconds + vertices: array of integers + Vertex indices (0 based) + data: 2D array + The data matrix (nvert * ntime) """ fid = open(filename, 'wb') # write start time in ms - fid.write(np.array(1000*stc['tmin'], dtype='>f4').tostring()) + fid.write(np.array(1000*tmin, dtype='>f4').tostring()) # write sampling rate in ms - fid.write(np.array(1000*stc['tstep'], dtype='>f4').tostring()) + fid.write(np.array(1000*tstep, dtype='>f4').tostring()) # write number of vertices - fid.write(np.array(stc['vertices'].shape[0], dtype='>I4').tostring()) + fid.write(np.array(vertices.shape[0], dtype='>I4').tostring()) # write the vertex indices - fid.write(np.array(stc['vertices'], dtype='>I4').tostring()) + fid.write(np.array(vertices, dtype='>I4').tostring()) # write the number of timepts - fid.write(np.array(stc['data'].shape[1], dtype='>I4').tostring()) + fid.write(np.array(data.shape[1], dtype='>I4').tostring()) # # write the data # - fid.write(np.array(stc['data'].T, dtype='>f4').tostring()) + fid.write(np.array(data.T, dtype='>f4').tostring()) # close the file 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
