[med-svn] [python-mne] 197/353: ENH + API : add support for sLORETA to apply_inverse, apply_inverse_raw, apply_inverse_epochs
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:24:57 UTC 2015
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 <alexandre.gramfort at inria.fr>
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
More information about the debian-med-commit
mailing list