From eb9126443ca867eff46cdaa9f66894ff217bba40 Mon Sep 17 00:00:00 2001 From: Alexandre Gramfort Date: Fri, 22 Jul 2011 16:39:00 -0400 Subject: [PATCH 1/5] 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(-) create mode 100644 examples/time_frequency/plot_source_label_time_frequency.py 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 00000000000..5c53aef314f --- /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 +# +# 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 b5bbcc474ef..d34b9e2471c 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 @@ 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 @@ # 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 da5fcadb3b6..1850fe892c4 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 273ffa899a6..dbf87664434 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 1496859e150..41bc7e8c6c3 100644 --- a/mne/minimum_norm/inverse.py +++ b/mne/minimum_norm/inverse.py @@ -21,7 +21,7 @@ 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 3622b901260..3059af8c9b9 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 @@ '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 b36adce5969..d46aaed2642 100644 --- a/mne/minimum_norm/time_frequency.py +++ b/mne/minimum_norm/time_frequency.py @@ -11,11 +11,48 @@ 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 From 11d22fd1affc3071ccb4dc9e5bbef9f87c2d2234 Mon Sep 17 00:00:00 2001 From: Alexandre Gramfort Date: Fri, 16 Sep 2011 14:48:39 -0400 Subject: [PATCH 2/5] 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 1850fe892c4..da5fcadb3b6 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 550e4942fff..6423e6769fa 100644 --- a/mne/minimum_norm/tests/test_time_frequency.py +++ b/mne/minimum_norm/tests/test_time_frequency.py @@ -21,7 +21,7 @@ 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 d46aaed2642..b596046c26c 100644 --- a/mne/minimum_norm/time_frequency.py +++ b/mne/minimum_norm/time_frequency.py @@ -6,95 +6,11 @@ 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 From e7d08744416003adf2564c2b36af0677cd7b102d Mon Sep 17 00:00:00 2001 From: Alexandre Gramfort Date: Fri, 16 Sep 2011 16:47:22 -0400 Subject: [PATCH 3/5] DOC: using scipy-lectures as tuto --- doc/source/getting_started.rst | 24 ++---------------------- 1 file changed, 2 insertions(+), 22 deletions(-) diff --git a/doc/source/getting_started.rst b/doc/source/getting_started.rst index 7ed6dc4b18d..dde14f99287 100755 --- a/doc/source/getting_started.rst +++ b/doc/source/getting_started.rst @@ -52,26 +52,6 @@ If you get a new prompt with no error messages, you should be good to go. Learning Python --------------- -If you are new to Python here are some online courses: +If you are new to Python here is a very good place to get started: - * http://docs.python.org/ - - * http://learnpythonthehardway.org/book/ - - * http://code.google.com/intl/fr/edu/languages/google-python-class/ - - * http://homepage.mac.com/s_lott/books/python.html - -More specific tutorials on scientific computing with Python using Numpy and Scipy see: - - * http://showmedo.com/videotutorials/numpy (video) - - * http://www.scipy.org/NumPy_for_Matlab_Users - - * http://mathesaurus.sourceforge.net/matlab-numpy.html - -Some discussions online are: - - * http://stackoverflow.com/questions/111857/what-did-you-use-to-teach-yourself-python - - * http://stackoverflow.com/questions/17988/how-to-learn-python + * http://scipy-lectures.github.com From 87c0dffa0521a91ed400829f017779bc79657280 Mon Sep 17 00:00:00 2001 From: Alexandre Gramfort Date: Fri, 16 Sep 2011 17:22:57 -0400 Subject: [PATCH 4/5] prettify example + adding test for source_induced_power --- .../plot_source_label_time_frequency.py | 4 ++-- mne/baseline.py | 7 +++++-- mne/minimum_norm/tests/test_time_frequency.py | 20 ++++++++++++++++--- 3 files changed, 24 insertions(+), 7 deletions(-) diff --git a/examples/time_frequency/plot_source_label_time_frequency.py b/examples/time_frequency/plot_source_label_time_frequency.py index 5c53aef314f..17621d49fc9 100644 --- a/examples/time_frequency/plot_source_label_time_frequency.py +++ b/examples/time_frequency/plot_source_label_time_frequency.py @@ -49,10 +49,10 @@ preload=True) # Compute a source estimate per frequency band -frequencies = np.arange(7, 30, 3) # define frequencies of interest +frequencies = np.arange(7, 30, 2) # 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', + label, baseline=(-0.1, 0), baseline_mode='percent', n_cycles=2, n_jobs=1) power = np.mean(power, axis=0) # average over sources diff --git a/mne/baseline.py b/mne/baseline.py index 6cb7c11a4f0..16435d4f32d 100644 --- a/mne/baseline.py +++ b/mne/baseline.py @@ -29,7 +29,7 @@ def rescale(data, times, baseline, mode, verbose=True, copy=True): If baseline is equal ot (None, None) all the time interval is used. - mode: 'logratio' | 'ratio' | 'zscore' | 'mean' + mode: 'logratio' | 'ratio' | 'zscore' | 'mean' | 'percent' 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, @@ -44,7 +44,7 @@ def rescale(data, times, baseline, mode, verbose=True, copy=True): if copy: data = data.copy() - valid_modes = ['logratio', 'ratio', 'zscore', 'mean'] + valid_modes = ['logratio', 'ratio', 'zscore', 'mean', 'percent'] if mode not in valid_modes: raise Exception('mode should be any of : %s' % valid_modes) @@ -73,6 +73,9 @@ def rescale(data, times, baseline, mode, verbose=True, copy=True): std = np.std(data[..., imin:imax], axis=-1)[..., None] data -= mean data /= std + elif mode == 'percent': + data -= mean + data /= mean elif verbose: print "No baseline correction applied..." diff --git a/mne/minimum_norm/tests/test_time_frequency.py b/mne/minimum_norm/tests/test_time_frequency.py index 6423e6769fa..c231db03780 100644 --- a/mne/minimum_norm/tests/test_time_frequency.py +++ b/mne/minimum_norm/tests/test_time_frequency.py @@ -8,7 +8,7 @@ from ... import fiff, find_events, Epochs from ...label import read_label from ..inverse import read_inverse_operator -from ..time_frequency import source_band_induced_power +from ..time_frequency import source_band_induced_power, source_induced_power examples_folder = op.join(op.dirname(__file__), '..', '..', '..', 'examples') @@ -39,8 +39,8 @@ def test_tfr_with_inverse_operator(): # 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, + events3 = events[:3] # take 3 events to keep the computation time low + epochs = Epochs(raw, events3, event_id, tmin, tmax, picks=picks, baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6), preload=True) @@ -60,3 +60,17 @@ def test_tfr_with_inverse_operator(): n_cycles=2, use_fft=False, pca=False, label=label) assert_array_almost_equal(stcs['alpha'].data, stcs_no_pca['alpha'].data) + + # Compute a source estimate per frequency band + events = find_events(raw) + epochs = Epochs(raw, events[:10], event_id, tmin, tmax, picks=picks, + baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6), + preload=True) + + frequencies = np.arange(7, 30, 2) # define frequencies of interest + power, phase_lock = source_induced_power(epochs, inverse_operator, + frequencies, label, baseline=(-0.1, 0), + baseline_mode='percent', n_cycles=2, n_jobs=1) + assert_true(np.all(phase_lock > 0)) + assert_true(np.all(phase_lock < 1)) + assert_true(np.max(power) > 10) From abac1f61c1a9454189f8278b400f6c34b78168ae Mon Sep 17 00:00:00 2001 From: Alexandre Gramfort Date: Fri, 16 Sep 2011 21:31:06 -0400 Subject: [PATCH 5/5] ENH: misc extra line --- mne/minimum_norm/tests/test_time_frequency.py | 1 - 1 file changed, 1 deletion(-) diff --git a/mne/minimum_norm/tests/test_time_frequency.py b/mne/minimum_norm/tests/test_time_frequency.py index c231db03780..c6a1d9c5ce4 100644 --- a/mne/minimum_norm/tests/test_time_frequency.py +++ b/mne/minimum_norm/tests/test_time_frequency.py @@ -62,7 +62,6 @@ def test_tfr_with_inverse_operator(): assert_array_almost_equal(stcs['alpha'].data, stcs_no_pca['alpha'].data) # Compute a source estimate per frequency band - events = find_events(raw) epochs = Epochs(raw, events[:10], event_id, tmin, tmax, picks=picks, baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6), preload=True)