This is an automated email from the git hooks/post-receive script. yoh pushed a commit to tag 0.4 in repository python-mne.
commit 99e8c0ba3a1fbec6438e814111493f01f403acd5 Author: Alexandre Gramfort <[email protected]> Date: Sun Jun 17 19:36:20 2012 +0300 ENH + API : add support for sLORETA to apply_inverse, apply_inverse_raw, apply_inverse_epochs --- doc/source/python_tutorial.rst | 12 +-- doc/source/whats_new.rst | 2 + examples/inverse/plot_compute_mne_inverse.py | 12 +-- .../plot_compute_mne_inverse_epochs_in_label.py | 7 +- .../plot_compute_mne_inverse_raw_in_label.py | 14 ++-- .../inverse/plot_compute_mne_inverse_volume.py | 11 ++- examples/inverse/plot_make_inverse_operator.py | 7 +- mne/epochs.py | 1 - mne/minimum_norm/inverse.py | 90 ++++++++++++++-------- mne/minimum_norm/tests/test_inverse.py | 31 ++++---- mne/minimum_norm/time_frequency.py | 35 +++++---- 11 files changed, 129 insertions(+), 93 deletions(-) diff --git a/doc/source/python_tutorial.rst b/doc/source/python_tutorial.rst index 1377d4d..d799646 100644 --- a/doc/source/python_tutorial.rst +++ b/doc/source/python_tutorial.rst @@ -126,6 +126,8 @@ First extract events: >>> events = mne.find_events(raw, stim_channel='STI 014') Reading 6450 ... 48149 = 42.956 ... 320.665 secs... [done] + 319 events found + Events id: [ 1 2 3 4 5 32] >>> print events[:5] [[6994 0 2] [7086 0 3] @@ -275,17 +277,17 @@ Define the inverse parameters: >>> snr = 3.0 >>> lambda2 = 1.0 / snr ** 2 - >>> dSPM = True + >>> method = "dSPM" Compute the inverse solution: - >>> stc = apply_inverse(evoked, inverse_operator, lambda2, dSPM) + >>> stc = apply_inverse(evoked, inverse_operator, lambda2, method) Preparing the inverse operator for use... Scaled noise and source covariance from nave = 1 to nave = 72 Created the regularized inverter Created an SSP operator (subspace dimension = 3) Created the whitener using a full noise covariance matrix (3 small eigenvalues omitted) - Computing noise-normalization factors... [done] + Computing noise-normalization factors (dSPM)... [done] Picked 305 channels from the data Computing inverse... (eigenleads need to be weighted)... combining the current components... (dSPM)... [done] @@ -303,13 +305,13 @@ Compute inverse solution during the first 15s: >>> from mne.minimum_norm import apply_inverse_raw >>> start, stop = raw.time_to_index(0, 15) # read the first 15s of data - >>> stc = apply_inverse_raw(raw, inverse_operator, lambda2, dSPM, label, start, stop) + >>> stc = apply_inverse_raw(raw, inverse_operator, lambda2, method, label, start, stop) Preparing the inverse operator for use... Scaled noise and source covariance from nave = 1 to nave = 1 Created the regularized inverter Created an SSP operator (subspace dimension = 3) Created the whitener using a full noise covariance matrix (3 small eigenvalues omitted) - Computing noise-normalization factors... [done] + Computing noise-normalization factors (dSPM)... [done] Picked 305 channels from the data Computing inverse... Reading 6450 ... 8701 = 42.956 ... 57.947 secs... [done] (eigenleads need to be weighted)... combining the current components... [done] diff --git a/doc/source/whats_new.rst b/doc/source/whats_new.rst index 4372f4b..97a53bc 100644 --- a/doc/source/whats_new.rst +++ b/doc/source/whats_new.rst @@ -9,6 +9,8 @@ Current Changelog ~~~~~~~~~ + - Add support for sLORETA in apply_inverse, apply_inverse_raw, apply_inverse_epochs (API Change) by `Alex Gramfort`_. + - Add method to regularize a noise covariance by `Alex Gramfort`_. - Read and write measurement info in forward and inverse operators for interactive visualization in mne_analyze by `Alex Gramfort`_. diff --git a/examples/inverse/plot_compute_mne_inverse.py b/examples/inverse/plot_compute_mne_inverse.py index c67b1fb..ca62f26 100644 --- a/examples/inverse/plot_compute_mne_inverse.py +++ b/examples/inverse/plot_compute_mne_inverse.py @@ -24,24 +24,24 @@ data_path = sample.data_path('..') fname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif' fname_evoked = data_path + '/MEG/sample/sample_audvis-ave.fif' -setno = 0 snr = 3.0 lambda2 = 1.0 / snr ** 2 -dSPM = True +method = "dSPM" # use dSPM method (could also be MNE or sLORETA) # Load data -evoked = Evoked(fname_evoked, setno=setno, baseline=(None, 0)) +evoked = Evoked(fname_evoked, setno=0, baseline=(None, 0)) inverse_operator = read_inverse_operator(fname_inv) # Compute inverse solution -stc = apply_inverse(evoked, inverse_operator, lambda2, dSPM, pick_normal=False) +stc = apply_inverse(evoked, inverse_operator, lambda2, method, + pick_normal=False) # Save result in stc files -stc.save('mne_dSPM_inverse') +stc.save('mne_%s_inverse' % method) ############################################################################### # View activation time-series pl.plot(1e3 * stc.times, stc.data[::100, :].T) pl.xlabel('time (ms)') -pl.ylabel('dSPM value') +pl.ylabel('%s value' % method) pl.show() diff --git a/examples/inverse/plot_compute_mne_inverse_epochs_in_label.py b/examples/inverse/plot_compute_mne_inverse_epochs_in_label.py index 473c64e..6ebcad8 100644 --- a/examples/inverse/plot_compute_mne_inverse_epochs_in_label.py +++ b/examples/inverse/plot_compute_mne_inverse_epochs_in_label.py @@ -30,9 +30,9 @@ label_name = 'Aud-lh' fname_label = data_path + '/MEG/sample/labels/%s.label' % label_name event_id, tmin, tmax = 1, -0.2, 0.5 -snr = 3.0 +snr = 1.0 # use smaller SNR or raw data lambda2 = 1.0 / snr ** 2 -dSPM = True +method = "dSPM" # use dSPM method (could also be MNE or sLORETA) # Load data inverse_operator = read_inverse_operator(fname_inv) @@ -53,7 +53,7 @@ epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks, eog=150e-6)) # Compute inverse solution and stcs for each epoch -stcs = apply_inverse_epochs(epochs, inverse_operator, lambda2, dSPM, label, +stcs = apply_inverse_epochs(epochs, inverse_operator, lambda2, method, label, pick_normal=True) data = sum(stc.data for stc in stcs) / len(stcs) @@ -66,6 +66,7 @@ label_mean_flip = np.mean(flip[:, np.newaxis] * data, axis=0) ############################################################################### # View activation time-series +pl.figure() h0 = pl.plot(1e3 * stcs[0].times, data.T, 'k') h1 = pl.plot(1e3 * stcs[0].times, label_mean, 'r', linewidth=3) h2 = pl.plot(1e3 * stcs[0].times, label_mean_flip, 'g', linewidth=3) diff --git a/examples/inverse/plot_compute_mne_inverse_raw_in_label.py b/examples/inverse/plot_compute_mne_inverse_raw_in_label.py index ee4bcbf..910b9d3 100644 --- a/examples/inverse/plot_compute_mne_inverse_raw_in_label.py +++ b/examples/inverse/plot_compute_mne_inverse_raw_in_label.py @@ -1,9 +1,9 @@ """ ============================================= -Compute MNE-dSPM inverse solution on raw data +Compute sLORETA inverse solution on raw data ============================================= -Compute dSPM inverse solution on raw dataset restricted +Compute sLORETA inverse solution on raw dataset restricted to a brain label and stores the solution in stc files for visualisation. @@ -28,9 +28,9 @@ fname_raw = data_path + '/MEG/sample/sample_audvis_raw.fif' label_name = 'Aud-lh' fname_label = data_path + '/MEG/sample/labels/%s.label' % label_name -snr = 3.0 +snr = 1.0 # use smaller SNR or raw data lambda2 = 1.0 / snr ** 2 -dSPM = True +method = "sLORETA" # use sLORETA method (could also be MNE or dSPM) # Load data raw = Raw(fname_raw) @@ -40,15 +40,15 @@ label = mne.read_label(fname_label) start, stop = raw.time_to_index(0, 15) # read the first 15s of data # Compute inverse solution -stc = apply_inverse_raw(raw, inverse_operator, lambda2, dSPM, label, +stc = apply_inverse_raw(raw, inverse_operator, lambda2, method, label, start, stop, pick_normal=False) # Save result in stc files -stc.save('mne_dSPM_raw_inverse_%s' % label_name) +stc.save('mne_%s_raw_inverse_%s' % (method, label_name)) ############################################################################### # View activation time-series pl.plot(1e3 * stc.times, stc.data[::100, :].T) pl.xlabel('time (ms)') -pl.ylabel('dSPM value') +pl.ylabel('%s value' % method) pl.show() diff --git a/examples/inverse/plot_compute_mne_inverse_volume.py b/examples/inverse/plot_compute_mne_inverse_volume.py index 884a43b..2ecd3db 100644 --- a/examples/inverse/plot_compute_mne_inverse_volume.py +++ b/examples/inverse/plot_compute_mne_inverse_volume.py @@ -25,23 +25,22 @@ data_path = sample.data_path('..') fname_inv = data_path + '/MEG/sample/sample_audvis-meg-vol-7-meg-inv.fif' fname_evoked = data_path + '/MEG/sample/sample_audvis-ave.fif' -setno = 0 snr = 3.0 lambda2 = 1.0 / snr ** 2 -dSPM = True +method = "dSPM" # use dSPM method (could also be MNE or sLORETA) # Load data -evoked = Evoked(fname_evoked, setno=setno, baseline=(None, 0)) +evoked = Evoked(fname_evoked, setno=0, baseline=(None, 0)) inverse_operator = read_inverse_operator(fname_inv) src = inverse_operator['src'] # Compute inverse solution -stc = apply_inverse(evoked, inverse_operator, lambda2, dSPM) +stc = apply_inverse(evoked, inverse_operator, lambda2, method) stc.crop(0.0, 0.2) # Save result in a 4D nifti file -img = mne.save_stc_as_volume('mne_dSPM_inverse.nii.gz', stc, src, - mri_resolution=False) # set to True for full MRI resolution +img = mne.save_stc_as_volume('mne_%s_inverse.nii.gz' % method, stc, + src, mri_resolution=False) # set to True for full MRI resolution data = img.get_data() # plot result (one slice) diff --git a/examples/inverse/plot_make_inverse_operator.py b/examples/inverse/plot_make_inverse_operator.py index c8b55f4..f10e3a4 100644 --- a/examples/inverse/plot_make_inverse_operator.py +++ b/examples/inverse/plot_make_inverse_operator.py @@ -27,13 +27,11 @@ fname_fwd = data_path + '/MEG/sample/sample_audvis-eeg-oct-6-fwd.fif' fname_cov = data_path + '/MEG/sample/sample_audvis-cov.fif' fname_evoked = data_path + '/MEG/sample/sample_audvis-ave.fif' -setno = 0 snr = 3.0 lambda2 = 1.0 / snr ** 2 -dSPM = True # Load data -evoked = Evoked(fname_evoked, setno=setno, baseline=(None, 0)) +evoked = Evoked(fname_evoked, setno=0, baseline=(None, 0)) forward = mne.read_forward_solution(fname_fwd, surf_ori=True) noise_cov = mne.read_cov(fname_cov) @@ -48,7 +46,8 @@ inverse_operator = make_inverse_operator(evoked.info, forward, noise_cov, write_inverse_operator('sample_audvis-eeg-oct-6-eeg-inv.fif', inverse_operator) # Compute inverse solution -stc = apply_inverse(evoked, inverse_operator, lambda2, dSPM, pick_normal=False) +stc = apply_inverse(evoked, inverse_operator, lambda2, "dSPM", + pick_normal=False) # Save result in stc files stc.save('mne_dSPM_inverse') diff --git a/mne/epochs.py b/mne/epochs.py index 1013d40..5d922b3 100644 --- a/mne/epochs.py +++ b/mne/epochs.py @@ -146,7 +146,6 @@ class Epochs(object): self.info['projs'] = activate_proj(self.info['projs'], copy=False) # Add EEG ref reference proj - print "Adding average EEG reference projection." eeg_sel = pick_types(self.info, meg=False, eeg=True) if len(eeg_sel) > 0: eeg_proj = make_eeg_average_ref_proj(self.info) diff --git a/mne/minimum_norm/inverse.py b/mne/minimum_norm/inverse.py index 1b86ac9..8dc74c2 100644 --- a/mne/minimum_norm/inverse.py +++ b/mne/minimum_norm/inverse.py @@ -414,7 +414,7 @@ def _chech_ch_names(inv, info): 'are not present in the data (%s)' % missing_ch_names) -def prepare_inverse_operator(orig, nave, lambda2, dSPM): +def prepare_inverse_operator(orig, nave, lambda2, method): """Prepare an inverse operator for actually computing the inverse Parameters @@ -425,15 +425,14 @@ def prepare_inverse_operator(orig, nave, lambda2, dSPM): Number of averages (scales the noise covariance) lambda2: float The regularization factor. Recommended to be 1 / SNR**2 - dSPM: bool - If True, compute the noise-normalization factors for dSPM. + method: "MNE" | "dSPM" | "sLORETA" + Use mininum norm, dSPM or sLORETA Returns ------- inv: dict Prepared inverse operator """ - if nave <= 0: raise ValueError('The number of averages should be positive') @@ -500,18 +499,24 @@ def prepare_inverse_operator(orig, nave, lambda2, dSPM): # # Finally, compute the noise-normalization factors # - if dSPM: - print ' Computing noise-normalization factors...', + if method in ["dSPM", 'sLORETA']: + if method == "dSPM": + print ' Computing noise-normalization factors (dSPM)...', + noise_weight = inv['reginv'] + else: + print ' Computing noise-normalization factors (sLORETA)...', + noise_weight = inv['reginv'] * \ + np.sqrt((1. + inv['sing'] ** 2 / lambda2)) noise_norm = np.zeros(inv['eigen_leads']['nrow']) nrm2, = linalg.get_blas_funcs(('nrm2',), (noise_norm,)) if inv['eigen_leads_weighted']: for k in range(inv['eigen_leads']['nrow']): - one = inv['eigen_leads']['data'][k, :] * inv['reginv'] + one = inv['eigen_leads']['data'][k, :] * noise_weight noise_norm[k] = nrm2(one) else: for k in range(inv['eigen_leads']['nrow']): one = sqrt(inv['source_cov']['data'][k]) * \ - inv['eigen_leads']['data'][k, :] * inv['reginv'] + inv['eigen_leads']['data'][k, :] * noise_weight noise_norm[k] = nrm2(one) # @@ -536,7 +541,7 @@ def prepare_inverse_operator(orig, nave, lambda2, dSPM): return inv -def _assemble_kernel(inv, label, dSPM, pick_normal): +def _assemble_kernel(inv, label, method, pick_normal): # # Simple matrix multiplication followed by combination of the # current components @@ -546,7 +551,7 @@ def _assemble_kernel(inv, label, dSPM, pick_normal): # eigen_leads = inv['eigen_leads']['data'] source_cov = inv['source_cov']['data'][:, None] - if dSPM: + if method != "MNE": noise_norm = inv['noisenorm'][:, None] src = inv['src'] @@ -555,7 +560,7 @@ def _assemble_kernel(inv, label, dSPM, pick_normal): if label is not None: vertno, src_sel = _get_label_sel(label, inv) - if dSPM: + if method != "MNE": noise_norm = noise_norm[src_sel] if inv['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI: @@ -601,7 +606,7 @@ def _assemble_kernel(inv, label, dSPM, pick_normal): print '(eigenleads need to be weighted)...', K = np.sqrt(source_cov) * np.dot(eigen_leads, trans) - if not dSPM: + if method == "MNE": noise_norm = None return K, noise_norm, vertno @@ -641,8 +646,28 @@ def _get_label_sel(label, inv): return vertno, src_sel -def apply_inverse(evoked, inverse_operator, lambda2, dSPM=True, - pick_normal=False): +def _check_method(method, dSPM): + if dSPM is not None: + warnings.warn('DEPRECATION: The dSPM parameter has been changed to ' + 'method. Please update your code') + method = dSPM + if method == True: + warnings.warn('DEPRECATION:Inverse method should now be "MNE" or ' + '"dSPM" or "sLORETA".') + method = "dSPM" + if method == False: + warnings.warn('DEPRECATION:Inverse method should now be "MNE" or ' + '"dSPM" or "sLORETA".') + method = "MNE" + + if method not in ["MNE", "dSPM", "sLORETA"]: + raise ValueError('method parameter should be "MNE" or "dSPM" ' + 'or "sLORETA".') + return method + + +def apply_inverse(evoked, inverse_operator, lambda2, method="dSPM", + pick_normal=False, dSPM=None): """Apply inverse operator to evoked data Computes a L2-norm inverse solution @@ -657,8 +682,8 @@ def apply_inverse(evoked, inverse_operator, lambda2, dSPM=True, Inverse operator read with mne.read_inverse_operator lambda2: float The regularization parameter - dSPM: bool - do dSPM ? + method: "MNE" | "dSPM" | "sLORETA" + Use mininum norm, dSPM or sLORETA pick_normal: bool If True, rather than pooling the orientations by taking the norm, only the radial component is kept. This is only implemented @@ -669,6 +694,7 @@ def apply_inverse(evoked, inverse_operator, lambda2, dSPM=True, stc: SourceEstimate The source estimates """ + method = _check_method(method, dSPM) # # Set up the inverse according to the parameters # @@ -676,7 +702,7 @@ def apply_inverse(evoked, inverse_operator, lambda2, dSPM=True, _chech_ch_names(inverse_operator, evoked.info) - inv = prepare_inverse_operator(inverse_operator, nave, lambda2, dSPM) + inv = prepare_inverse_operator(inverse_operator, nave, lambda2, method) # # Pick the correct channels from the data # @@ -684,7 +710,7 @@ def apply_inverse(evoked, inverse_operator, lambda2, dSPM=True, print 'Picked %d channels from the data' % len(sel) print 'Computing inverse...', - K, noise_norm, _ = _assemble_kernel(inv, None, dSPM, pick_normal) + K, noise_norm, _ = _assemble_kernel(inv, None, method, pick_normal) sol = np.dot(K, evoked.data[sel]) # apply imaging kernel is_free_ori = (inverse_operator['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI @@ -707,10 +733,10 @@ def apply_inverse(evoked, inverse_operator, lambda2, dSPM=True, return stc -def apply_inverse_raw(raw, inverse_operator, lambda2, dSPM=True, +def apply_inverse_raw(raw, inverse_operator, lambda2, method="dSPM", label=None, start=None, stop=None, nave=1, time_func=None, pick_normal=False, - buffer_size=None): + buffer_size=None, dSPM=None): """Apply inverse operator to Raw data Computes a L2-norm inverse solution @@ -725,8 +751,8 @@ def apply_inverse_raw(raw, inverse_operator, lambda2, dSPM=True, Inverse operator read with mne.read_inverse_operator lambda2: float The regularization parameter - dSPM: bool - do dSPM ? + method: "MNE" | "dSPM" | "sLORETA" + Use mininum norm, dSPM or sLORETA label: Label Restricts the source estimates to a given label start: int @@ -755,12 +781,14 @@ def apply_inverse_raw(raw, inverse_operator, lambda2, dSPM=True, stc: SourceEstimate The source estimates """ + method = _check_method(method, dSPM) + _chech_ch_names(inverse_operator, raw.info) # # Set up the inverse according to the parameters # - inv = prepare_inverse_operator(inverse_operator, nave, lambda2, dSPM) + inv = prepare_inverse_operator(inverse_operator, nave, lambda2, method) # # Pick the correct channels from the data # @@ -773,7 +801,7 @@ def apply_inverse_raw(raw, inverse_operator, lambda2, dSPM=True, if time_func is not None: data = time_func(data) - K, noise_norm, vertno = _assemble_kernel(inv, label, dSPM, pick_normal) + K, noise_norm, vertno = _assemble_kernel(inv, label, method, pick_normal) is_free_ori = (inverse_operator['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI and not pick_normal) @@ -811,8 +839,8 @@ def apply_inverse_raw(raw, inverse_operator, lambda2, dSPM=True, return stc -def apply_inverse_epochs(epochs, inverse_operator, lambda2, dSPM=True, - label=None, nave=1, pick_normal=False): +def apply_inverse_epochs(epochs, inverse_operator, lambda2, method="dSPM", + label=None, nave=1, pick_normal=False, dSPM=None): """Apply inverse operator to Epochs Computes a L2-norm inverse solution on each epochs and returns @@ -826,8 +854,8 @@ def apply_inverse_epochs(epochs, inverse_operator, lambda2, dSPM=True, Inverse operator read with mne.read_inverse_operator lambda2: float The regularization parameter - dSPM: bool - do dSPM ? + method: "MNE" | "dSPM" | "sLORETA" + Use mininum norm, dSPM or sLORETA label: Label Restricts the source estimates to a given label nave: int @@ -843,12 +871,14 @@ def apply_inverse_epochs(epochs, inverse_operator, lambda2, dSPM=True, stc: list of SourceEstimate The source estimates for all epochs """ + method = _check_method(method, dSPM) + _chech_ch_names(inverse_operator, epochs.info) # # Set up the inverse according to the parameters # - inv = prepare_inverse_operator(inverse_operator, nave, lambda2, dSPM) + inv = prepare_inverse_operator(inverse_operator, nave, lambda2, method) # # Pick the correct channels from the data # @@ -856,7 +886,7 @@ def apply_inverse_epochs(epochs, inverse_operator, lambda2, dSPM=True, print 'Picked %d channels from the data' % len(sel) print 'Computing inverse...', - K, noise_norm, vertno = _assemble_kernel(inv, label, dSPM, pick_normal) + K, noise_norm, vertno = _assemble_kernel(inv, label, method, pick_normal) stcs = list() tstep = 1.0 / epochs.info['sfreq'] diff --git a/mne/minimum_norm/tests/test_inverse.py b/mne/minimum_norm/tests/test_inverse.py index 540f4e3..b396ee0 100644 --- a/mne/minimum_norm/tests/test_inverse.py +++ b/mne/minimum_norm/tests/test_inverse.py @@ -47,7 +47,6 @@ noise_cov = read_cov(fname_cov) raw = fiff.Raw(fname_raw) snr = 3.0 lambda2 = 1.0 / snr ** 2 -dSPM = True def _compare(a, b): @@ -90,17 +89,17 @@ def test_apply_inverse_operator(): """ evoked = fiff.Evoked(fname_data, setno=0, baseline=(None, 0)) - stc = apply_inverse(evoked, inverse_operator, lambda2, dSPM=False) - + stc = apply_inverse(evoked, inverse_operator, lambda2, "MNE") assert_true(stc.data.min() > 0) assert_true(stc.data.max() < 10e-10) assert_true(stc.data.mean() > 1e-11) - stc = apply_inverse(evoked, inverse_operator, lambda2, dSPM=True) - - assert_true(np.all(stc.data > 0)) - assert_true(np.all(stc.data < 35)) + stc = apply_inverse(evoked, inverse_operator, lambda2, "sLORETA") + assert_true(stc.data.min() > 0) + assert_true(stc.data.max() < 9.0) + assert_true(stc.data.mean() > 0.1) + stc = apply_inverse(evoked, inverse_operator, lambda2, "dSPM") assert_true(stc.data.min() > 0) assert_true(stc.data.max() < 35) assert_true(stc.data.mean() > 0.1) @@ -115,7 +114,7 @@ def test_apply_inverse_operator(): read_my_inv_op = read_inverse_operator('test-inv.fif') _compare(my_inv_op, read_my_inv_op) - my_stc = apply_inverse(evoked, my_inv_op, lambda2, dSPM) + my_stc = apply_inverse(evoked, my_inv_op, lambda2, "dSPM") assert_true('dev_head_t' in my_inv_op['info']) assert_true('mri_head_t' in my_inv_op) @@ -141,8 +140,8 @@ def test_make_inverse_operator_fixed(): assert_array_almost_equal(inverse_operator_fixed['source_cov']['data'], inv_op['source_cov']['data']) - stc_fixed = apply_inverse(evoked, inverse_operator_fixed, lambda2, dSPM) - my_stc = apply_inverse(evoked, inv_op, lambda2, dSPM) + stc_fixed = apply_inverse(evoked, inverse_operator_fixed, lambda2, "dSPM") + my_stc = apply_inverse(evoked, inv_op, lambda2, "dSPM") assert_equal(stc_fixed.times, my_stc.times) assert_array_almost_equal(stc_fixed.data, my_stc.data, 2) @@ -157,7 +156,7 @@ def test_inverse_operator_volume(): """Test MNE inverse computation on volume source space""" evoked = fiff.Evoked(fname_data, setno=0, baseline=(None, 0)) inverse_operator_vol = read_inverse_operator(fname_vol_inv) - stc = apply_inverse(evoked, inverse_operator_vol, lambda2, dSPM) + stc = apply_inverse(evoked, inverse_operator_vol, lambda2, "dSPM") stc.save('tmp-vl.stc') stc2 = SourceEstimate('tmp-vl.stc') assert_true(np.all(stc.data > 0)) @@ -172,11 +171,11 @@ def test_apply_mne_inverse_raw(): stop = 10 _, times = raw[0, start:stop] for pick_normal in [False, True]: - stc = apply_inverse_raw(raw, inverse_operator, lambda2, dSPM=True, + stc = apply_inverse_raw(raw, inverse_operator, lambda2, "dSPM", label=label, start=start, stop=stop, nave=1, pick_normal=pick_normal, buffer_size=None) - stc2 = apply_inverse_raw(raw, inverse_operator, lambda2, dSPM=True, + stc2 = apply_inverse_raw(raw, inverse_operator, lambda2, "dSPM", label=label, start=start, stop=stop, nave=1, pick_normal=pick_normal, buffer_size=3) @@ -201,11 +200,11 @@ def test_apply_mne_inverse_fixed_raw(): inv_op = make_inverse_operator(raw.info, fwd, noise_cov, loose=None, depth=0.8) - stc = apply_inverse_raw(raw, inv_op, lambda2, dSPM=True, + stc = apply_inverse_raw(raw, inv_op, lambda2, "dSPM", label=label, start=start, stop=stop, nave=1, pick_normal=False, buffer_size=None) - stc2 = apply_inverse_raw(raw, inv_op, lambda2, dSPM=True, + stc2 = apply_inverse_raw(raw, inv_op, lambda2, "dSPM", label=label, start=start, stop=stop, nave=1, pick_normal=False, buffer_size=3) @@ -227,7 +226,7 @@ def test_apply_mne_inverse_epochs(): events = read_events(fname_event)[:15] epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks, baseline=(None, 0), reject=reject, flat=flat) - stcs = apply_inverse_epochs(epochs, inverse_operator, lambda2, dSPM, + stcs = apply_inverse_epochs(epochs, inverse_operator, lambda2, "dSPM", label=label, pick_normal=True) assert_true(len(stcs) == 4) diff --git a/mne/minimum_norm/time_frequency.py b/mne/minimum_norm/time_frequency.py index 95f2667..ffc421c 100644 --- a/mne/minimum_norm/time_frequency.py +++ b/mne/minimum_norm/time_frequency.py @@ -9,15 +9,15 @@ from ..fiff.constants import FIFF from ..time_frequency.tfr import cwt, morlet from ..baseline import rescale from .inverse import combine_xyz, prepare_inverse_operator, _assemble_kernel, \ - _make_stc, _pick_channels_inverse_operator + _make_stc, _pick_channels_inverse_operator, _check_method from ..parallel import parallel_func def source_band_induced_power(epochs, inverse_operator, bands, label=None, - lambda2=1.0 / 9.0, dSPM=True, n_cycles=5, df=1, + lambda2=1.0 / 9.0, method="dSPM", n_cycles=5, df=1, use_fft=False, decim=1, baseline=None, baseline_mode='logratio', pca=True, - n_jobs=1): + n_jobs=1, dSPM=None): """Compute source space induced power in given frequency bands Parameters @@ -32,8 +32,8 @@ def source_band_induced_power(epochs, inverse_operator, bands, label=None, Restricts the source estimates to a given label lambda2: float The regularization parameter of the minimum norm - dSPM: bool - Do dSPM or not? + method: "MNE" | "dSPM" | "sLORETA" + Use mininum norm, dSPM or sLORETA n_cycles: int Number of cycles df: float @@ -62,13 +62,15 @@ def source_band_induced_power(epochs, inverse_operator, bands, label=None, n_jobs: int Number of jobs to run in parallel """ + method = _check_method(method, dSPM) + frequencies = np.concatenate([np.arange(band[0], band[1] + df / 2.0, df) for _, band in bands.iteritems()]) powers, _, vertno = _source_induced_power(epochs, inverse_operator, frequencies, label=label, - lambda2=lambda2, dSPM=dSPM, + lambda2=lambda2, method=method, n_cycles=n_cycles, decim=decim, use_fft=use_fft, pca=pca, n_jobs=n_jobs, with_plv=False) @@ -162,7 +164,7 @@ def _compute_pow_plv(data, K, sel, Ws, source_ori, use_fft, Vh, with_plv, def _source_induced_power(epochs, inverse_operator, frequencies, label=None, - lambda2=1.0 / 9.0, dSPM=True, n_cycles=5, decim=1, + lambda2=1.0 / 9.0, method="dSPM", n_cycles=5, decim=1, use_fft=False, pca=True, pick_normal=True, n_jobs=1, with_plv=True): """Aux function for source_induced_power @@ -176,7 +178,7 @@ def _source_induced_power(epochs, inverse_operator, frequencies, label=None, nave = len(epochs_data) # XXX : can do better when no preload - inv = prepare_inverse_operator(inverse_operator, nave, lambda2, dSPM) + inv = prepare_inverse_operator(inverse_operator, nave, lambda2, method) # # Pick the correct channels from the data # @@ -190,7 +192,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, noise_norm, vertno = _assemble_kernel(inv, label, dSPM, pick_normal) + K, noise_norm, vertno = _assemble_kernel(inv, label, method, pick_normal) if pca: U, s, Vh = linalg.svd(K, full_matrices=False) @@ -222,16 +224,17 @@ def _source_induced_power(epochs, inverse_operator, frequencies, label=None, else: plv = None - if dSPM: + if method != "MNE": power *= noise_norm.ravel()[:, None, None] ** 2 return power, plv, vertno def source_induced_power(epochs, inverse_operator, frequencies, label=None, - lambda2=1.0 / 9.0, dSPM=True, n_cycles=5, decim=1, + lambda2=1.0 / 9.0, method="dSPM", n_cycles=5, decim=1, use_fft=False, pick_normal=False, baseline=None, - baseline_mode='logratio', pca=True, n_jobs=1): + baseline_mode='logratio', pca=True, n_jobs=1, + dSPM=None): """Compute induced power and phase lock Computation can optionaly be restricted in a label. @@ -248,8 +251,8 @@ def source_induced_power(epochs, inverse_operator, frequencies, label=None, Array of frequencies of interest lambda2: float The regularization parameter of the minimum norm - dSPM: bool - Do dSPM or not? + method: "MNE" | "dSPM" | "sLORETA" + Use mininum norm, dSPM or sLORETA n_cycles: int Number of cycles decim: int @@ -280,9 +283,11 @@ def source_induced_power(epochs, inverse_operator, frequencies, label=None, n_jobs: int Number of jobs to run in parallel """ + method = _check_method(method, dSPM) + power, plv, vertno = _source_induced_power(epochs, inverse_operator, frequencies, - label, lambda2, dSPM, n_cycles, decim, + label, lambda2, method, n_cycles, decim, use_fft, pick_normal=pick_normal, pca=pca, n_jobs=n_jobs) -- 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
