[med-svn] [python-mne] 328/376: API : s/source_induced_power/source_band_induced_power ENH : new source_induced_power to computer power and phase lock in a label
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:23:16 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 eb9126443ca867eff46cdaa9f66894ff217bba40
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Fri Jul 22 16:39:00 2011 -0400
API : s/source_induced_power/source_band_induced_power
ENH : new source_induced_power to computer power and phase lock in a label
---
.../plot_source_label_time_frequency.py | 84 +++++++
.../plot_source_space_time_frequency.py | 6 +-
mne/label.py | 20 ++
mne/minimum_norm/__init__.py | 2 +-
mne/minimum_norm/inverse.py | 14 +-
mne/minimum_norm/tests/test_time_frequency.py | 15 +-
mne/minimum_norm/time_frequency.py | 251 +++++++++++++++++----
7 files changed, 323 insertions(+), 69 deletions(-)
diff --git a/examples/time_frequency/plot_source_label_time_frequency.py b/examples/time_frequency/plot_source_label_time_frequency.py
new file mode 100644
index 0000000..5c53aef
--- /dev/null
+++ b/examples/time_frequency/plot_source_label_time_frequency.py
@@ -0,0 +1,84 @@
+"""
+=========================================================
+Compute power and phase lock in label of the source space
+=========================================================
+
+Returns time-frequency maps of induced power and phase lock
+in the source space. The inverse method is linear based on dSPM inverse operator.
+
+"""
+# Authors: Alexandre Gramfort <gramfort at nmr.mgh.harvard.edu>
+#
+# License: BSD (3-clause)
+
+print __doc__
+
+import numpy as np
+
+import mne
+from mne import fiff
+from mne.datasets import sample
+from mne.minimum_norm import read_inverse_operator, source_induced_power
+
+###############################################################################
+# Set parameters
+data_path = sample.data_path('..')
+raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
+fname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'
+label_name = 'Aud-lh'
+fname_label = data_path + '/MEG/sample/labels/%s.label' % label_name
+
+tmin, tmax, event_id = -0.2, 0.5, 1
+
+# Setup for reading the raw data
+raw = fiff.Raw(raw_fname)
+events = mne.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
+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
+frequencies = np.arange(7, 30, 3) # define frequencies of interest
+label = mne.read_label(fname_label)
+power, phase_lock = source_induced_power(epochs, inverse_operator, frequencies,
+ label, baseline=(None, 0), baseline_mode='logratio',
+ n_cycles=2, n_jobs=1)
+
+power = np.mean(power, axis=0) # average over sources
+phase_lock = np.mean(phase_lock, axis=0) # average over sources
+times = epochs.times
+
+###############################################################################
+# View time-frequency plots
+import pylab as pl
+pl.clf()
+pl.subplots_adjust(0.1, 0.08, 0.96, 0.94, 0.2, 0.43)
+pl.subplot(2, 1, 1)
+pl.imshow(20 * power, extent=[times[0], times[-1],
+ frequencies[0], frequencies[-1]],
+ aspect='auto', origin='lower')
+pl.xlabel('Time (s)')
+pl.ylabel('Frequency (Hz)')
+pl.title('Induced power in %s' % label_name)
+pl.colorbar()
+
+pl.subplot(2, 1, 2)
+pl.imshow(phase_lock, extent=[times[0], times[-1],
+ frequencies[0], frequencies[-1]],
+ aspect='auto', origin='lower')
+pl.xlabel('Time (s)')
+pl.ylabel('Frequency (Hz)')
+pl.title('Phase-lock in %s' % label_name)
+pl.colorbar()
+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 b5bbcc4..d34b9e2 100644
--- a/examples/time_frequency/plot_source_space_time_frequency.py
+++ b/examples/time_frequency/plot_source_space_time_frequency.py
@@ -17,7 +17,7 @@ print __doc__
import mne
from mne import fiff
from mne.datasets import sample
-from mne.minimum_norm import read_inverse_operator, source_induced_power
+from mne.minimum_norm import read_inverse_operator, source_band_induced_power
###############################################################################
# Set parameters
@@ -48,8 +48,8 @@ epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
# Compute a source estimate per frequency band
bands = dict(alpha=[9, 11], beta=[18, 22])
-stcs = source_induced_power(epochs, inverse_operator, bands, n_cycles=2,
- use_fft=False, n_jobs=-1)
+stcs = source_band_induced_power(epochs, inverse_operator, bands, n_cycles=2,
+ use_fft=False, n_jobs=1)
for b, stc in stcs.iteritems():
stc.save('induced_power_%s' % b)
diff --git a/mne/label.py b/mne/label.py
index da5fcad..1850fe8 100644
--- a/mne/label.py
+++ b/mne/label.py
@@ -80,3 +80,23 @@ 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/__init__.py b/mne/minimum_norm/__init__.py
index 273ffa8..dbf8766 100644
--- a/mne/minimum_norm/__init__.py
+++ b/mne/minimum_norm/__init__.py
@@ -1,3 +1,3 @@
from .inverse import read_inverse_operator, apply_inverse, minimum_norm, \
apply_inverse_raw
-from .time_frequency import source_induced_power
+from .time_frequency import source_band_induced_power, source_induced_power
diff --git a/mne/minimum_norm/inverse.py b/mne/minimum_norm/inverse.py
index 1496859..41bc7e8 100644
--- a/mne/minimum_norm/inverse.py
+++ b/mne/minimum_norm/inverse.py
@@ -21,7 +21,7 @@ from ..source_space import read_source_spaces_from_tree, find_source_space_hemi
from ..forward import _block_diag
from ..transforms import invert_transform, transform_source_space_to
from ..source_estimate import SourceEstimate
-
+from ..label import source_space_vertices
def read_inverse_operator(fname):
"""Read the inverse operator decomposition from a FIF file
@@ -551,17 +551,7 @@ def apply_inverse_raw(raw, inverse_operator, lambda2, dSPM=True,
noise_norm = inv['noisenorm'][:, None]
if label is not None:
- 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
-
+ 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:
diff --git a/mne/minimum_norm/tests/test_time_frequency.py b/mne/minimum_norm/tests/test_time_frequency.py
index 3622b90..3059af8 100644
--- a/mne/minimum_norm/tests/test_time_frequency.py
+++ b/mne/minimum_norm/tests/test_time_frequency.py
@@ -1,12 +1,13 @@
import os.path as op
import numpy as np
-from numpy.testing import assert_array_almost_equal, assert_equal
+from numpy.testing import assert_array_almost_equal
from ...datasets import sample
from ... import fiff, find_events, Epochs
+from ...label import read_label
from ..inverse import read_inverse_operator
-from ..time_frequency import source_induced_power
+from ..time_frequency import source_band_induced_power
examples_folder = op.join(op.dirname(__file__), '..', '..', '..', 'examples')
@@ -15,6 +16,7 @@ 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')
+fname_label = op.join(data_path, 'MEG', 'sample', 'labels', 'Aud-lh.label')
def test_tfr_with_inverse_operator():
@@ -43,16 +45,17 @@ def test_tfr_with_inverse_operator():
# Compute a source estimate per frequency band
bands = dict(alpha=[10, 10])
+ label = read_label(fname_label)
- stcs = source_induced_power(epochs, inverse_operator, bands, n_cycles=2,
- use_fft=False, pca=True)
+ stcs = source_band_induced_power(epochs, inverse_operator, bands, n_cycles=2,
+ use_fft=False, pca=True, label=label)
stc = stcs['alpha']
assert len(stcs) == len(bands.keys())
assert np.all(stc.data > 0)
assert_array_almost_equal(stc.times, epochs.times)
- stcs_no_pca = source_induced_power(epochs, inverse_operator, bands, n_cycles=2,
- use_fft=False, pca=False)
+ stcs_no_pca = source_band_induced_power(epochs, inverse_operator, bands,
+ n_cycles=2, use_fft=False, pca=False, label=label)
assert_array_almost_equal(stcs['alpha'].data, stcs_no_pca['alpha'].data)
diff --git a/mne/minimum_norm/time_frequency.py b/mne/minimum_norm/time_frequency.py
index b36adce..d46aaed 100644
--- a/mne/minimum_norm/time_frequency.py
+++ b/mne/minimum_norm/time_frequency.py
@@ -11,11 +11,48 @@ from ..time_frequency.tfr import cwt, morlet
from ..baseline import rescale
from .inverse import combine_xyz, prepare_inverse_operator
from ..parallel import parallel_func
+from ..label import source_space_vertices
-def _compute_power(data, K, sel, Ws, source_ori, use_fft, Vh):
+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"""
- power = 0
+ 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
@@ -23,30 +60,49 @@ def _compute_power(data, K, sel, Ws, source_ori, use_fft, Vh):
if Vh is not None:
e = np.dot(Vh, e) # reducing data rank
- for w in Ws:
+ for f, w in enumerate(Ws):
tfr = cwt(e, [w], use_fft=use_fft)
tfr = np.asfortranarray(tfr.reshape(len(e), -1))
- for t in [np.real(tfr), np.imag(tfr)]:
+ # 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)
-
- power += sol
+ pow_f += sol
del sol
- return power
+ 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, bands, lambda2=1.0 / 9.0,
- dSPM=True, n_cycles=5, df=1, use_fft=False,
- baseline=None, baseline_mode='logratio', pca=True,
- subtract_evoked=False, n_jobs=1):
- """Compute source space induced power
+
+def source_band_induced_power(epochs, inverse_operator, bands, label=None,
+ lambda2=1.0 / 9.0, dSPM=True, n_cycles=5, df=1,
+ use_fft=False, baseline=None,
+ baseline_mode='logratio', pca=True,
+ n_jobs=1):
+ """Compute source space induced power in given frequency bands
Parameters
----------
@@ -56,6 +112,8 @@ def source_induced_power(epochs, inverse_operator, bands, lambda2=1.0 / 9.0,
The inverse operator
bands: dict
Example : bands = dict(alpha=[8, 9])
+ label: Label
+ Restricts the source estimates to a given label
lambda2: float
The regularization parameter of the minimum norm
dSPM: bool
@@ -83,23 +141,61 @@ def source_induced_power(epochs, inverse_operator, bands, lambda2=1.0 / 9.0,
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.
n_jobs: int
Number of jobs to run in parallel
"""
- parallel, my_compute_power, n_jobs = parallel_func(_compute_power, n_jobs)
+ frequencies = np.concatenate([np.arange(band[0], band[1] + df / 2.0, df)
+ for _, band in bands.iteritems()])
+
+ powers, _, lh_vertno, rh_vertno = _source_induced_power(epochs,
+ inverse_operator, frequencies,
+ label=label,
+ lambda2=lambda2, dSPM=dSPM,
+ n_cycles=n_cycles,
+ use_fft=use_fft, pca=pca, n_jobs=n_jobs,
+ with_plv=False)
+
+ Fs = epochs.info['sfreq'] # sampling in Hz
+ stcs = dict()
+
+ for name, band in bands.iteritems():
+ idx = [k for k, f in enumerate(frequencies) if band[0] <= f <= band[1]]
+
+ # average power in band + mean over epochs
+ power = np.mean(powers[:, idx, :], axis=1)
+
+ # Run baseline correction
+ 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()
+
+ stcs[name] = stc
+
+ print '[done]'
+
+ return stcs
+
+
+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):
+ """Aux function for source_induced_power
+ """
+ parallel, my_compute_pow_plv, n_jobs = parallel_func(_compute_pow_plv, n_jobs)
#
# Set up the inverse according to the parameters
#
epochs_data = epochs.get_data()
- if subtract_evoked: # subtract with a copy not to touch epochs
- epochs_data = epochs_data - np.mean(epochs_data, axis=0)
-
nave = len(epochs_data) # XXX : can do better when no preload
inv = prepare_inverse_operator(inverse_operator, nave, lambda2, dSPM)
@@ -148,43 +244,104 @@ def source_induced_power(epochs, inverse_operator, bands, lambda2=1.0 / 9.0,
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
- stcs = dict()
- src = inv['src']
+ print 'Computing source power ...'
- for name, band in bands.iteritems():
- print 'Computing power in band %s [%s, %s] Hz...' % (name, band[0],
- band[1])
+ Ws = morlet(Fs, frequencies, n_cycles=n_cycles)
- freqs = np.arange(band[0], band[1] + df / 2.0, df) # frequencies
- Ws = morlet(Fs, freqs, n_cycles=n_cycles)
+ 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)
+ 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
- power = sum(parallel(my_compute_power(data, K, sel, Ws,
- inv['source_ori'], use_fft, Vh)
- for data in np.array_split(epochs_data, n_jobs)))
+ if with_plv:
+ 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:
- # print '(dSPM)...',
- power *= inv['noisenorm'][:, None] ** 2
+ if dSPM:
+ power *= noise_norm[:, None, None] ** 2
- # average power in band + mean over epochs
- power /= len(epochs_data) * len(freqs)
+ return power, plv, lh_vertno, rh_vertno
- # Run baseline correction
- 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 = src[0]['vertno']
- stc.rh_vertno = src[1]['vertno']
- stc._init_times()
+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,
+ baseline_mode='logratio', pca=True, n_jobs=1):
+ """Compute induced power and phase lock
- stcs[name] = stc
+ Computation can optionaly be restricted in a label.
- print '[done]'
+ Parameters
+ ----------
+ epochs: instance of Epochs
+ The epochs
+ inverse_operator: instance of inverse operator
+ The inverse operator
+ label: Label
+ Restricts the source estimates to a given label
+ frequencies: array
+ Array of frequencies of interest
+ lambda2: float
+ The regularization parameter of the minimum norm
+ dSPM: bool
+ Do dSPM or not?
+ n_cycles: int
+ Number of cycles
+ use_fft: bool
+ Do convolutions in time or frequency domain with FFT
+ 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)
+ the interval is between "a (s)" and "b (s)".
+ If a is None the beginning of the data is used
+ and if b is None then b is set to the end of the interval.
+ If baseline is equal ot (None, None) all the time
+ interval is used.
+ baseline_mode: None | 'logratio' | 'zscore'
+ Do baseline correction with ratio (power is divided by mean
+ 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)
+ n_jobs: int
+ Number of jobs to run in parallel
+ """
- return stcs
+ 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)
+
+ # Run baseline correction
+ 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