[med-svn] [python-mne] 218/353: ENH : speedup + reduce memory of morphing
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:25:00 UTC 2015
This is an automated email from the git hooks/post-receive script.
yoh pushed a commit to tag 0.4
in repository python-mne.
commit 5c6205b566c1097b115e14bf467c51664060cf71
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Thu Jun 21 17:34:05 2012 +0200
ENH : speedup + reduce memory of morphing
---
doc/source/whats_new.rst | 2 +
examples/inverse/plot_morph_data.py | 2 +-
mne/source_estimate.py | 84 ++++++++++++++++++++++++-------------
mne/tests/test_source_estimate.py | 10 +++--
4 files changed, 66 insertions(+), 32 deletions(-)
diff --git a/doc/source/whats_new.rst b/doc/source/whats_new.rst
index 97d87d7..6dc2680 100644
--- a/doc/source/whats_new.rst
+++ b/doc/source/whats_new.rst
@@ -7,6 +7,8 @@ Current
Changelog
~~~~~~~~~
+ - speedup + reduce memory of mne.morph_data by `Alex Gramfort`_.
+
- Backporting scipy.signa.firwin2 so filtering works with old scipy by `Alex Gramfort`_.
- LCMV Beamformer by `Alex Gramfort`_.
diff --git a/examples/inverse/plot_morph_data.py b/examples/inverse/plot_morph_data.py
index 749018f..0554ab0 100644
--- a/examples/inverse/plot_morph_data.py
+++ b/examples/inverse/plot_morph_data.py
@@ -28,7 +28,7 @@ src_fname = data_path + '/MEG/sample/sample_audvis-meg-oct-6-fwd.fif'
# Read input stc file
stc_from = mne.SourceEstimate(fname)
# Morph
-stc_to = mne.morph_data(subject_from, subject_to, stc_from)
+stc_to = mne.morph_data(subject_from, subject_to, stc_from, n_jobs=1)
stc_to.save('%s_audvis-meg' % subject_to)
diff --git a/mne/source_estimate.py b/mne/source_estimate.py
index 865fb67..cccb049 100644
--- a/mne/source_estimate.py
+++ b/mne/source_estimate.py
@@ -6,9 +6,12 @@
import os
import copy
+from math import ceil
import numpy as np
from scipy import sparse
+from .parallel import parallel_func
+
def read_stc(filename):
"""Read an STC file
@@ -549,8 +552,38 @@ def mesh_edges(tris):
return edges
+def _morph_buffer(data, idx_use, e, smooth, n_vertices, nearest, maps):
+ n_iter = 100 # max nb of smoothing iterations
+ for k in range(n_iter):
+ e_use = e[:, idx_use]
+ data1 = e_use * np.ones(len(idx_use))
+ data = e_use * data
+ idx_use = np.where(data1)[0]
+ if smooth is None:
+ if ((k == (n_iter - 1)) or (len(idx_use) >= n_vertices)):
+ # stop when source space in filled with non-zeros
+ break
+ elif k == (smooth - 1):
+ break
+ data = data[idx_use, :] / data1[idx_use][:, None]
+
+ data[idx_use, :] /= data1[idx_use][:, None]
+ print ' %d smooth iterations done.' % (k + 1)
+ data_morphed = maps[nearest, :] * data
+ return data_morphed
+
+
+def _compute_nearest(xhs, rr):
+ nearest = np.zeros(len(rr), dtype=np.int)
+ dr = 32
+ for k in range(0, len(rr), dr):
+ dots = np.dot(rr[k:k + dr], xhs.T)
+ nearest[k:k + dr] = np.argmax(dots, axis=1)
+ return nearest
+
+
def morph_data(subject_from, subject_to, stc_from, grade=5, smooth=None,
- subjects_dir=None):
+ subjects_dir=None, buffer_size=64, n_jobs=1, verbose=0):
"""Morph a source estimate from one subject to another
The functions requires to set MNE_ROOT and SUBJECTS_DIR variables.
@@ -571,6 +604,13 @@ def morph_data(subject_from, subject_to, stc_from, grade=5, smooth=None,
with non-zero values.
subjects_dir : string
Path to SUBJECTS_DIR is not set in the environment
+ buffer_size : int
+ Morph data in chunks of `buffer_size` time instants.
+ Saves memory when morphing long time intervals.
+ n_jobs: int
+ Number of jobs to run in parallel
+ verbose: int
+ Verbosity level.
Returns
-------
@@ -610,18 +650,14 @@ def morph_data(subject_from, subject_to, stc_from, grade=5, smooth=None,
ico = s
break
- nearest = np.zeros((2, ico['np']), dtype=np.int)
-
lhs /= np.sqrt(np.sum(lhs ** 2, axis=1))[:, None]
rhs /= np.sqrt(np.sum(rhs ** 2, axis=1))[:, None]
- rr = ico['rr']
- dr = 16
- for k in range(0, len(rr), dr):
- dots = np.dot(rr[k:k + dr], lhs.T)
- nearest[0, k:k + dr] = np.argmax(dots, axis=1)
- dots = np.dot(rr[k:k + dr], rhs.T)
- nearest[1, k:k + dr] = np.argmax(dots, axis=1)
+ # Compute nearest vertices in high dim mesh
+ parallel, my_compute_nearest, _ = \
+ parallel_func(_compute_nearest, n_jobs, verbose)
+ lhs, rhs, rr = [a.astype(np.float32) for a in [lhs, rhs, ico['rr']]]
+ nearest = parallel(my_compute_nearest(xhs, rr) for xhs in [lhs, rhs])
# morph the data
maps = read_morph_map(subject_from, subject_to, subjects_dir)
@@ -631,6 +667,11 @@ def morph_data(subject_from, subject_to, stc_from, grade=5, smooth=None,
data = [lh_data, rh_data]
data_morphed = [None, None]
+ n_chunks = ceil(stc_from.data.shape[1] / float(buffer_size))
+
+ parallel, my_morph_buffer, _ = \
+ parallel_func(_morph_buffer, n_jobs, verbose)
+
for hemi in [0, 1]:
e = mesh_edges(tris[hemi])
e.data[e.data == 2] = 1
@@ -639,24 +680,11 @@ def morph_data(subject_from, subject_to, stc_from, grade=5, smooth=None,
idx_use = stc_from.vertno[hemi]
if len(idx_use) == 0:
continue
- n_iter = 100 # max nb of smoothing iterations
- for k in range(n_iter):
- 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 smooth is None:
- if ((k == (n_iter - 1)) or (len(idx_use) >= n_vertices)):
- # stop when source space in filled with non-zeros
- break
- elif k == (smooth - 1):
- break
- data[hemi] = data[hemi][idx_use, :] / data1[idx_use][:, None]
-
- data[hemi][idx_use, :] /= data1[idx_use][:, None]
-
- print ' %d smooth iterations done.' % (k + 1)
- data_morphed[hemi] = maps[hemi][nearest[hemi], :] * data[hemi]
+ data_morphed[hemi] = np.concatenate(
+ parallel(my_morph_buffer(data_buffer, idx_use, e, smooth,
+ n_vertices, nearest[hemi], maps[hemi])
+ for data_buffer
+ in np.array_split(data[hemi], n_chunks, axis=1)), axis=1)
stc_to = copy.deepcopy(stc_from)
stc_to.vertno = [nearest[0], nearest[1]]
diff --git a/mne/tests/test_source_estimate.py b/mne/tests/test_source_estimate.py
index 0d3e4bf..5fbaf12 100644
--- a/mne/tests/test_source_estimate.py
+++ b/mne/tests/test_source_estimate.py
@@ -82,14 +82,18 @@ def test_morph_data():
subject_to = 'fsaverage'
fname = op.join(data_path, 'MEG', 'sample', 'sample_audvis-meg')
stc_from = SourceEstimate(fname)
+ stc_from.crop(0.09, 0.1) # for faster computation
stc_to = morph_data(subject_from, subject_to, stc_from,
- grade=3, smooth=12)
-
+ grade=3, smooth=12, buffer_size=1000)
stc_to.save('%s_audvis-meg' % subject_to)
+ stc_to2 = morph_data(subject_from, subject_to, stc_from,
+ grade=3, smooth=12, buffer_size=3)
+ assert_array_almost_equal(stc_to.data, stc_to2.data)
+
mean_from = stc_from.data.mean(axis=0)
mean_to = stc_to.data.mean(axis=0)
- assert_true(np.corrcoef(mean_to, mean_from).min() > 0.99)
+ assert_true(np.corrcoef(mean_to, mean_from).min() > 0.999)
def test_spatio_temporal_tris_connectivity():
--
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