[med-svn] [python-mne] 31/376: ENH : now computing MNE solution
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:22:01 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 b475cbc0ea0bf927df8c48dbc0cd9585f48b7377
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
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
More information about the debian-med-commit
mailing list