[med-svn] [python-mne] 356/376: ENH : using pick_normal for PLV in source space + cleanup
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:23:21 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 11d22fd1affc3071ccb4dc9e5bbef9f87c2d2234
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Fri Sep 16 14:48:39 2011 -0400
ENH : using pick_normal for PLV in source space + cleanup
---
mne/label.py | 20 ---
mne/minimum_norm/tests/test_time_frequency.py | 2 +-
mne/minimum_norm/time_frequency.py | 230 ++++++++++----------------
3 files changed, 87 insertions(+), 165 deletions(-)
diff --git a/mne/label.py b/mne/label.py
index 1850fe8..da5fcad 100644
--- a/mne/label.py
+++ b/mne/label.py
@@ -80,23 +80,3 @@ def label_time_courses(labelfile, stcfile):
times = stc['tmin'] + stc['tstep'] * np.arange(stc['data'].shape[1])
return values, times, vertices
-
-
-def source_space_vertices(label, src):
- """XXX
- """
- lh_vertno = src[0]['vertno']
- rh_vertno = src[1]['vertno']
-
- if label['hemi'] == 'lh':
- vertno_sel = np.intersect1d(lh_vertno, label['vertices'])
- src_sel = np.searchsorted(lh_vertno, vertno_sel)
- lh_vertno = vertno_sel
- rh_vertno = np.array([])
- elif label['hemi'] == 'rh':
- vertno_sel = np.intersect1d(rh_vertno, label['vertices'])
- src_sel = np.searchsorted(rh_vertno, vertno_sel) + len(lh_vertno)
- lh_vertno = np.array([])
- rh_vertno = vertno_sel
-
- return src_sel, lh_vertno, rh_vertno
\ No newline at end of file
diff --git a/mne/minimum_norm/tests/test_time_frequency.py b/mne/minimum_norm/tests/test_time_frequency.py
index 550e494..6423e67 100644
--- a/mne/minimum_norm/tests/test_time_frequency.py
+++ b/mne/minimum_norm/tests/test_time_frequency.py
@@ -21,7 +21,7 @@ fname_label = op.join(data_path, 'MEG', 'sample', 'labels', 'Aud-lh.label')
def test_tfr_with_inverse_operator():
- """Test MNE inverse computation"""
+ """Test time freq with MNE inverse computation"""
tmin, tmax, event_id = -0.2, 0.5, 1
diff --git a/mne/minimum_norm/time_frequency.py b/mne/minimum_norm/time_frequency.py
index d46aaed..b596046 100644
--- a/mne/minimum_norm/time_frequency.py
+++ b/mne/minimum_norm/time_frequency.py
@@ -6,95 +6,11 @@ import numpy as np
from scipy import linalg
from ..fiff.constants import FIFF
-from ..source_estimate import SourceEstimate
from ..time_frequency.tfr import cwt, morlet
from ..baseline import rescale
-from .inverse import combine_xyz, prepare_inverse_operator
+from .inverse import combine_xyz, prepare_inverse_operator, _assemble_kernel, \
+ _make_stc
from ..parallel import parallel_func
-from ..label import source_space_vertices
-
-
-def _mean_xyz(vec):
- """Compute mean of the three Cartesian components of a matrix
-
- Parameters
- ----------
- vec : 2d array of shape [3 n x p]
- Input [ x1 y1 z1 ... x_n y_n z_n ] where x1 ... z_n
- can be vectors
-
- Returns
- -------
- comb : array
- Output vector [(x_1 + y_1 + z_1) / 3, ..., (x_n + y_n + z_n) / 3]
- """
- if (vec.shape[0] % 3) != 0:
- raise Exception('Input must have 3N rows')
-
- comb = vec[0::3]
- comb += vec[1::3]
- comb += vec[2::3]
- comb /= 3
- return comb
-
-
-def _compute_pow_plv(data, K, sel, Ws, source_ori, use_fft, Vh, with_plv):
- """Aux function for source_induced_power"""
- n_times = data.shape[2]
- n_freqs = len(Ws)
- n_sources = K.shape[0]
- if source_ori == FIFF.FIFFV_MNE_FREE_ORI:
- n_sources /= 3
-
- shape = (n_sources, n_freqs, n_times)
- power = np.zeros(shape, dtype=np.float) # power
- if with_plv:
- shape = (K.shape[0], n_freqs, n_times)
- plv = np.zeros(shape, dtype=np.complex) # phase lock
- else:
- plv = None
-
- 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 f, w in enumerate(Ws):
- tfr = cwt(e, [w], use_fft=use_fft)
- tfr = np.asfortranarray(tfr.reshape(len(e), -1))
-
- # phase lock and power at freq f
- if with_plv:
- plv_f = np.zeros((K.shape[0], n_times), dtype=np.complex)
- pow_f = np.zeros((n_sources, n_times), dtype=np.float)
-
- for k, t in enumerate([np.real(tfr), np.imag(tfr)]):
- sol = np.dot(K, t)
-
- if with_plv:
- if k == 0: # real
- plv_f += sol
- else: # imag
- plv_f += 1j * sol
-
- if source_ori == FIFF.FIFFV_MNE_FREE_ORI:
- # print 'combining the current components...',
- sol = combine_xyz(sol, square=True)
- else:
- np.power(sol, 2, sol)
- pow_f += sol
- del sol
-
- power[:, f, :] += pow_f
- del pow_f
-
- if with_plv:
- plv_f /= np.abs(plv_f)
- plv[:, f, :] += plv_f
- del plv_f
-
- return power, plv
def source_band_induced_power(epochs, inverse_operator, bands, label=None,
@@ -144,7 +60,6 @@ def source_band_induced_power(epochs, inverse_operator, bands, label=None,
n_jobs: int
Number of jobs to run in parallel
"""
-
frequencies = np.concatenate([np.arange(band[0], band[1] + df / 2.0, df)
for _, band in bands.iteritems()])
@@ -169,14 +84,7 @@ def source_band_induced_power(epochs, inverse_operator, bands, label=None,
power = rescale(power, epochs.times, baseline, baseline_mode,
verbose=True, copy=False)
- stc = SourceEstimate(None)
- stc.data = power
- stc.tmin = epochs.times[0]
- stc.tstep = 1.0 / Fs
- stc.lh_vertno = lh_vertno
- stc.rh_vertno = rh_vertno
- stc._init_times()
-
+ stc = _make_stc(power, epochs.times[0], 1.0 / Fs, lh_vertno, rh_vertno)
stcs[name] = stc
print '[done]'
@@ -184,13 +92,80 @@ def source_band_induced_power(epochs, inverse_operator, bands, label=None,
return stcs
+def _compute_pow_plv(data, K, sel, Ws, source_ori, use_fft, Vh, with_plv,
+ pick_normal):
+ """Aux function for source_induced_power"""
+ n_times = data.shape[2]
+ n_freqs = len(Ws)
+ n_sources = K.shape[0]
+ is_free_ori = False
+ if source_ori == FIFF.FIFFV_MNE_FREE_ORI:
+ is_free_ori = True
+ n_sources /= 3
+
+ shape = (n_sources, n_freqs, n_times)
+ power = np.zeros(shape, dtype=np.float) # power
+ if with_plv:
+ shape = (n_sources, n_freqs, n_times)
+ plv = np.zeros(shape, dtype=np.complex) # phase lock
+ else:
+ plv = None
+
+ 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 f, w in enumerate(Ws):
+ tfr = cwt(e, [w], use_fft=use_fft)
+ tfr = np.asfortranarray(tfr.reshape(len(e), -1))
+
+ # phase lock and power at freq f
+ if with_plv:
+ plv_f = np.zeros((n_sources, n_times), dtype=np.complex)
+ pow_f = np.zeros((n_sources, n_times), dtype=np.float)
+
+ for k, t in enumerate([np.real(tfr), np.imag(tfr)]):
+ sol = np.dot(K, t)
+
+ sol_pick_normal = sol
+ if is_free_ori:
+ sol_pick_normal = sol[2::3]
+
+ if with_plv:
+ if k == 0: # real
+ plv_f += sol_pick_normal
+ else: # imag
+ plv_f += 1j * sol_pick_normal
+
+ if is_free_ori:
+ # print 'combining the current components...',
+ sol = combine_xyz(sol, square=True)
+ else:
+ np.power(sol, 2, sol)
+ pow_f += sol
+ del sol
+
+ power[:, f, :] += pow_f
+ del pow_f
+
+ if with_plv:
+ plv_f /= np.abs(plv_f)
+ plv[:, f, :] += plv_f
+ del plv_f
+
+ return power, plv
+
+
def _source_induced_power(epochs, inverse_operator, frequencies, label=None,
lambda2=1.0 / 9.0, dSPM=True, n_cycles=5,
- use_fft=False, pca=True, n_jobs=1, with_plv=True):
+ use_fft=False, pca=True, pick_normal=True,
+ n_jobs=1, with_plv=True):
"""Aux function for source_induced_power
"""
- parallel, my_compute_pow_plv, n_jobs = parallel_func(_compute_pow_plv, n_jobs)
-
+ parallel, my_compute_pow_plv, n_jobs = parallel_func(_compute_pow_plv,
+ n_jobs)
#
# Set up the inverse according to the parameters
#
@@ -212,10 +187,7 @@ def _source_induced_power(epochs, inverse_operator, frequencies, label=None,
# This does all the data transformations to compute the weights for the
# eigenleads
#
- K = inv['reginv'][:, None] * reduce(np.dot,
- [inv['eigen_fields']['data'],
- inv['whitener'],
- inv['proj']])
+ K, noise_norm, lh_vertno, rh_vertno = _assemble_kernel(inv, label, dSPM)
if pca:
U, s, Vh = linalg.svd(K)
@@ -226,39 +198,6 @@ def _source_induced_power(epochs, inverse_operator, frequencies, label=None,
else:
Vh = None
- #
- # 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)...',
- K = np.dot(inv['eigen_leads']['data'], K)
- else:
- #
- # R^0.5 has to factored in
- #
- # print '(eigenleads need to be weighted)...',
- K = np.sqrt(inv['source_cov']['data'])[:, None] * \
- np.dot(inv['eigen_leads']['data'], K)
-
- noise_norm = inv['noisenorm']
- src = inverse_operator['src']
- lh_vertno = src[0]['vertno']
- rh_vertno = src[1]['vertno']
- if label is not None:
- src_sel, lh_vertno, rh_vertno = source_space_vertices(label, src)
- noise_norm = noise_norm[src_sel]
-
- if inv['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI:
- src_sel = 3 * src_sel
- src_sel = np.c_[src_sel, src_sel + 1, src_sel + 2]
- src_sel = src_sel.ravel()
-
- K = K[src_sel, :]
-
Fs = epochs.info['sfreq'] # sampling in Hz
print 'Computing source power ...'
@@ -268,7 +207,7 @@ def _source_induced_power(epochs, inverse_operator, frequencies, label=None,
n_jobs = min(n_jobs, len(epochs_data))
out = parallel(my_compute_pow_plv(data, K, sel, Ws,
inv['source_ori'], use_fft, Vh,
- with_plv)
+ with_plv, pick_normal)
for data in np.array_split(epochs_data, n_jobs))
power = sum(o[0] for o in out)
power /= len(epochs_data) # average power over epochs
@@ -277,20 +216,18 @@ def _source_induced_power(epochs, inverse_operator, frequencies, label=None,
plv = sum(o[1] for o in out)
plv = np.abs(plv)
plv /= len(epochs_data) # average power over epochs
- if with_plv and (inv['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI):
- plv = _mean_xyz(plv)
else:
plv = None
if dSPM:
- power *= noise_norm[:, None, None] ** 2
+ power *= noise_norm.ravel()[:, None, None] ** 2
return power, plv, lh_vertno, rh_vertno
def source_induced_power(epochs, inverse_operator, frequencies, label=None,
lambda2=1.0 / 9.0, dSPM=True, n_cycles=5,
- use_fft=False, baseline=None,
+ use_fft=False, pick_normal=False, baseline=None,
baseline_mode='logratio', pca=True, n_jobs=1):
"""Compute induced power and phase lock
@@ -314,6 +251,10 @@ def source_induced_power(epochs, inverse_operator, frequencies, label=None,
Number of cycles
use_fft: bool
Do convolutions in time or frequency domain with FFT
+ pick_normal: bool
+ If True, rather than pooling the orientations by taking the norm,
+ only the radial component is kept. This is only implemented
+ when working with loose orientations.
baseline: None (default) or tuple of length 2
The time interval to apply baseline correction.
If None do not apply it. If baseline is (a, b)
@@ -334,14 +275,15 @@ def source_induced_power(epochs, inverse_operator, frequencies, label=None,
n_jobs: int
Number of jobs to run in parallel
"""
-
power, plv, lh_vertno, rh_vertno = _source_induced_power(epochs,
inverse_operator, frequencies,
label, lambda2, dSPM, n_cycles,
- use_fft, pca=True, n_jobs=1)
+ use_fft, pick_normal=pick_normal, pca=pca,
+ n_jobs=n_jobs)
# Run baseline correction
- power = rescale(power, epochs.times, baseline, baseline_mode,
- verbose=True, copy=False)
+ if baseline is not None:
+ power = rescale(power, epochs.times, baseline, baseline_mode,
+ verbose=True, copy=False)
return power, plv
--
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