FIX: prevent raw._data overwrite from _add_exg and _add_chpi (mne-too…

…ls#8633)

* FIX: prevent raw._data overwrite from _add_exg and _add_chpi

* Test function for simulation bug fix

* Fixed pep8 issues

* Import ordering changed according to coding guidelines

* Cascading simulation test moved to simulation dir

* Cleanup

* Pull update for latest.inc

* Update - contrib doc

* WIP - failing test for chpi simulation

* FIX: Faster test

* FIX: latest

* FIX: Missed one

* FIX: Missed

Co-authored-by: Jeff Stout <[email protected]>
Co-authored-by: Eric Larson <[email protected]>

doc/changes/latest.inc

@@ -30,6 +30,8 @@ Current (0.22.dev0)
 
 .. |Qianliang Li| replace:: **Qianliang Li**
 
+.. |Jeff Stout| replace:: **Jeff Stout**
+
 
 Enhancements
 ~~~~~~~~~~~~
@@ -103,6 +105,8 @@ Bugs
 ~~~~
 - Fix orthogonalization of power envelopes in :func:`mne.connectivity.envelope_correlation` (:gh:`8658` **by new contributor** |Qianliang Li|_ and `Eric Larson`_)
 
+- Fix data overwrite of cascading simulation operations :`mne.simulation.simulate_raw` (:gh:`8633` **by new contributor** |Jeff Stout|_)
+
 - Fix a transpose issue of :func:`mne.decoding.CSP.plot_filters` (:gh:`8580` **by new contributor** |Hongjiang Ye|_)
 
 - Fix :func:`mne.io.read_raw_curry` to deal with Curry datasets that have channels that are listed in the labels file, but which are absent from the saved data file (e.g. 'Ref' channel).  Also now populates info['meas_date'] if possible (:gh:`8400` **by new contributor** |Tod Flak|_)

doc/changes/names.inc

@@ -342,4 +342,6 @@
 
 .. _Hongjiang Ye: https://github.com/rubyyhj
 
+.. _Jeff Stout: https://megcore.nih.gov/index.php/Staff
+
 .. _Qianliang Li: https://www.dtu.dk/english/service/phonebook/person?id=126774

mne/simulation/raw.py

@@ -524,7 +524,7 @@ def _add_exg(raw, kind, head_pos, interp, n_jobs, random_state):
     proc_lims = np.concatenate([np.arange(0, len(used), 10000), [len(used)]])
     for start, stop in zip(proc_lims[:-1], proc_lims[1:]):
         fwd, _ = interper.feed(stop - start)
-        data[picks, start:stop] = einsum(
+        data[picks, start:stop] += einsum(
             'svt,vt->st', fwd, exg_data[:, start:stop])
         assert not used[start:stop].any()
         used[start:stop] = True
@@ -580,7 +580,7 @@ def add_chpi(raw, head_pos=None, interp='cos2', n_jobs=1, verbose=None):
     lims = np.concatenate([offsets, [len(raw.times)]])
     for start, stop in zip(lims[:-1], lims[1:]):
         fwd, = interper.feed(stop - start)
-        data[meg_picks, start:stop] = einsum(
+        data[meg_picks, start:stop] += einsum(
             'svt,vt->st', fwd, sinusoids[:, start:stop])
         assert not used[start:stop].any()
         used[start:stop] = True

mne/simulation/tests/test_raw.py

@@ -24,13 +24,13 @@
 from mne.tests.test_chpi import _assert_quats
 from mne.datasets import testing
 from mne.simulation import (simulate_sparse_stc, simulate_raw, add_eog,
-                            add_ecg, add_chpi)
+                            add_ecg, add_chpi, add_noise)
 from mne.source_space import _compare_source_spaces
 from mne.surface import _get_ico_surface
 from mne.io import read_raw_fif, RawArray
 from mne.io.constants import FIFF
 from mne.time_frequency import psd_welch
-from mne.utils import run_tests_if_main, catch_logging, check_version
+from mne.utils import catch_logging, check_version
 
 base_path = op.join(op.dirname(__file__), '..', '..', 'io', 'tests', 'data')
 raw_fname_short = op.join(base_path, 'test_raw.fif')
@@ -479,4 +479,39 @@ def test_simulate_raw_chpi():
                   vel_atol=0.03)  # velicity huge because of t_step_min above
 
 
-run_tests_if_main()
+@testing.requires_testing_data
+def test_simulation_cascade():
+    """Test that cascading operations do not overwrite data."""
+    # Create 10 second raw dataset with zeros in the data matrix
+    raw_null = read_raw_fif(raw_chpi_fname, allow_maxshield='yes')
+    raw_null.crop(0, 1).pick_types(meg=True).load_data()
+    raw_null.apply_function(lambda x: np.zeros_like(x))
+    assert_array_equal(raw_null.get_data(), 0.)
+
+    # Calculate independent signal additions
+    raw_eog = raw_null.copy()
+    add_eog(raw_eog, random_state=0)
+
+    raw_ecg = raw_null.copy()
+    add_ecg(raw_ecg, random_state=0)
+
+    raw_noise = raw_null.copy()
+    cov = make_ad_hoc_cov(raw_null.info)
+    add_noise(raw_noise, cov, random_state=0)
+
+    raw_chpi = raw_null.copy()
+    add_chpi(raw_chpi)
+
+    # Calculate Cascading signal additions
+    raw_cascade = raw_null.copy()
+    add_eog(raw_cascade, random_state=0)
+    add_ecg(raw_cascade, random_state=0)
+    add_chpi(raw_cascade)
+    add_noise(raw_cascade, cov, random_state=0)
+
+    cascade_data = raw_cascade.get_data()
+    serial_data = 0.
+    for raw_other in (raw_eog, raw_ecg, raw_noise, raw_chpi):
+        serial_data += raw_other.get_data()
+
+    assert_allclose(cascade_data, serial_data, atol=1e-20)