From f3c1577981faa5746a73f38c76c7b2120a2879ad Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 23 Jan 2019 17:28:18 -0500 Subject: [PATCH 1/4] ENH: Add add_noise to simulation --- doc/python_reference.rst | 1 + doc/whats_new.rst | 2 + mne/simulation/__init__.py | 2 +- mne/simulation/evoked.py | 66 ++++++++++++++++++++++++++--- mne/simulation/tests/test_evoked.py | 39 +++++++++++++++-- 5 files changed, 101 insertions(+), 9 deletions(-) diff --git a/doc/python_reference.rst b/doc/python_reference.rst index c08c7f5f140..307a7238f20 100644 --- a/doc/python_reference.rst +++ b/doc/python_reference.rst @@ -876,6 +876,7 @@ Simulation .. autosummary:: :toctree: generated/ + add_noise simulate_evoked simulate_raw simulate_stc diff --git a/doc/whats_new.rst b/doc/whats_new.rst index 15e60f20944..0e97245da0f 100644 --- a/doc/whats_new.rst +++ b/doc/whats_new.rst @@ -31,6 +31,8 @@ Changelog - :func:`mne.io.read_raw_edf` now detects analog stim channels labeled ``'STATUS'`` and sets them as stim channel. :func:`mne.io.read_raw_edf` no longer synthesize TAL annotations into stim channel but stores them in ``raw.annotations`` when reading by `Joan Massich`_ +- Add `mne.simulation.add_noise` for ad-hoc noise addition to `io.Raw`, `Epochs`, and `Evoked` instances, by `Eric Larson`_ + - Add ``drop_refs=True`` parameter to :func:`set_bipolar_reference` to drop/keep anode and cathode channels after applying the reference by `Cristóbal Moënne-Loccoz`_. - Add use of :func:`scipy.signal.windows.dpss` for faster multitaper window computations in PSD functions by `Eric Larson`_ diff --git a/mne/simulation/__init__.py b/mne/simulation/__init__.py index 2958f5041c2..4fdd8e809c8 100644 --- a/mne/simulation/__init__.py +++ b/mne/simulation/__init__.py @@ -1,6 +1,6 @@ """Data simulation code.""" -from .evoked import simulate_evoked, simulate_noise_evoked +from .evoked import simulate_evoked, simulate_noise_evoked, add_noise from .raw import simulate_raw from .source import select_source_in_label, simulate_stc, simulate_sparse_stc from .metrics import source_estimate_quantification diff --git a/mne/simulation/evoked.py b/mne/simulation/evoked.py index cf593c0d310..e181118e30b 100644 --- a/mne/simulation/evoked.py +++ b/mne/simulation/evoked.py @@ -8,9 +8,9 @@ import numpy as np -from ..io.pick import pick_channels_cov +from ..io.pick import pick_channels_cov, pick_info from ..forward import apply_forward -from ..utils import check_random_state, verbose +from ..utils import logger, verbose, check_random_state, _check_preload @verbose @@ -107,9 +107,65 @@ def simulate_noise_evoked(evoked, cov, iir_filter=None, random_state=None): .. versionadded:: 0.10.0 """ noise = evoked.copy() - noise.data = _generate_noise(evoked.info, cov, iir_filter, random_state, - evoked.data.shape[1])[0] - return noise + noise.data[:] = 0 + return _add_noise(noise, cov, iir_filter, random_state) + + +@verbose +def add_noise(inst, cov, iir_filter=None, random_state=None, + verbose=None): + """Create noise as a multivariate Gaussian. + + The spatial covariance of the noise is given from the cov matrix. + + Parameters + ---------- + inst : instance of Evoked, Epochs, or Raw + Instance to which to add noise. + cov : instance of Covariance + The noise covariance. + iir_filter : None | array-like + IIR filter coefficients (denominator). + random_state : None | int | np.random.RandomState + To specify the random generator state. + + Returns + ------- + noise : evoked object + an instance of evoked + + Notes + ----- + Only channels in both ``inst.info['ch_names']`` and + ``cov['names']`` will have noise added to them. + + .. versionadded:: 0.18.0 + """ + # We always allow subselection here + return _add_noise( + inst, cov, iir_filter, random_state, allow_subselection=True) + + +def _add_noise(inst, cov, iir_filter, random_state, allow_subselection=False): + """Add noise, possibly with channel subselection.""" + _check_preload(inst, 'Adding noise') + data = inst._data + assert data.ndim in (2, 3) + if data.ndim == 2: + data = data[np.newaxis] + # Subselect if necessary + info = inst.info + picks = slice(None) + if allow_subselection: + use_chs = list(set(info['ch_names']) & set(cov['names'])) + picks = np.where(np.in1d(info['ch_names'], use_chs))[0] + logger.info('Adding noise to %d/%d channels (%d channels in cov)' + % (len(picks), len(info['chs']), len(cov['names']))) + info = pick_info(inst.info, picks) + for epoch in data: + epoch[picks] = _generate_noise(info, cov, iir_filter, random_state, + epoch.shape[1])[0] + return inst def _generate_noise(info, cov, iir_filter, random_state, n_samples, zi=None): diff --git a/mne/simulation/tests/test_evoked.py b/mne/simulation/tests/test_evoked.py index 3538ba56e9d..3c94042c373 100644 --- a/mne/simulation/tests/test_evoked.py +++ b/mne/simulation/tests/test_evoked.py @@ -10,12 +10,13 @@ import pytest from mne import (read_cov, read_forward_solution, convert_forward_solution, - pick_types_forward, read_evokeds) + pick_types_forward, read_evokeds, pick_types, EpochsArray, + compute_covariance, compute_raw_covariance) from mne.datasets import testing -from mne.simulation import simulate_sparse_stc, simulate_evoked +from mne.simulation import simulate_sparse_stc, simulate_evoked, add_noise from mne.io import read_raw_fif from mne.cov import regularize -from mne.utils import run_tests_if_main +from mne.utils import run_tests_if_main, catch_logging data_path = testing.data_path(download=False) fwd_fname = op.join(data_path, 'MEG', 'sample', @@ -79,4 +80,36 @@ def test_simulate_evoked(): cov, nave=nave, iir_filter=None) +@testing.requires_testing_data +def test_add_noise(): + """Test noise addition.""" + data_path = testing.data_path() + raw = read_raw_fif(data_path + '/MEG/sample/sample_audvis_trunc_raw.fif') + raw.del_proj() + picks = pick_types(raw.info, eeg=True, exclude=()) + cov = compute_raw_covariance(raw, picks=picks) + with pytest.raises(RuntimeError, match='to be loaded'): + add_noise(raw, cov) + raw.crop(0, 1).load_data() + raw._data[:] = 0. + epochs = EpochsArray(np.zeros((1, len(raw.ch_names), 100)), + raw.info.copy()) + epochs.info['bads'] = [] + for inst in (raw, epochs): + with catch_logging() as log: + add_noise(inst, cov, verbose=True) + log = log.getvalue() + want = ('to {0}/{1} channels ({0}' + .format(len(cov['names']), len(raw.ch_names))) + assert want in log + if inst is raw: + cov_new = compute_raw_covariance(raw, picks=picks, + verbose='error') # samples + else: + cov_new = compute_covariance(epochs, verbose='error') # avg ref + assert cov['names'] == cov_new['names'] + r = np.corrcoef(cov['data'].ravel(), cov_new['data'].ravel())[0, 1] + assert r > 0.99 + + run_tests_if_main() From 02d26dbc0390ba6ed58318359a48f3a4343edba1 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Thu, 24 Jan 2019 16:05:38 -0500 Subject: [PATCH 2/4] FIX: Minor fixes --- mne/simulation/evoked.py | 28 ++++++++++++++++++-------- mne/tests/test_docstring_parameters.py | 1 - 2 files changed, 20 insertions(+), 9 deletions(-) diff --git a/mne/simulation/evoked.py b/mne/simulation/evoked.py index e181118e30b..d96c77369b2 100644 --- a/mne/simulation/evoked.py +++ b/mne/simulation/evoked.py @@ -10,7 +10,8 @@ from ..io.pick import pick_channels_cov, pick_info from ..forward import apply_forward -from ..utils import logger, verbose, check_random_state, _check_preload +from ..utils import (logger, verbose, check_random_state, _check_preload, + deprecated) @verbose @@ -75,12 +76,14 @@ def simulate_evoked(fwd, stc, info, cov, nave=30, iir_filter=None, """ evoked = apply_forward(fwd, stc, info, use_cps=use_cps) if nave < np.inf: - noise = simulate_noise_evoked(evoked, cov, iir_filter, random_state) + noise = _simulate_noise_evoked(evoked, cov, iir_filter, random_state) evoked.data += noise.data / math.sqrt(nave) evoked.nave = np.int(nave) return evoked +@deprecated('simulate_noise_evoked is deprecated in 0.18 and will be removed ' + 'in 0.19, use add_noise instead') def simulate_noise_evoked(evoked, cov, iir_filter=None, random_state=None): """Create noise as a multivariate Gaussian. @@ -106,9 +109,14 @@ def simulate_noise_evoked(evoked, cov, iir_filter=None, random_state=None): ----- .. versionadded:: 0.10.0 """ + return _simulate_noise_evoked(evoked, cov, iir_filter, random_state) + + +def _simulate_noise_evoked(evoked, cov, iir_filter, random_state): noise = evoked.copy() noise.data[:] = 0 - return _add_noise(noise, cov, iir_filter, random_state) + return _add_noise(noise, cov, iir_filter, random_state, + allow_subselection=False) @verbose @@ -128,25 +136,29 @@ def add_noise(inst, cov, iir_filter=None, random_state=None, IIR filter coefficients (denominator). random_state : None | int | np.random.RandomState To specify the random generator state. + verbose : bool, str, int, or None + If not None, override default verbose level (see :func:`mne.verbose` + and :ref:`Logging documentation ` for more). Returns ------- - noise : evoked object - an instance of evoked + inst : instance of Evoked, Epochs, or Raw + The instance, modified to have additional noise. Notes ----- Only channels in both ``inst.info['ch_names']`` and ``cov['names']`` will have noise added to them. + This function operates inplace on ``inst``. + .. versionadded:: 0.18.0 """ # We always allow subselection here - return _add_noise( - inst, cov, iir_filter, random_state, allow_subselection=True) + return _add_noise(inst, cov, iir_filter, random_state) -def _add_noise(inst, cov, iir_filter, random_state, allow_subselection=False): +def _add_noise(inst, cov, iir_filter, random_state, allow_subselection=True): """Add noise, possibly with channel subselection.""" _check_preload(inst, 'Adding noise') data = inst._data diff --git a/mne/tests/test_docstring_parameters.py b/mne/tests/test_docstring_parameters.py index 44e1924d5c2..147d5131ce9 100644 --- a/mne/tests/test_docstring_parameters.py +++ b/mne/tests/test_docstring_parameters.py @@ -242,7 +242,6 @@ def test_tabs(): read_tag requires_sample_data rescale -simulate_noise_evoked source_estimate_quantification whiten_evoked write_fiducials From c94a153a794dfe1a4d4e7d16dc3347b60946af99 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Tue, 29 Jan 2019 12:26:04 -0500 Subject: [PATCH 3/4] FIX: Actually add values --- mne/simulation/evoked.py | 13 ++++++++++--- mne/simulation/tests/test_evoked.py | 23 +++++++++++++++++++---- 2 files changed, 29 insertions(+), 7 deletions(-) diff --git a/mne/simulation/evoked.py b/mne/simulation/evoked.py index d96c77369b2..346dea4cab2 100644 --- a/mne/simulation/evoked.py +++ b/mne/simulation/evoked.py @@ -11,7 +11,7 @@ from ..io.pick import pick_channels_cov, pick_info from ..forward import apply_forward from ..utils import (logger, verbose, check_random_state, _check_preload, - deprecated) + deprecated, _validate_type) @verbose @@ -160,6 +160,13 @@ def add_noise(inst, cov, iir_filter=None, random_state=None, def _add_noise(inst, cov, iir_filter, random_state, allow_subselection=True): """Add noise, possibly with channel subselection.""" + from ..cov import Covariance + from ..io import BaseRaw + from ..epochs import BaseEpochs + from ..evoked import Evoked + _validate_type(cov, Covariance, 'cov') + _validate_type(inst, (BaseRaw, BaseEpochs, Evoked), + 'inst', 'Raw, Epochs, or Evoked') _check_preload(inst, 'Adding noise') data = inst._data assert data.ndim in (2, 3) @@ -175,8 +182,8 @@ def _add_noise(inst, cov, iir_filter, random_state, allow_subselection=True): % (len(picks), len(info['chs']), len(cov['names']))) info = pick_info(inst.info, picks) for epoch in data: - epoch[picks] = _generate_noise(info, cov, iir_filter, random_state, - epoch.shape[1])[0] + epoch[picks] += _generate_noise(info, cov, iir_filter, random_state, + epoch.shape[1])[0] return inst diff --git a/mne/simulation/tests/test_evoked.py b/mne/simulation/tests/test_evoked.py index 3c94042c373..d6396151ef7 100644 --- a/mne/simulation/tests/test_evoked.py +++ b/mne/simulation/tests/test_evoked.py @@ -6,7 +6,7 @@ import numpy as np from numpy.testing import (assert_array_almost_equal, assert_array_equal, - assert_equal) + assert_equal, assert_allclose) import pytest from mne import (read_cov, read_forward_solution, convert_forward_solution, @@ -91,22 +91,37 @@ def test_add_noise(): with pytest.raises(RuntimeError, match='to be loaded'): add_noise(raw, cov) raw.crop(0, 1).load_data() + with pytest.raises(TypeError, match='Raw, Epochs, or Evoked'): + add_noise(0., cov) + with pytest.raises(TypeError, match='Covariance'): + add_noise(raw, 0.) + # test a no-op (data preserved) + orig_data = raw[:][0] + zero_cov = cov.copy() + zero_cov['data'].fill(0) + add_noise(raw, zero_cov) + new_data = raw[:][0] + assert_allclose(orig_data, new_data, atol=1e-30) + # set to zero to make comparisons easier raw._data[:] = 0. epochs = EpochsArray(np.zeros((1, len(raw.ch_names), 100)), raw.info.copy()) epochs.info['bads'] = [] - for inst in (raw, epochs): + evoked = epochs.average(picks=np.arange(len(raw.ch_names))) + for inst in (raw, epochs, evoked): with catch_logging() as log: add_noise(inst, cov, verbose=True) log = log.getvalue() want = ('to {0}/{1} channels ({0}' .format(len(cov['names']), len(raw.ch_names))) assert want in log + if inst is evoked: + inst = EpochsArray(inst.data[np.newaxis], inst.info) if inst is raw: - cov_new = compute_raw_covariance(raw, picks=picks, + cov_new = compute_raw_covariance(inst, picks=picks, verbose='error') # samples else: - cov_new = compute_covariance(epochs, verbose='error') # avg ref + cov_new = compute_covariance(inst, verbose='error') # avg ref assert cov['names'] == cov_new['names'] r = np.corrcoef(cov['data'].ravel(), cov_new['data'].ravel())[0, 1] assert r > 0.99 From 178e5b3afbca5855ba2aaf4ee76d8be3a69cf2ec Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Tue, 29 Jan 2019 13:13:05 -0500 Subject: [PATCH 4/4] FIX: Heisenbug --- mne/simulation/tests/test_evoked.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/mne/simulation/tests/test_evoked.py b/mne/simulation/tests/test_evoked.py index d6396151ef7..124d3395cec 100644 --- a/mne/simulation/tests/test_evoked.py +++ b/mne/simulation/tests/test_evoked.py @@ -83,6 +83,7 @@ def test_simulate_evoked(): @testing.requires_testing_data def test_add_noise(): """Test noise addition.""" + rng = np.random.RandomState(0) data_path = testing.data_path() raw = read_raw_fif(data_path + '/MEG/sample/sample_audvis_trunc_raw.fif') raw.del_proj() @@ -110,7 +111,7 @@ def test_add_noise(): evoked = epochs.average(picks=np.arange(len(raw.ch_names))) for inst in (raw, epochs, evoked): with catch_logging() as log: - add_noise(inst, cov, verbose=True) + add_noise(inst, cov, random_state=rng, verbose=True) log = log.getvalue() want = ('to {0}/{1} channels ({0}' .format(len(cov['names']), len(raw.ch_names)))