[med-svn] [python-mne] 46/376: ENH: make MNE-dSPM computation work
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:22:04 UTC 2015
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 <alexandre.gramfort at inria.fr>
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 <gramfort at nmr.mgh.harvard.edu>
@@ -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
More information about the debian-med-commit
mailing list