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 <[email protected]> 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 <[email protected]> +# +# 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 _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
