[med-svn] [python-mne] 227/376: ENH : speeding up source_induced_power using true data dimension

[Yaroslav Halchenko]

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 a63f8b10811d6068f658f45c2ee77384824a4b45
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date:   Thu Apr 28 11:57:14 2011 -0400

    ENH : speeding up source_induced_power using true data dimension
---
 .../plot_source_space_time_frequency.py            | 12 +++++++++--
 mne/minimum_norm/time_frequency.py                 | 23 +++++++++++++++++++---
 2 files changed, 30 insertions(+), 5 deletions(-)

diff --git a/examples/time_frequency/plot_source_space_time_frequency.py b/examples/time_frequency/plot_source_space_time_frequency.py
index f26f9aa..b7412d8 100644
--- a/examples/time_frequency/plot_source_space_time_frequency.py
+++ b/examples/time_frequency/plot_source_space_time_frequency.py
@@ -55,5 +55,13 @@ stcs = source_induced_power(epochs, inverse_operator, bands, n_cycles=2,
 for b, stc in stcs.iteritems():
     stc.save('induced_power_%s' % b)
 
-# View sources
-plot_source_estimate(inverse_operator['src'], stcs['alpha'])
+###############################################################################
+# plot mean power
+import pylab as pl
+pl.plot(stcs['alpha'].times, stcs['alpha'].data.mean(axis=0), label='Alpha')
+pl.plot(stcs['beta'].times, stcs['beta'].data.mean(axis=0), label='Beta')
+pl.xlabel('Time (ms)')
+pl.ylabel('Power')
+pl.legend()
+pl.title('Mean source induced power')
+pl.show()
diff --git a/mne/minimum_norm/time_frequency.py b/mne/minimum_norm/time_frequency.py
index d70c7b3..6c015da 100644
--- a/mne/minimum_norm/time_frequency.py
+++ b/mne/minimum_norm/time_frequency.py
@@ -3,6 +3,7 @@
 # License: BSD (3-clause)
 
 import numpy as np
+from scipy import linalg
 
 from ..fiff.constants import FIFF
 from ..stc import SourceEstimate
@@ -10,13 +11,16 @@ from ..time_frequency.tfr import cwt, morlet
 from .inverse import combine_xyz, prepare_inverse_operator
 
 
-def _compute_power(data, K, sel, Ws, source_ori, use_fft):
+def _compute_power(data, K, sel, Ws, source_ori, use_fft, Vh):
     """Aux function for source_induced_power"""
     power = 0
 
     for e in data:
         e = e[sel]  # keep only selected channels
 
+        if Vh is not None:
+            e = np.dot(Vh, e)  # reducing data rank
+
         for w in Ws:
             tfr = cwt(e, [w], use_fft=use_fft)
             tfr = np.asfortranarray(tfr.reshape(len(e), -1))
@@ -38,7 +42,7 @@ def _compute_power(data, K, sel, Ws, source_ori, use_fft):
 
 def source_induced_power(epochs, inverse_operator, bands, lambda2=1.0 / 9.0,
                          dSPM=True, n_cycles=5, df=1, use_fft=False,
-                         baseline=None, baseline_mode='logratio',
+                         baseline=None, baseline_mode='logratio', pca=True,
                          subtract_evoked=False, n_jobs=1):
     """Compute source space induced power
 
@@ -73,6 +77,10 @@ def source_induced_power(epochs, inverse_operator, bands, lambda2=1.0 / 9.0,
         power during baseline) or zscore (power is divided by standard
         deviatio of power during baseline after substracting the mean,
         power = [power - mean(power_baseline)] / std(power_baseline))
+    pca: bool
+        If True, the true dimension of data is estimated before running
+        the time frequency transforms. It reduces the computation times
+        e.g. with a dataset that was maxfiltered (true dim is 64)
     subtract_evoked: bool
         If True, the evoked component (average of all epochs) if subtracted
         from each epochs.
@@ -127,6 +135,15 @@ def source_induced_power(epochs, inverse_operator, bands, lambda2=1.0 / 9.0,
                                            inv['whitener'],
                                            inv['proj']])
 
+    if pca:
+        U, s, Vh = linalg.svd(K)
+        rank = np.sum(s > 1e-8*s[0])
+        K = np.dot(K, s[:rank] * U[:, :rank])
+        Vh = Vh[:rank]
+        print 'Reducing data rank to %d' % rank
+    else:
+        Vh = None
+
     #
     #   Transformation into current distributions by weighting the
     #   eigenleads with the weights computed above
@@ -158,7 +175,7 @@ def source_induced_power(epochs, inverse_operator, bands, lambda2=1.0 / 9.0,
         Ws = morlet(Fs, freqs, n_cycles=n_cycles)
 
         power = sum(parallel(my_compute_power(data, K, sel, Ws,
-                                                inv['source_ori'], use_fft)
+                                            inv['source_ori'], use_fft, Vh)
                             for data in np.array_split(epochs_data, n_jobs)))
 
         if dSPM:

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git