From 7c53b5b1a802e5bce25209e9d8ba4cf2578c1a1e Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Wed, 24 Apr 2019 19:02:31 +0200 Subject: [PATCH 01/16] WIP: DICS whitening --- mne/beamformer/_dics.py | 26 ++++++++++++++++++-------- mne/time_frequency/csd.py | 17 ++++++++++++++--- mne/utils/check.py | 20 +++++++++++--------- 3 files changed, 43 insertions(+), 20 deletions(-) diff --git a/mne/beamformer/_dics.py b/mne/beamformer/_dics.py index 38dfe86a8e0..27093c198ad 100644 --- a/mne/beamformer/_dics.py +++ b/mne/beamformer/_dics.py @@ -21,8 +21,8 @@ @verbose -def make_dics(info, forward, csd, reg=0.05, label=None, pick_ori=None, - rank=None, inversion='single', weight_norm=None, +def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, + pick_ori=None, rank=None, inversion='single', weight_norm=None, normalize_fwd=True, real_filter=False, reduce_rank=False, verbose=None): """Compute a Dynamic Imaging of Coherent Sources (DICS) spatial filter. @@ -48,6 +48,12 @@ def make_dics(info, forward, csd, reg=0.05, label=None, pick_ori=None, reg : float The regularization to apply to the cross-spectral density before computing the inverse. + noise_csd : instance of CrossSpectralDensity | None + Noise cross-spectral density (CSD) matrices. If provided, whitening + will be done. The noise CSDs need to have been computed for the same + frequencies as the data CSDs. Providing noise CSDs is mandatory if you + mix sensor types, e.g. gradiometers with magnetometers or EEG with + MEG. label : Label | None Restricts the solution to a given label. pick_ori : None | 'normal' | 'max-power' @@ -218,15 +224,10 @@ def make_dics(info, forward, csd, reg=0.05, label=None, pick_ori=None, exp = None # turn off depth weighting entirely combine_xyz = False - _check_one_ch_type('dics', info, forward) + _check_one_ch_type('dics', info, forward, csd, noise_csd) # pick info, get gain matrix, etc. - _, info, proj, vertices, G, _, nn, orient_std = _prepare_beamformer_input( - info, forward, label, pick_ori, - combine_xyz=combine_xyz, exp=exp) subject = _subject_from_forward(forward) - src_type = _get_src_type(forward['src'], vertices) - del forward ch_names = list(info['ch_names']) csd_picks = [csd.ch_names.index(ch) for ch in ch_names] @@ -239,6 +240,14 @@ def make_dics(info, forward, csd, reg=0.05, label=None, pick_ori=None, (freq, i + 1, n_freqs)) Cm = csd.get_data(index=i) + if noise_csd is not None: + noise_csd_freq = noise_csd.get_data(index=i, as_cov=True) + else: + noise_csd_freq = None + + _, _, proj, vertices, G, _, nn, orient_std = _prepare_beamformer_input( + info, forward, label, pick_ori, noise_cov=noise_csd_freq, + combine_xyz=combine_xyz, exp=exp) if real_filter: Cm = Cm.real @@ -254,6 +263,7 @@ def make_dics(info, forward, csd, reg=0.05, label=None, pick_ori=None, Ws = np.array(Ws) + src_type = _get_src_type(forward['src'], vertices) filters = Beamformer( kind='DICS', weights=Ws, csd=csd, ch_names=ch_names, proj=proj, vertices=vertices, subject=subject, pick_ori=pick_ori, diff --git a/mne/time_frequency/csd.py b/mne/time_frequency/csd.py index 068a4a70d78..f91717ef5b4 100644 --- a/mne/time_frequency/csd.py +++ b/mne/time_frequency/csd.py @@ -295,7 +295,7 @@ def pick_frequency(self, freq=None, index=None): return self[index] - def get_data(self, frequency=None, index=None): + def get_data(self, frequency=None, index=None, as_cov=False): """Get the CSD matrix for a given frequency as NumPy array. If there is only one matrix defined in the CSD object, calling this @@ -311,10 +311,13 @@ def get_data(self, frequency=None, index=None): index : int | None Return the CSD matrix for the frequency or frequency-bin with the given index. + as_cov : bool + Whether to return the data as a numpy array (`False`, the default), + or pack it in a :class:`mne.Covariance` object (`True`). Returns ------- - csd : ndarray, shape (n_channels, n_channels) + csd : ndarray, shape (n_channels, n_channels) | instance of Covariance The CSD matrix corresponding to the requested frequency. See Also @@ -332,7 +335,15 @@ def get_data(self, frequency=None, index=None): raise ValueError('Cannot specify both a frequency and index.') index = self._get_frequency_index(frequency) - return _vector_to_sym_mat(self._data[:, index]) + data = _vector_to_sym_mat(self._data[:, index]) + if as_cov: + # Pack the data into a Covariance object + from ..cov import Covariance # to avoid circular import + return Covariance(data, self.ch_names, bads=[], projs=self.projs, + nfree=self.n_fft) + return None + else: + return data @copy_function_doc_to_method_doc(plot_csd) def plot(self, info=None, mode='csd', colorbar=True, cmap='viridis', diff --git a/mne/utils/check.py b/mne/utils/check.py index db65d10bcc0..25331fbba10 100644 --- a/mne/utils/check.py +++ b/mne/utils/check.py @@ -446,29 +446,31 @@ def _check_rank(rank): def _check_one_ch_type(method, info, forward, data_cov=None, noise_cov=None): """Check number of sensor types and presence of noise covariance matrix.""" from ..cov import make_ad_hoc_cov, Covariance + from ..time_frequency.csd import CrossSpectralDensity from ..io.pick import pick_info from ..channels.channels import _contains_ch_type - picks = _check_info_inv(info, forward, data_cov=data_cov, - noise_cov=noise_cov) - info_pick = pick_info(info, picks) + if isinstance(data_cov, CrossSpectralDensity): + # FIXME + picks = list(range(len(data_cov.ch_names))) + info_pick = info + else: + _validate_type(noise_cov, Covariance, 'noise_cov') + picks = _check_info_inv(info, forward, data_cov=data_cov, + noise_cov=noise_cov) + info_pick = pick_info(info, picks) ch_types =\ [_contains_ch_type(info_pick, tt) for tt in ('mag', 'grad', 'eeg')] if sum(ch_types) > 1: - if method == 'lcmv' and noise_cov is None: + if noise_cov is None: raise ValueError('Source reconstruction with several sensor types' ' requires a noise covariance matrix to be ' 'able to apply whitening.') - if method == 'dics': - raise RuntimeError( - 'The use of several sensor types with the DICS beamformer is ' - 'not supported yet.') if noise_cov is None: noise_cov = make_ad_hoc_cov(info_pick, std=1.) else: noise_cov = noise_cov.copy() if 'estimator' in noise_cov: del noise_cov['estimator'] - _validate_type(noise_cov, Covariance, 'noise_cov') return noise_cov, picks From 34a1578f717bf3d8f277b48fc78b32b8827af901 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Tue, 5 Nov 2019 10:29:05 +0200 Subject: [PATCH 02/16] Hacky way of making things work --- examples/inverse/plot_dics_source_power.py | 9 ++-- mne/beamformer/_dics.py | 10 ++-- mne/io/pick.py | 57 ++++++++++++++++++++++ mne/io/tests/test_pick.py | 20 +++++++- mne/time_frequency/csd.py | 5 +- mne/utils/check.py | 3 +- 6 files changed, 91 insertions(+), 13 deletions(-) diff --git a/examples/inverse/plot_dics_source_power.py b/examples/inverse/plot_dics_source_power.py index 0e6704ee27f..10c97cb15df 100644 --- a/examples/inverse/plot_dics_source_power.py +++ b/examples/inverse/plot_dics_source_power.py @@ -41,13 +41,9 @@ raw = mne.io.read_raw_fif(raw_fname) -# Set picks, use a single sensor type -picks = mne.pick_types(raw.info, meg='grad', exclude='bads') - # Read epochs events = mne.find_events(raw) -epochs = mne.Epochs(raw, events, event_id=1, tmin=-1.5, tmax=2, picks=picks, - preload=True) +epochs = mne.Epochs(raw, events, event_id=1, tmin=-1.5, tmax=2, preload=True) # Read forward operator and point to freesurfer subject directory fname_fwd = op.join(data_path, 'derivatives', 'sub-{}'.format(subject), @@ -73,7 +69,8 @@ ############################################################################### # Computing DICS spatial filters using the CSD that was computed on the entire # timecourse. -filters = make_dics(epochs.info, fwd, csd.mean(), pick_ori='max-power') +filters = make_dics(epochs.info, fwd, csd.mean(), noise_csd=csd_baseline.mean(), + pick_ori='max-power') ############################################################################### # Applying DICS spatial filters separately to the CSD computed using the diff --git a/mne/beamformer/_dics.py b/mne/beamformer/_dics.py index 27093c198ad..ae54b189ba2 100644 --- a/mne/beamformer/_dics.py +++ b/mne/beamformer/_dics.py @@ -8,9 +8,10 @@ # License: BSD (3-clause) import numpy as np +from ..io.pick import pick_info from ..utils import (logger, verbose, warn, _check_one_ch_type, _check_channels_spatial_filter, _check_rank, - _check_option) + _check_option, _check_info_inv) from ..forward import _subject_from_forward from ..minimum_norm.inverse import combine_xyz, _check_reference from ..source_estimate import _make_stc, _get_src_type @@ -19,7 +20,6 @@ _compute_beamformer, _check_src_type, Beamformer, _compute_power) - @verbose def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, pick_ori=None, rank=None, inversion='single', weight_norm=None, @@ -192,6 +192,9 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, _check_option('inversion', inversion, ['single', 'matrix']) _check_option('weight_norm', weight_norm, ['unit-noise-gain', 'nai', None]) + picks = _check_info_inv(info, forward, csd, noise_csd) + info = pick_info(info, picks) + # Leadfield rank and optional rank reduction # (to deal with problems with complex eigenvalues within the computation # of the optimal orientation when using pinv if the leadfield was only @@ -230,7 +233,8 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, subject = _subject_from_forward(forward) ch_names = list(info['ch_names']) - csd_picks = [csd.ch_names.index(ch) for ch in ch_names] + csd_picks = pick_types(info, meg=True, eeg=True, seeg=True) + #csd_picks = [csd.ch_names.index(ch) for ch in ch_names] logger.info('Computing DICS spatial filters...') Ws = [] diff --git a/mne/io/pick.py b/mne/io/pick.py index 2f5f45bad96..3702a57b3ab 100644 --- a/mne/io/pick.py +++ b/mne/io/pick.py @@ -10,6 +10,9 @@ import numpy as np +from ..base import BaseRaw +from ..epochs import BaseEpochs +from ..io import Info from .constants import FIFF from ..utils import (logger, verbose, _validate_type, fill_doc, _ensure_int, _check_option) @@ -1095,3 +1098,57 @@ def _get_channel_types(info, picks=None, unique=True, ch_types = [ch_type for ch_type in ch_types if ch_type in _DATA_CH_TYPES_SPLIT] return set(ch_types) if unique is True else ch_types + + +@verbose +def pick_channels_intersection(instances, verbose=None): + """Pick channels present in all instances. + + Goes through all given instances and picks the channels that are present in + all of them. + + Parameters + ---------- + instances : list of (inst | list of str) + List of instances that define a list of channel names. Each entry can + be one of:: + - :class:`AverageTFR` + - :class:`Covariance` + - :class:`CrossSpectralDensity` + - :class:`Epochs` + - :class:`EpochsTFR` + - :class:`Evoked` + - :class:`Forward` + - :class:`Info` + - :class:`Raw` + - List of string channel names + %(verbose)s + + Returns + ------- + ch_sel : list of str + A list of channel names that are present in all given instances. + """ + ch_sel = None + for inst in instances: + if isinstance(inst, (BaseRaw, BaseEpochs)): + channels = inst.ch_names + elif isinstance(inst, Info): + channels = inst['ch_names'] + elif isinstance(inst, (list, tuple, set, np.ndarray)): + if len(inst) > 0 and not isinstance(inst[0], str): + raise TypeError('When giving a list of channel names, the ' + ' elements of the list must be strings.') + channels = inst + + if ch_sel is None: + ch_sel = set(channels) + else: + ch_sel &= set(channels) + + if ch_sel is None: + raise TypeError('No instances define a list of channel names.') + elif len(ch_sel) == 0: + raise ValueError('Given instances have no channels in common.') + + return list(ch_sel) diff --git a/mne/io/tests/test_pick.py b/mne/io/tests/test_pick.py index 8ba2bdbeda9..50d17205984 100644 --- a/mne/io/tests/test_pick.py +++ b/mne/io/tests/test_pick.py @@ -6,7 +6,7 @@ import pytest import numpy as np -from mne import (pick_channels_regexp, pick_types, Epochs, +from mne import (pick_channels_regexp, pick_types, Epochs, EpochsArray, read_forward_solution, rename_channels, pick_info, pick_channels, create_info) from mne import __file__ as _root_init_fname @@ -15,7 +15,7 @@ from mne.io.pick import (channel_indices_by_type, channel_type, pick_types_forward, _picks_by_type, _picks_to_idx, get_channel_types, _DATA_CH_TYPES_SPLIT, - _contains_ch_type) + _contains_ch_type, pick_channels_intersection) from mne.io.constants import FIFF from mne.datasets import testing from mne.utils import run_tests_if_main, catch_logging, assert_object_equal @@ -507,4 +507,20 @@ def test_picks_to_idx(): assert_array_equal([0], _picks_to_idx(info, 'data')) +def test_pick_channels_intersection(): + raw = RawArray(np.zeros((1, 1)), create_info(['raw', 'CH2', 'CH3'], 1024.)) + epochs = EpochsArray(np.zeros((1, 1, 1)), + create_info(['CH1', 'info', 'CH2'], 1024.)) + info = create_info(['CH1', 'CH2', 'info'], 1024.) + + pytest.raises(TypeError, pick_channels_intersection, 1) + pytest.raises(TypeError, pick_channels_intersection, [1]) + pytest.raises(ValueError, pick_channels_intersection, []) + pytest.raises(TypeError, pick_channels_intersection, [raw, [1]]) + pytest.raises(ValueError, pick_channels_intersection, [raw, ['nooverlap']]) + assert pick_channels_intersection([raw, ['CH1', 'CH2']]) == ['CH2'] + assert pick_channels_intersection([raw, epochs]) == ['CH2'] + assert pick_channels_intersection([raw, epochs, info]) == ['CH2'] + + run_tests_if_main() diff --git a/mne/time_frequency/csd.py b/mne/time_frequency/csd.py index f91717ef5b4..ba8897d3b70 100644 --- a/mne/time_frequency/csd.py +++ b/mne/time_frequency/csd.py @@ -92,7 +92,10 @@ def __init__(self, data, ch_names, frequencies, n_fft, tmin=None, self.frequencies = frequencies self.n_fft = n_fft - self.projs = cp.deepcopy(projs) + if projs is None: + self.projs = [] + else: + self.projs = cp.deepcopy(projs) @property def n_channels(self): diff --git a/mne/utils/check.py b/mne/utils/check.py index 25331fbba10..48944db11b0 100644 --- a/mne/utils/check.py +++ b/mne/utils/check.py @@ -451,6 +451,7 @@ def _check_one_ch_type(method, info, forward, data_cov=None, noise_cov=None): from ..channels.channels import _contains_ch_type if isinstance(data_cov, CrossSpectralDensity): # FIXME + _validate_type(noise_cov, [None, CrossSpectralDensity], 'noise_cov') picks = list(range(len(data_cov.ch_names))) info_pick = info else: @@ -469,7 +470,7 @@ def _check_one_ch_type(method, info, forward, data_cov=None, noise_cov=None): noise_cov = make_ad_hoc_cov(info_pick, std=1.) else: noise_cov = noise_cov.copy() - if 'estimator' in noise_cov: + if isinstance(noise_cov, Covariance) and 'estimator' in noise_cov: del noise_cov['estimator'] return noise_cov, picks From 56320f4ea8e5e4aafeed423055856a3c5bc2ee32 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Tue, 5 Nov 2019 10:40:26 +0200 Subject: [PATCH 03/16] Undo broken equalize_channels implementation --- mne/io/pick.py | 57 --------------------------------------- mne/io/tests/test_pick.py | 18 +------------ 2 files changed, 1 insertion(+), 74 deletions(-) diff --git a/mne/io/pick.py b/mne/io/pick.py index 3702a57b3ab..2f5f45bad96 100644 --- a/mne/io/pick.py +++ b/mne/io/pick.py @@ -10,9 +10,6 @@ import numpy as np -from ..base import BaseRaw -from ..epochs import BaseEpochs -from ..io import Info from .constants import FIFF from ..utils import (logger, verbose, _validate_type, fill_doc, _ensure_int, _check_option) @@ -1098,57 +1095,3 @@ def _get_channel_types(info, picks=None, unique=True, ch_types = [ch_type for ch_type in ch_types if ch_type in _DATA_CH_TYPES_SPLIT] return set(ch_types) if unique is True else ch_types - - -@verbose -def pick_channels_intersection(instances, verbose=None): - """Pick channels present in all instances. - - Goes through all given instances and picks the channels that are present in - all of them. - - Parameters - ---------- - instances : list of (inst | list of str) - List of instances that define a list of channel names. Each entry can - be one of:: - - :class:`AverageTFR` - - :class:`Covariance` - - :class:`CrossSpectralDensity` - - :class:`Epochs` - - :class:`EpochsTFR` - - :class:`Evoked` - - :class:`Forward` - - :class:`Info` - - :class:`Raw` - - List of string channel names - %(verbose)s - - Returns - ------- - ch_sel : list of str - A list of channel names that are present in all given instances. - """ - ch_sel = None - for inst in instances: - if isinstance(inst, (BaseRaw, BaseEpochs)): - channels = inst.ch_names - elif isinstance(inst, Info): - channels = inst['ch_names'] - elif isinstance(inst, (list, tuple, set, np.ndarray)): - if len(inst) > 0 and not isinstance(inst[0], str): - raise TypeError('When giving a list of channel names, the ' - ' elements of the list must be strings.') - channels = inst - - if ch_sel is None: - ch_sel = set(channels) - else: - ch_sel &= set(channels) - - if ch_sel is None: - raise TypeError('No instances define a list of channel names.') - elif len(ch_sel) == 0: - raise ValueError('Given instances have no channels in common.') - - return list(ch_sel) diff --git a/mne/io/tests/test_pick.py b/mne/io/tests/test_pick.py index 50d17205984..e7d2e848d96 100644 --- a/mne/io/tests/test_pick.py +++ b/mne/io/tests/test_pick.py @@ -15,7 +15,7 @@ from mne.io.pick import (channel_indices_by_type, channel_type, pick_types_forward, _picks_by_type, _picks_to_idx, get_channel_types, _DATA_CH_TYPES_SPLIT, - _contains_ch_type, pick_channels_intersection) + _contains_ch_type) from mne.io.constants import FIFF from mne.datasets import testing from mne.utils import run_tests_if_main, catch_logging, assert_object_equal @@ -507,20 +507,4 @@ def test_picks_to_idx(): assert_array_equal([0], _picks_to_idx(info, 'data')) -def test_pick_channels_intersection(): - raw = RawArray(np.zeros((1, 1)), create_info(['raw', 'CH2', 'CH3'], 1024.)) - epochs = EpochsArray(np.zeros((1, 1, 1)), - create_info(['CH1', 'info', 'CH2'], 1024.)) - info = create_info(['CH1', 'CH2', 'info'], 1024.) - - pytest.raises(TypeError, pick_channels_intersection, 1) - pytest.raises(TypeError, pick_channels_intersection, [1]) - pytest.raises(ValueError, pick_channels_intersection, []) - pytest.raises(TypeError, pick_channels_intersection, [raw, [1]]) - pytest.raises(ValueError, pick_channels_intersection, [raw, ['nooverlap']]) - assert pick_channels_intersection([raw, ['CH1', 'CH2']]) == ['CH2'] - assert pick_channels_intersection([raw, epochs]) == ['CH2'] - assert pick_channels_intersection([raw, epochs, info]) == ['CH2'] - - run_tests_if_main() From 1e65b1948d71ec1d602fd23f81c7ce771ce2dfa4 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Wed, 6 Nov 2019 15:50:06 +0200 Subject: [PATCH 04/16] Update mne/beamformer/_dics.py Co-Authored-By: Eric Larson --- mne/beamformer/_dics.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/mne/beamformer/_dics.py b/mne/beamformer/_dics.py index ae54b189ba2..43ab9c6e03f 100644 --- a/mne/beamformer/_dics.py +++ b/mne/beamformer/_dics.py @@ -54,6 +54,8 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, frequencies as the data CSDs. Providing noise CSDs is mandatory if you mix sensor types, e.g. gradiometers with magnetometers or EEG with MEG. + + .. versionadded:: 0.20 label : Label | None Restricts the solution to a given label. pick_ori : None | 'normal' | 'max-power' From 275f68ba4b7bfbb08afd461e98b780a4e29e5d77 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Wed, 6 Nov 2019 15:50:15 +0200 Subject: [PATCH 05/16] Update mne/time_frequency/csd.py Co-Authored-By: Eric Larson --- mne/time_frequency/csd.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/mne/time_frequency/csd.py b/mne/time_frequency/csd.py index ba8897d3b70..942d8c5c7a0 100644 --- a/mne/time_frequency/csd.py +++ b/mne/time_frequency/csd.py @@ -318,6 +318,8 @@ def get_data(self, frequency=None, index=None, as_cov=False): Whether to return the data as a numpy array (`False`, the default), or pack it in a :class:`mne.Covariance` object (`True`). + .. versionadded:: 0.20 + Returns ------- csd : ndarray, shape (n_channels, n_channels) | instance of Covariance From 46efd8d8ee406dade8ca9ecdcbdc4675026ac2d4 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Wed, 6 Nov 2019 16:40:59 +0200 Subject: [PATCH 06/16] remove stray line --- mne/time_frequency/csd.py | 1 - 1 file changed, 1 deletion(-) diff --git a/mne/time_frequency/csd.py b/mne/time_frequency/csd.py index 942d8c5c7a0..1b2c1b56d52 100644 --- a/mne/time_frequency/csd.py +++ b/mne/time_frequency/csd.py @@ -346,7 +346,6 @@ def get_data(self, frequency=None, index=None, as_cov=False): from ..cov import Covariance # to avoid circular import return Covariance(data, self.ch_names, bads=[], projs=self.projs, nfree=self.n_fft) - return None else: return data From da95615a0d42d2cfa9992d44c13c135b8c021a0b Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Thu, 7 Nov 2019 15:11:38 +0200 Subject: [PATCH 07/16] Make DICS example work again, now with whitening --- examples/inverse/plot_dics_source_power.py | 13 ++++++++++--- mne/beamformer/_dics.py | 6 +++--- 2 files changed, 13 insertions(+), 6 deletions(-) diff --git a/examples/inverse/plot_dics_source_power.py b/examples/inverse/plot_dics_source_power.py index 10c97cb15df..ea47975be45 100644 --- a/examples/inverse/plot_dics_source_power.py +++ b/examples/inverse/plot_dics_source_power.py @@ -66,17 +66,24 @@ # ERS activity starts at 0.5 seconds after stimulus onset csd_ers = csd_morlet(epochs, freqs, tmin=0.5, tmax=1.5, decim=20) +############################################################################### +# To compute the source power for a frequency band, rather than each frequency +# separately, we average the CSD objects across frequencies. +csd = csd.mean() +csd_baseline = csd_baseline.mean() +csd_ers = csd_ers.mean() + ############################################################################### # Computing DICS spatial filters using the CSD that was computed on the entire # timecourse. -filters = make_dics(epochs.info, fwd, csd.mean(), noise_csd=csd_baseline.mean(), +filters = make_dics(epochs.info, fwd, csd, noise_csd=csd_baseline, pick_ori='max-power') ############################################################################### # Applying DICS spatial filters separately to the CSD computed using the # baseline and the CSD computed during the ERS activity. -baseline_source_power, freqs = apply_dics_csd(csd_baseline.mean(), filters) -beta_source_power, freqs = apply_dics_csd(csd_ers.mean(), filters) +baseline_source_power, freqs = apply_dics_csd(csd_baseline, filters) +beta_source_power, freqs = apply_dics_csd(csd_ers, filters) ############################################################################### # Visualizing source power during ERS activity relative to the baseline power. diff --git a/mne/beamformer/_dics.py b/mne/beamformer/_dics.py index 43ab9c6e03f..ff043ca579c 100644 --- a/mne/beamformer/_dics.py +++ b/mne/beamformer/_dics.py @@ -8,7 +8,7 @@ # License: BSD (3-clause) import numpy as np -from ..io.pick import pick_info +from ..io.pick import pick_info, pick_types from ..utils import (logger, verbose, warn, _check_one_ch_type, _check_channels_spatial_filter, _check_rank, _check_option, _check_info_inv) @@ -52,7 +52,7 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, Noise cross-spectral density (CSD) matrices. If provided, whitening will be done. The noise CSDs need to have been computed for the same frequencies as the data CSDs. Providing noise CSDs is mandatory if you - mix sensor types, e.g. gradiometers with magnetometers or EEG with + mix sensor types, e.g. gradiometers with magnetometers or EEG with MEG. .. versionadded:: 0.20 @@ -194,7 +194,7 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, _check_option('inversion', inversion, ['single', 'matrix']) _check_option('weight_norm', weight_norm, ['unit-noise-gain', 'nai', None]) - picks = _check_info_inv(info, forward, csd, noise_csd) + picks = _check_info_inv(info, forward) info = pick_info(info, picks) # Leadfield rank and optional rank reduction From 99650588f679d50adc510e3bda864d74982a83a0 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Thu, 7 Nov 2019 16:54:19 +0200 Subject: [PATCH 08/16] Also whiten the CSD --- examples/inverse/plot_dics_source_power.py | 4 +++- mne/beamformer/_dics.py | 14 ++++++++++---- mne/cov.py | 16 ++++++++-------- mne/minimum_norm/inverse.py | 2 +- mne/utils/check.py | 2 +- 5 files changed, 23 insertions(+), 15 deletions(-) diff --git a/examples/inverse/plot_dics_source_power.py b/examples/inverse/plot_dics_source_power.py index ea47975be45..1ed7ec85439 100644 --- a/examples/inverse/plot_dics_source_power.py +++ b/examples/inverse/plot_dics_source_power.py @@ -45,6 +45,8 @@ events = mne.find_events(raw) epochs = mne.Epochs(raw, events, event_id=1, tmin=-1.5, tmax=2, preload=True) +epochs.pick_types(meg='grad') + # Read forward operator and point to freesurfer subject directory fname_fwd = op.join(data_path, 'derivatives', 'sub-{}'.format(subject), 'sub-{}_task-{}-fwd.fif'.format(subject, task)) @@ -77,7 +79,7 @@ # Computing DICS spatial filters using the CSD that was computed on the entire # timecourse. filters = make_dics(epochs.info, fwd, csd, noise_csd=csd_baseline, - pick_ori='max-power') + pick_ori='max-power', real_filter=True) ############################################################################### # Applying DICS spatial filters separately to the CSD computed using the diff --git a/mne/beamformer/_dics.py b/mne/beamformer/_dics.py index ff043ca579c..b6ffdd0bd30 100644 --- a/mne/beamformer/_dics.py +++ b/mne/beamformer/_dics.py @@ -246,21 +246,27 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, (freq, i + 1, n_freqs)) Cm = csd.get_data(index=i) + if real_filter: + Cm = Cm.real + if noise_csd is not None: noise_csd_freq = noise_csd.get_data(index=i, as_cov=True) + if real_filter: + noise_csd_freq['data'] = noise_csd_freq['data'].real else: noise_csd_freq = None - _, _, proj, vertices, G, _, nn, orient_std = _prepare_beamformer_input( - info, forward, label, pick_ori, noise_cov=noise_csd_freq, + _, _, proj, vertices, G, whitener, nn, orient_std = _prepare_beamformer_input( + info, forward, label, pick_ori, noise_csd_freq, combine_xyz=combine_xyz, exp=exp) - if real_filter: - Cm = Cm.real # Ensure the CSD is in the same order as the leadfield Cm = Cm[csd_picks, :][:, csd_picks] + # Whiten the CSD + Cm = np.dot(whitener, np.dot(Cm, whitener.conj().T)) + # compute spatial filter W = _compute_beamformer(G, Cm, reg, n_orient, weight_norm, pick_ori, reduce_rank, rank=rank, inversion=inversion, diff --git a/mne/cov.py b/mne/cov.py index 3f56ccbf289..a02338490db 100644 --- a/mne/cov.py +++ b/mne/cov.py @@ -891,9 +891,9 @@ def _check_scalings_user(scalings): def _eigvec_subspace(eig, eigvec, mask): """Compute the subspace from a subset of eigenvectors.""" # We do the same thing we do with projectors: - P = np.eye(len(eigvec)) - np.dot(eigvec[~mask].T, eigvec[~mask]) + P = np.eye(len(eigvec)) - np.dot(eigvec[~mask].conj().T, eigvec[~mask]) eig, eigvec = eigh(P) - eigvec = eigvec.T + eigvec = eigvec.conj().T return eig, eigvec @@ -1292,7 +1292,7 @@ def _get_ch_whitener(A, pca, ch_type, rank): """Get whitener params for a set of channels.""" # whitening operator eig, eigvec = eigh(A, overwrite_a=True) - eigvec = eigvec.T + eigvec = eigvec.conj().T mask = np.ones(len(eig), bool) eig[:-rank] = 0.0 mask[:-rank] = False @@ -1399,8 +1399,8 @@ def _smart_eigh(C, info, rank, scalings=None, projs=None, # time saving short-circuit if proj_subspace and sum(rank.values()) == C.shape[0]: return np.ones(n_chan), np.eye(n_chan), np.ones(n_chan, bool) - eig = np.zeros(n_chan) - eigvec = np.zeros((n_chan, n_chan)) + eig = np.zeros(n_chan, dtype=C.dtype) + eigvec = np.zeros((n_chan, n_chan), dtype=C.dtype) mask = np.zeros(n_chan, bool) for ch_type, picks in _picks_by_type(info, meg_combined=True, ref_meg=False, exclude='bads'): @@ -1717,11 +1717,11 @@ def compute_whitener(noise_cov, info=None, picks=None, rank=None, nzero = (eig > 0) eig[~nzero] = 0. # get rid of numerical noise (negative) ones - W = np.zeros((n_chan, 1), dtype=np.float) + W = np.zeros((n_chan, 1), dtype=noise_cov['eigvec'].dtype) W[nzero, 0] = 1.0 / np.sqrt(eig[nzero]) # Rows of eigvec are the eigenvectors W = W * noise_cov['eigvec'] # C ** -0.5 - C = np.sqrt(eig) * noise_cov['eigvec'].T # C ** 0.5 + C = np.sqrt(eig) * noise_cov['eigvec'].conj().T # C ** 0.5 n_nzero = nzero.sum() logger.info(' Created the whitener using a noise covariance matrix ' 'with rank %d (%d small eigenvalues omitted)' @@ -1732,7 +1732,7 @@ def compute_whitener(noise_cov, info=None, picks=None, rank=None, W = W[nzero] C = C[:, nzero] elif pca is False: - W = np.dot(noise_cov['eigvec'].T, W) + W = np.dot(noise_cov['eigvec'].conj().T, W) C = np.dot(C, noise_cov['eigvec']) # Triage return diff --git a/mne/minimum_norm/inverse.py b/mne/minimum_norm/inverse.py index c994cf508df..c2645606df4 100644 --- a/mne/minimum_norm/inverse.py +++ b/mne/minimum_norm/inverse.py @@ -1379,7 +1379,7 @@ def _prepare_forward(forward, info, noise_cov, fixed, loose, rank, pca, noise_cov, info, info_picked['ch_names'], rank) whitener, _ = compute_whitener( noise_cov, info, info_picked['ch_names'], pca=pca, verbose=False) - gain = np.dot(whitener, forward['sol']['data']) + gain = np.dot(whitener.real, forward['sol']['data']) logger.info('Creating the source covariance matrix') source_std = np.ones(gain.shape[1]) diff --git a/mne/utils/check.py b/mne/utils/check.py index 48944db11b0..2cee1e7a595 100644 --- a/mne/utils/check.py +++ b/mne/utils/check.py @@ -450,8 +450,8 @@ def _check_one_ch_type(method, info, forward, data_cov=None, noise_cov=None): from ..io.pick import pick_info from ..channels.channels import _contains_ch_type if isinstance(data_cov, CrossSpectralDensity): - # FIXME _validate_type(noise_cov, [None, CrossSpectralDensity], 'noise_cov') + # FIXME picks = list(range(len(data_cov.ch_names))) info_pick = info else: From 444ec4c09e7b947b18124d405f389afc335d622f Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Thu, 7 Nov 2019 16:54:19 +0200 Subject: [PATCH 09/16] Also whiten the CSD --- examples/inverse/plot_dics_source_power.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/examples/inverse/plot_dics_source_power.py b/examples/inverse/plot_dics_source_power.py index 1ed7ec85439..27e94ae06eb 100644 --- a/examples/inverse/plot_dics_source_power.py +++ b/examples/inverse/plot_dics_source_power.py @@ -45,8 +45,6 @@ events = mne.find_events(raw) epochs = mne.Epochs(raw, events, event_id=1, tmin=-1.5, tmax=2, preload=True) -epochs.pick_types(meg='grad') - # Read forward operator and point to freesurfer subject directory fname_fwd = op.join(data_path, 'derivatives', 'sub-{}'.format(subject), 'sub-{}_task-{}-fwd.fif'.format(subject, task)) From 0dc4b004c164690315c1f43f517b7f448129169f Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Mon, 11 Nov 2019 09:58:49 +0200 Subject: [PATCH 10/16] Make unit tests pass again --- mne/beamformer/_dics.py | 4 ++-- mne/beamformer/tests/test_dics.py | 8 ++++---- mne/utils/check.py | 2 +- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/mne/beamformer/_dics.py b/mne/beamformer/_dics.py index b6ffdd0bd30..abf7c267f7e 100644 --- a/mne/beamformer/_dics.py +++ b/mne/beamformer/_dics.py @@ -235,8 +235,8 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, subject = _subject_from_forward(forward) ch_names = list(info['ch_names']) - csd_picks = pick_types(info, meg=True, eeg=True, seeg=True) - #csd_picks = [csd.ch_names.index(ch) for ch in ch_names] + #csd_picks = pick_types(info, meg=True, eeg=True, seeg=True) + csd_picks = [csd.ch_names.index(ch) for ch in ch_names] logger.info('Computing DICS spatial filters...') Ws = [] diff --git a/mne/beamformer/tests/test_dics.py b/mne/beamformer/tests/test_dics.py index 6def11ac09b..2ea0b29828a 100644 --- a/mne/beamformer/tests/test_dics.py +++ b/mne/beamformer/tests/test_dics.py @@ -114,7 +114,7 @@ def test_make_dics(tmpdir, _load_forward): fwd_free, fwd_surf, fwd_fixed, fwd_vol, label = _load_forward epochs, _, csd, _ = _simulate_data(fwd_fixed) - with pytest.raises(RuntimeError, match='several sensor types'): + with pytest.raises(ValueError, match='several sensor types'): make_dics(epochs.info, fwd_surf, csd, label=label, pick_ori=None) epochs.pick_types(meg='grad') @@ -257,7 +257,7 @@ def test_apply_dics_csd(_load_forward): source_ind = vertices.tolist().index(source_vertno) reg = 1 # Lots of regularization for our toy dataset - with pytest.raises(RuntimeError, match='several sensor types'): + with pytest.raises(ValueError, match='several sensor types'): make_dics(epochs.info, fwd_free, csd) epochs.pick_types(meg='grad') @@ -385,7 +385,7 @@ def test_apply_dics_timeseries(_load_forward): source_ind = vertices.tolist().index(source_vertno) reg = 5 # Lots of regularization for our toy dataset - with pytest.raises(RuntimeError, match='several sensor types'): + with pytest.raises(ValueError, match='several sensor types'): make_dics(evoked.info, fwd_surf, csd) evoked.pick_types(meg='grad') @@ -487,7 +487,7 @@ def test_tf_dics(_load_forward): frequencies = [10, 20] freq_bins = [(8, 12), (18, 22)] - with pytest.raises(RuntimeError, match='several sensor types'): + with pytest.raises(ValueError, match='several sensor types'): stcs = tf_dics(epochs, fwd_surf, None, tmin, tmax, tstep, win_lengths, freq_bins=freq_bins, frequencies=frequencies, decim=10, reg=reg, label=label) diff --git a/mne/utils/check.py b/mne/utils/check.py index 2cee1e7a595..1a03c94e295 100644 --- a/mne/utils/check.py +++ b/mne/utils/check.py @@ -455,7 +455,7 @@ def _check_one_ch_type(method, info, forward, data_cov=None, noise_cov=None): picks = list(range(len(data_cov.ch_names))) info_pick = info else: - _validate_type(noise_cov, Covariance, 'noise_cov') + _validate_type(noise_cov, [None, Covariance], 'noise_cov') picks = _check_info_inv(info, forward, data_cov=data_cov, noise_cov=noise_cov) info_pick = pick_info(info, picks) From b1120816ac6d0448bc08883b41447460061b95f6 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Mon, 11 Nov 2019 11:32:09 +0200 Subject: [PATCH 11/16] Undo --- mne/minimum_norm/inverse.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mne/minimum_norm/inverse.py b/mne/minimum_norm/inverse.py index c2645606df4..c994cf508df 100644 --- a/mne/minimum_norm/inverse.py +++ b/mne/minimum_norm/inverse.py @@ -1379,7 +1379,7 @@ def _prepare_forward(forward, info, noise_cov, fixed, loose, rank, pca, noise_cov, info, info_picked['ch_names'], rank) whitener, _ = compute_whitener( noise_cov, info, info_picked['ch_names'], pca=pca, verbose=False) - gain = np.dot(whitener.real, forward['sol']['data']) + gain = np.dot(whitener, forward['sol']['data']) logger.info('Creating the source covariance matrix') source_std = np.ones(gain.shape[1]) From 6556ff6457f4056c6bd8fdc5ec8b7c48384c6d92 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Mon, 11 Nov 2019 11:56:48 +0200 Subject: [PATCH 12/16] Be careful not to use a 32-bit value as dtype --- mne/cov.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/mne/cov.py b/mne/cov.py index a02338490db..5146e6d09a2 100644 --- a/mne/cov.py +++ b/mne/cov.py @@ -1399,8 +1399,10 @@ def _smart_eigh(C, info, rank, scalings=None, projs=None, # time saving short-circuit if proj_subspace and sum(rank.values()) == C.shape[0]: return np.ones(n_chan), np.eye(n_chan), np.ones(n_chan, bool) - eig = np.zeros(n_chan, dtype=C.dtype) - eigvec = np.zeros((n_chan, n_chan), dtype=C.dtype) + + dtype = np.complex if C.dtype == np.complex else np.float + eig = np.zeros(n_chan, dtype) + eigvec = np.zeros((n_chan, n_chan), dtype) mask = np.zeros(n_chan, bool) for ch_type, picks in _picks_by_type(info, meg_combined=True, ref_meg=False, exclude='bads'): @@ -1717,7 +1719,8 @@ def compute_whitener(noise_cov, info=None, picks=None, rank=None, nzero = (eig > 0) eig[~nzero] = 0. # get rid of numerical noise (negative) ones - W = np.zeros((n_chan, 1), dtype=noise_cov['eigvec'].dtype) + dtype = np.complex if noise_cov['eigvec'].dtype == np.complex else np.float + W = np.zeros((n_chan, 1), dtype) W[nzero, 0] = 1.0 / np.sqrt(eig[nzero]) # Rows of eigvec are the eigenvectors W = W * noise_cov['eigvec'] # C ** -0.5 From a132bb3b3854cde7ffef3bb6791db4e6e943e7e8 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Mon, 11 Nov 2019 12:00:10 +0200 Subject: [PATCH 13/16] PEP8 --- mne/beamformer/_dics.py | 12 ++++++------ mne/io/tests/test_pick.py | 2 +- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/mne/beamformer/_dics.py b/mne/beamformer/_dics.py index abf7c267f7e..f3a3e3b85ed 100644 --- a/mne/beamformer/_dics.py +++ b/mne/beamformer/_dics.py @@ -8,7 +8,7 @@ # License: BSD (3-clause) import numpy as np -from ..io.pick import pick_info, pick_types +from ..io.pick import pick_info from ..utils import (logger, verbose, warn, _check_one_ch_type, _check_channels_spatial_filter, _check_rank, _check_option, _check_info_inv) @@ -20,6 +20,7 @@ _compute_beamformer, _check_src_type, Beamformer, _compute_power) + @verbose def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, pick_ori=None, rank=None, inversion='single', weight_norm=None, @@ -235,7 +236,6 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, subject = _subject_from_forward(forward) ch_names = list(info['ch_names']) - #csd_picks = pick_types(info, meg=True, eeg=True, seeg=True) csd_picks = [csd.ch_names.index(ch) for ch in ch_names] logger.info('Computing DICS spatial filters...') @@ -256,10 +256,10 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, else: noise_csd_freq = None - _, _, proj, vertices, G, whitener, nn, orient_std = _prepare_beamformer_input( - info, forward, label, pick_ori, noise_csd_freq, - combine_xyz=combine_xyz, exp=exp) - + _, _, proj, vertices, G, whitener, nn, orient_std = \ + _prepare_beamformer_input( + info, forward, label, pick_ori, noise_csd_freq, + combine_xyz=combine_xyz, exp=exp) # Ensure the CSD is in the same order as the leadfield Cm = Cm[csd_picks, :][:, csd_picks] diff --git a/mne/io/tests/test_pick.py b/mne/io/tests/test_pick.py index e7d2e848d96..8ba2bdbeda9 100644 --- a/mne/io/tests/test_pick.py +++ b/mne/io/tests/test_pick.py @@ -6,7 +6,7 @@ import pytest import numpy as np -from mne import (pick_channels_regexp, pick_types, Epochs, EpochsArray, +from mne import (pick_channels_regexp, pick_types, Epochs, read_forward_solution, rename_channels, pick_info, pick_channels, create_info) from mne import __file__ as _root_init_fname From 60066c85af4214bd44a36f100016bbe33dd4dfc0 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Mon, 6 Jan 2020 10:24:44 +0200 Subject: [PATCH 14/16] Fix the whitening! --- mne/beamformer/_dics.py | 9 ++++++++- mne/beamformer/_lcmv.py | 8 ++++++-- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/mne/beamformer/_dics.py b/mne/beamformer/_dics.py index f3a3e3b85ed..424a3fb8fe2 100644 --- a/mne/beamformer/_dics.py +++ b/mne/beamformer/_dics.py @@ -240,6 +240,7 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, logger.info('Computing DICS spatial filters...') Ws = [] + whiteners = [] for i, freq in enumerate(frequencies): if n_freqs > 1: logger.info(' computing DICS spatial filter at %sHz (%d/%d)' % @@ -272,8 +273,10 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, reduce_rank, rank=rank, inversion=inversion, nn=nn, orient_std=orient_std) Ws.append(W) + whiteners.append(whitener) Ws = np.array(Ws) + whiteners = np.array(whiteners) src_type = _get_src_type(forward['src'], vertices) filters = Beamformer( @@ -281,7 +284,7 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, vertices=vertices, subject=subject, pick_ori=pick_ori, inversion=inversion, weight_norm=weight_norm, normalize_fwd=bool(normalize_fwd), src_type=src_type, - n_orient=n_orient if pick_ori is None else 1) + n_orient=n_orient if pick_ori is None else 1, whiteners=whiteners) return filters @@ -505,6 +508,10 @@ def apply_dics_csd(csd, filters, verbose=None): Cm = Cm[csd_picks, :][:, csd_picks] W = filters['weights'][i] + # Whiten the CSD + whitener = filters['whiteners'][i] + Cm = np.dot(whitener, np.dot(Cm, whitener.conj().T)) + source_power[:, i] = _compute_power(Cm, W, n_orient) logger.info('[done]') diff --git a/mne/beamformer/_lcmv.py b/mne/beamformer/_lcmv.py index 29fd42625a1..2ac75ab2342 100644 --- a/mne/beamformer/_lcmv.py +++ b/mne/beamformer/_lcmv.py @@ -429,8 +429,12 @@ def apply_lcmv_cov(data_cov, filters, verbose=None): data_cov = pick_channels_cov(data_cov, sel_names) n_orient = filters['weights'].shape[0] // filters['nsource'] - source_power = _compute_power(data_cov['data'], filters['weights'], - n_orient) + + # Whiten the CSD + whitener = filters['whitener'] + Cm = np.dot(whitener, np.dot(data_cov['data'], whitener.conj().T)) + + source_power = _compute_power(Cm, filters['weights'], n_orient) # compatibility with 0.16, add src_type as None if not present: filters, warn_text = _check_src_type(filters) From ce8bd08ade778bf4cedf737f9929f5808f629e39 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Wed, 8 Jan 2020 13:05:42 +0200 Subject: [PATCH 15/16] Only have one noise CSD --- mne/beamformer/_dics.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/mne/beamformer/_dics.py b/mne/beamformer/_dics.py index 424a3fb8fe2..d8396628079 100644 --- a/mne/beamformer/_dics.py +++ b/mne/beamformer/_dics.py @@ -238,6 +238,18 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, csd_picks = [csd.ch_names.index(ch) for ch in ch_names] + if noise_csd is not None: + if len(noise_csd.frequencies) > 1: + noise_csd = noise_csd.mean() + noise_csd = noise_csd.get_data(as_cov=True) + if real_filter: + noise_csd['data'] = noise_csd['data'].real + + _, _, proj, vertices, G, whitener, nn, orient_std = \ + _prepare_beamformer_input( + info, forward, label, pick_ori, noise_csd, + combine_xyz=combine_xyz, exp=exp) + logger.info('Computing DICS spatial filters...') Ws = [] whiteners = [] @@ -250,18 +262,6 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, if real_filter: Cm = Cm.real - if noise_csd is not None: - noise_csd_freq = noise_csd.get_data(index=i, as_cov=True) - if real_filter: - noise_csd_freq['data'] = noise_csd_freq['data'].real - else: - noise_csd_freq = None - - _, _, proj, vertices, G, whitener, nn, orient_std = \ - _prepare_beamformer_input( - info, forward, label, pick_ori, noise_csd_freq, - combine_xyz=combine_xyz, exp=exp) - # Ensure the CSD is in the same order as the leadfield Cm = Cm[csd_picks, :][:, csd_picks] From 1ffebafcfb3fb7ec096f3261a4170b9910efbc72 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Thu, 13 Feb 2020 12:14:13 +0000 Subject: [PATCH 16/16] Bit of cleanup --- mne/beamformer/_dics.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/mne/beamformer/_dics.py b/mne/beamformer/_dics.py index d8396628079..b45458a5053 100644 --- a/mne/beamformer/_dics.py +++ b/mne/beamformer/_dics.py @@ -239,6 +239,7 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, csd_picks = [csd.ch_names.index(ch) for ch in ch_names] if noise_csd is not None: + # Use the same noise CSD for all frequencies if len(noise_csd.frequencies) > 1: noise_csd = noise_csd.mean() noise_csd = noise_csd.get_data(as_cov=True) @@ -247,12 +248,11 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, _, _, proj, vertices, G, whitener, nn, orient_std = \ _prepare_beamformer_input( - info, forward, label, pick_ori, noise_csd, - combine_xyz=combine_xyz, exp=exp) + info, forward, label, pick_ori, noise_cov=noise_csd, rank=rank, + pca=False, combine_xyz=combine_xyz, exp=exp) logger.info('Computing DICS spatial filters...') Ws = [] - whiteners = [] for i, freq in enumerate(frequencies): if n_freqs > 1: logger.info(' computing DICS spatial filter at %sHz (%d/%d)' % @@ -273,10 +273,8 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, reduce_rank, rank=rank, inversion=inversion, nn=nn, orient_std=orient_std) Ws.append(W) - whiteners.append(whitener) Ws = np.array(Ws) - whiteners = np.array(whiteners) src_type = _get_src_type(forward['src'], vertices) filters = Beamformer( @@ -284,7 +282,7 @@ def make_dics(info, forward, csd, reg=0.05, noise_csd=None, label=None, vertices=vertices, subject=subject, pick_ori=pick_ori, inversion=inversion, weight_norm=weight_norm, normalize_fwd=bool(normalize_fwd), src_type=src_type, - n_orient=n_orient if pick_ori is None else 1, whiteners=whiteners) + n_orient=n_orient if pick_ori is None else 1, whitener=whitener) return filters @@ -487,6 +485,7 @@ def apply_dics_csd(csd, filters, verbose=None): vertices = filters['vertices'] n_orient = filters['n_orient'] subject = filters['subject'] + whitener = filters['whitener'] n_sources = np.sum([len(v) for v in vertices]) # If CSD is summed over multiple frequencies, take the average frequency @@ -509,7 +508,6 @@ def apply_dics_csd(csd, filters, verbose=None): W = filters['weights'][i] # Whiten the CSD - whitener = filters['whiteners'][i] Cm = np.dot(whitener, np.dot(Cm, whitener.conj().T)) source_power[:, i] = _compute_power(Cm, W, n_orient)