[med-svn] [python-mne] 215/376: ENH : speed improvement e.g. in source_induced_power
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:22:40 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 097c81207af3f0a31d4cadac63a3ebd38d49bc4b
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Sun Apr 24 16:34:25 2011 -0400
ENH : speed improvement e.g. in source_induced_power
---
examples/plot_morph_data.py | 3 +-
.../plot_source_space_time_frequency.py | 10 ++--
mne/minimum_norm/inverse.py | 11 ++++-
mne/minimum_norm/tests/test_time_frequency.py | 54 ++++++++++++++++++++++
mne/minimum_norm/time_frequency.py | 32 ++++++-------
mne/source_estimate.py | 6 +--
6 files changed, 89 insertions(+), 27 deletions(-)
diff --git a/examples/plot_morph_data.py b/examples/plot_morph_data.py
index 6a8a3e3..9b097a3 100644
--- a/examples/plot_morph_data.py
+++ b/examples/plot_morph_data.py
@@ -28,7 +28,7 @@ src_fname = data_path + '/MEG/sample/sample_audvis-meg-oct-6-fwd.fif'
stc_from = mne.SourceEstimate(fname)
src_from = mne.read_source_spaces(src_fname)
-stc_to = mne.morph_data(subject_from, subject_to, src_from, stc_from, 3)
+stc_to = mne.morph_data(subject_from, subject_to, src_from, stc_from)
stc_to.save('%s_audvis-meg' % subject_to)
@@ -40,4 +40,3 @@ pl.xlabel('time (ms)')
pl.ylabel('Mean Source amplitude')
pl.legend()
pl.show()
-
diff --git a/examples/time_frequency/plot_source_space_time_frequency.py b/examples/time_frequency/plot_source_space_time_frequency.py
index b3a3418..eb48a6b 100644
--- a/examples/time_frequency/plot_source_space_time_frequency.py
+++ b/examples/time_frequency/plot_source_space_time_frequency.py
@@ -3,7 +3,9 @@
Compute induced power in the source space with dSPM
===================================================
-XXX
+Returns STC files ie source estimates of induced power
+for different bands in the source space. The inverse method
+is linear based on dSPM inverse operator.
"""
# Authors: Alexandre Gramfort <gramfort at nmr.mgh.harvard.edu>
@@ -39,14 +41,14 @@ picks = fiff.pick_types(raw.info, meg=True, eeg=False, eog=True,
# Load condition 1
event_id = 1
-# events = events[:5]
+events = events[:10] # take 10 events to keep the computation time low
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6),
preload=True)
# Compute a source estimate per frequency band
-# bands = dict(alpha=[9, 12], beta=[13, 25])
-bands = dict(alpha=[8, 9])
+bands = dict(alpha=[9, 11], beta=[18, 22])
+
stcs = source_induced_power(epochs, inverse_operator, bands, n_cycles=2,
use_fft=False)
diff --git a/mne/minimum_norm/inverse.py b/mne/minimum_norm/inverse.py
index 5e10929..48060b9 100755
--- a/mne/minimum_norm/inverse.py
+++ b/mne/minimum_norm/inverse.py
@@ -261,7 +261,7 @@ def read_inverse_operator(fname):
# Compute inverse solution
-def combine_xyz(vec):
+def combine_xyz(vec, square=False):
"""Compute the three Cartesian components of a vector or matrix together
Parameters
@@ -281,7 +281,14 @@ def combine_xyz(vec):
raise ValueError('Input must have 3N rows')
n, p = vec.shape
- return np.sqrt((np.abs(vec).reshape(n / 3, 3, p) ** 2).sum(axis=1))
+ if np.iscomplexobj(vec):
+ vec = np.abs(vec)
+ comb = vec[0::3] ** 2
+ comb += vec[1::3] ** 2
+ comb += vec[2::3] ** 2
+ if not square:
+ comb = np.sqrt(comb)
+ return comb
def prepare_inverse_operator(orig, nave, lambda2, dSPM):
diff --git a/mne/minimum_norm/tests/test_time_frequency.py b/mne/minimum_norm/tests/test_time_frequency.py
new file mode 100644
index 0000000..e153129
--- /dev/null
+++ b/mne/minimum_norm/tests/test_time_frequency.py
@@ -0,0 +1,54 @@
+import os.path as op
+
+import numpy as np
+from numpy.testing import assert_array_almost_equal, assert_equal
+
+from ...datasets import sample
+from ... import fiff, find_events, Epochs
+from ..inverse import read_inverse_operator
+from ..time_frequency import source_induced_power
+
+
+examples_folder = op.join(op.dirname(__file__), '..', '..', '..', 'examples')
+data_path = sample.data_path(examples_folder)
+fname_inv = op.join(data_path, 'MEG', 'sample',
+ 'sample_audvis-meg-oct-6-meg-inv.fif')
+fname_data = op.join(data_path, 'MEG', 'sample',
+ 'sample_audvis_raw.fif')
+
+
+def test_tfr_with_inverse_operator():
+ """Test MNE inverse computation"""
+
+ tmin, tmax, event_id = -0.2, 0.5, 1
+
+ # Setup for reading the raw data
+ raw = fiff.Raw(fname_data)
+ events = find_events(raw)
+ inverse_operator = read_inverse_operator(fname_inv)
+
+ include = []
+ exclude = raw.info['bads'] + ['MEG 2443', 'EEG 053'] # bads + 2 more
+
+ # picks MEG gradiometers
+ picks = fiff.pick_types(raw.info, meg=True, eeg=False, eog=True,
+ stim=False, include=include, exclude=exclude)
+
+ # Load condition 1
+ event_id = 1
+ events = events[:3] # take 3 events to keep the computation time low
+ epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks,
+ baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6),
+ preload=True)
+
+ # Compute a source estimate per frequency band
+ bands = dict(alpha=[10, 10])
+
+ stcs = source_induced_power(epochs, inverse_operator, bands, n_cycles=2,
+ use_fft=False)
+
+ stc = stcs['alpha']
+ assert len(stcs) == len(bands.keys())
+ assert np.all(stc.data > 0)
+ assert_array_almost_equal(stc.times, epochs.times)
+
diff --git a/mne/minimum_norm/time_frequency.py b/mne/minimum_norm/time_frequency.py
index 1392e5e..1dbf124 100644
--- a/mne/minimum_norm/time_frequency.py
+++ b/mne/minimum_norm/time_frequency.py
@@ -94,39 +94,39 @@ def source_induced_power(epochs, inverse_operator, bands, lambda2=1.0 / 9.0,
stcs = dict()
src = inv['src']
- n_times = len(epochs.times)
for name, band in bands.iteritems():
print 'Computing power in band %s [%s, %s] Hz...' % (name, band[0],
band[1])
freqs = np.arange(band[0], band[1] + df / 2.0, df) # frequencies
- n_freqs = len(freqs)
Ws = morlet(Fs, freqs, n_cycles=n_cycles)
power = 0
for e in epochs_data:
e = e[sel] # keep only selected channels
- tfr = cwt(e, Ws, use_fft=use_fft)
- tfr = tfr.reshape(len(e), -1)
- sol = np.dot(K, tfr)
+ for w in Ws:
+ tfr = cwt(e, [w], use_fft=use_fft)
+ tfr = np.asfortranarray(tfr.reshape(len(e), -1))
- if inv['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI:
- # print 'combining the current components...',
- sol = combine_xyz(sol)
+ for t in [np.real(tfr), np.imag(tfr)]:
+ sol = np.dot(K, t)
- if dSPM:
- # print '(dSPM)...',
- sol *= inv['noisenorm'][:, None]
+ if inv['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI:
+ # print 'combining the current components...',
+ sol = combine_xyz(sol, square=True)
- # average power in band
- sol = np.mean(np.reshape(sol ** 2, (-1, n_freqs, n_times)), axis=1)
- power += sol
- del sol
+ if dSPM:
+ # print '(dSPM)...',
+ sol *= inv['noisenorm'][:, None] ** 2
- power /= len(epochs_data)
+ power += sol
+ del sol
+
+ # average power in band + mean over epochs
+ power /= len(epochs_data) * len(freqs)
# Run baseline correction
if baseline is not None:
diff --git a/mne/source_estimate.py b/mne/source_estimate.py
index 453b318..671dea6 100755
--- a/mne/source_estimate.py
+++ b/mne/source_estimate.py
@@ -263,7 +263,6 @@ def mesh_edges(tris):
edges = edges + edges.T
return edges
-
def morph_data(subject_from, subject_to, src_from, stc_from, grade=5,
subjects_dir=None):
"""Morph a source estimate from one subject to another
@@ -313,8 +312,9 @@ def morph_data(subject_from, subject_to, src_from, stc_from, grade=5,
idx_use = np.where(src_from[hemi]['inuse'])[0] # XXX
n_iter = 100 # max nb of smoothing iterations
for k in range(n_iter):
- data1 = e[:, idx_use] * np.ones(len(idx_use))
- data[hemi] = e[:, idx_use] * data[hemi]
+ e_use = e[:, idx_use]
+ data1 = e_use * np.ones(len(idx_use))
+ data[hemi] = e_use * data[hemi]
idx_use = np.where(data1)[0]
if (k == (n_iter - 1)) or (len(idx_use) >= n_vertices):
data[hemi][idx_use, :] /= data1[idx_use][:, None]
--
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