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 097c81207af3f0a31d4cadac63a3ebd38d49bc4b Author: Alexandre Gramfort <[email protected]> Date: Sun Apr 24 16:34:25 2011 -0400 ENH : speed improvement e.g. in source_induced_power --- examples/plot_morph_data.py | 3 +- .../plot_source_space_time_frequency.py | 10 ++-- mne/minimum_norm/inverse.py | 11 ++++- mne/minimum_norm/tests/test_time_frequency.py | 54 ++++++++++++++++++++++ mne/minimum_norm/time_frequency.py | 32 ++++++------- mne/source_estimate.py | 6 +-- 6 files changed, 89 insertions(+), 27 deletions(-) diff --git a/examples/plot_morph_data.py b/examples/plot_morph_data.py index 6a8a3e3..9b097a3 100644 --- a/examples/plot_morph_data.py +++ b/examples/plot_morph_data.py @@ -28,7 +28,7 @@ src_fname = data_path + '/MEG/sample/sample_audvis-meg-oct-6-fwd.fif' stc_from = mne.SourceEstimate(fname) src_from = mne.read_source_spaces(src_fname) -stc_to = mne.morph_data(subject_from, subject_to, src_from, stc_from, 3) +stc_to = mne.morph_data(subject_from, subject_to, src_from, stc_from) stc_to.save('%s_audvis-meg' % subject_to) @@ -40,4 +40,3 @@ pl.xlabel('time (ms)') pl.ylabel('Mean Source amplitude') pl.legend() 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 b3a3418..eb48a6b 100644 --- a/examples/time_frequency/plot_source_space_time_frequency.py +++ b/examples/time_frequency/plot_source_space_time_frequency.py @@ -3,7 +3,9 @@ Compute induced power in the source space with dSPM =================================================== -XXX +Returns STC files ie source estimates of induced power +for different bands in the source space. The inverse method +is linear based on dSPM inverse operator. """ # Authors: Alexandre Gramfort <[email protected]> @@ -39,14 +41,14 @@ picks = fiff.pick_types(raw.info, meg=True, eeg=False, eog=True, # Load condition 1 event_id = 1 -# events = events[:5] +events = events[:10] # take 10 events to keep the computation time low 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 -# bands = dict(alpha=[9, 12], beta=[13, 25]) -bands = dict(alpha=[8, 9]) +bands = dict(alpha=[9, 11], beta=[18, 22]) + stcs = source_induced_power(epochs, inverse_operator, bands, n_cycles=2, use_fft=False) diff --git a/mne/minimum_norm/inverse.py b/mne/minimum_norm/inverse.py index 5e10929..48060b9 100755 --- a/mne/minimum_norm/inverse.py +++ b/mne/minimum_norm/inverse.py @@ -261,7 +261,7 @@ def read_inverse_operator(fname): # Compute inverse solution -def combine_xyz(vec): +def combine_xyz(vec, square=False): """Compute the three Cartesian components of a vector or matrix together Parameters @@ -281,7 +281,14 @@ def combine_xyz(vec): raise ValueError('Input must have 3N rows') n, p = vec.shape - return np.sqrt((np.abs(vec).reshape(n / 3, 3, p) ** 2).sum(axis=1)) + if np.iscomplexobj(vec): + vec = np.abs(vec) + comb = vec[0::3] ** 2 + comb += vec[1::3] ** 2 + comb += vec[2::3] ** 2 + if not square: + comb = np.sqrt(comb) + return comb def prepare_inverse_operator(orig, nave, lambda2, dSPM): diff --git a/mne/minimum_norm/tests/test_time_frequency.py b/mne/minimum_norm/tests/test_time_frequency.py new file mode 100644 index 0000000..e153129 --- /dev/null +++ b/mne/minimum_norm/tests/test_time_frequency.py @@ -0,0 +1,54 @@ +import os.path as op + +import numpy as np +from numpy.testing import assert_array_almost_equal, assert_equal + +from ...datasets import sample +from ... import fiff, find_events, Epochs +from ..inverse import read_inverse_operator +from ..time_frequency import source_induced_power + + +examples_folder = op.join(op.dirname(__file__), '..', '..', '..', 'examples') +data_path = sample.data_path(examples_folder) +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') + + +def test_tfr_with_inverse_operator(): + """Test MNE inverse computation""" + + tmin, tmax, event_id = -0.2, 0.5, 1 + + # Setup for reading the raw data + raw = fiff.Raw(fname_data) + events = 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 + events = events[:3] # take 3 events to keep the computation time low + epochs = 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 + bands = dict(alpha=[10, 10]) + + stcs = source_induced_power(epochs, inverse_operator, bands, n_cycles=2, + use_fft=False) + + stc = stcs['alpha'] + assert len(stcs) == len(bands.keys()) + assert np.all(stc.data > 0) + assert_array_almost_equal(stc.times, epochs.times) + diff --git a/mne/minimum_norm/time_frequency.py b/mne/minimum_norm/time_frequency.py index 1392e5e..1dbf124 100644 --- a/mne/minimum_norm/time_frequency.py +++ b/mne/minimum_norm/time_frequency.py @@ -94,39 +94,39 @@ def source_induced_power(epochs, inverse_operator, bands, lambda2=1.0 / 9.0, stcs = dict() src = inv['src'] - n_times = len(epochs.times) for name, band in bands.iteritems(): print 'Computing power in band %s [%s, %s] Hz...' % (name, band[0], band[1]) freqs = np.arange(band[0], band[1] + df / 2.0, df) # frequencies - n_freqs = len(freqs) Ws = morlet(Fs, freqs, n_cycles=n_cycles) power = 0 for e in epochs_data: e = e[sel] # keep only selected channels - tfr = cwt(e, Ws, use_fft=use_fft) - tfr = tfr.reshape(len(e), -1) - sol = np.dot(K, tfr) + for w in Ws: + tfr = cwt(e, [w], use_fft=use_fft) + tfr = np.asfortranarray(tfr.reshape(len(e), -1)) - if inv['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI: - # print 'combining the current components...', - sol = combine_xyz(sol) + for t in [np.real(tfr), np.imag(tfr)]: + sol = np.dot(K, t) - if dSPM: - # print '(dSPM)...', - sol *= inv['noisenorm'][:, None] + if inv['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI: + # print 'combining the current components...', + sol = combine_xyz(sol, square=True) - # average power in band - sol = np.mean(np.reshape(sol ** 2, (-1, n_freqs, n_times)), axis=1) - power += sol - del sol + if dSPM: + # print '(dSPM)...', + sol *= inv['noisenorm'][:, None] ** 2 - power /= len(epochs_data) + power += sol + del sol + + # average power in band + mean over epochs + power /= len(epochs_data) * len(freqs) # Run baseline correction if baseline is not None: diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 453b318..671dea6 100755 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -263,7 +263,6 @@ def mesh_edges(tris): edges = edges + edges.T return edges - def morph_data(subject_from, subject_to, src_from, stc_from, grade=5, subjects_dir=None): """Morph a source estimate from one subject to another @@ -313,8 +312,9 @@ def morph_data(subject_from, subject_to, src_from, stc_from, grade=5, idx_use = np.where(src_from[hemi]['inuse'])[0] # XXX n_iter = 100 # max nb of smoothing iterations for k in range(n_iter): - data1 = e[:, idx_use] * np.ones(len(idx_use)) - data[hemi] = e[:, idx_use] * data[hemi] + e_use = e[:, idx_use] + data1 = e_use * np.ones(len(idx_use)) + data[hemi] = e_use * data[hemi] idx_use = np.where(data1)[0] if (k == (n_iter - 1)) or (len(idx_use) >= n_vertices): data[hemi][idx_use, :] /= data1[idx_use][:, None] -- 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
