FIX: Fix simulate_evoked and apply_forward

doc/changes/latest.inc

@@ -99,6 +99,8 @@ Bugs
 
 - Fixed bug in :meth:`mne.Epochs.drop_bad` where subsequent rejections failed if they only specified thresholds for a subset of the channel types used in a previous rejection (:gh:`9485` by `Richard Höchenberger`_).
 
+- Fix bug with `mne.simulation.simulate_evoked`, `mne.apply_forward`, and `mne.apply_forward_raw` where systems with EEG channels that come before MEG channels would have them mixed up in the output evoked or raw object (:gh:`#9513` by `Eric Larson`_)
+
 - In :func:`mne.viz.plot_ica_scores` and :meth:`mne.preprocessing.ICA.plot_scores`, the figure and axis titles no longer overlap when plotting only a single EOG or ECG channel (:gh:`9489` by `Richard Höchenberger`_).
 
 API changes

mne/forward/forward.py

@@ -1298,7 +1298,7 @@ def _stc_src_sel(src, stc, on_missing='raise',
     return src_sel, stc_sel, out_vertices
 
 
-def _fill_measurement_info(info, fwd, sfreq):
+def _fill_measurement_info(info, fwd, sfreq, data):
     """Fill the measurement info of a Raw or Evoked object."""
     sel = pick_channels(info['ch_names'], fwd['sol']['row_names'])
     info = pick_info(info, sel)
@@ -1316,7 +1316,11 @@ def _fill_measurement_info(info, fwd, sfreq):
                 lowpass=sfreq / 2., sfreq=sfreq, projs=[])
     info._check_consistency()
 
-    return info
+    # reorder data (which is in fwd order) to match that of info
+    order = [fwd['sol']['row_names'].index(name) for name in info['ch_names']]
+    data = data[order]
+
+    return info, data
 
 
 @verbose
@@ -1418,9 +1422,9 @@ def apply_forward(fwd, stc, info, start=None, stop=None, use_cps=True,
 
     # fill the measurement info
     sfreq = float(1.0 / stc.tstep)
-    info_out = _fill_measurement_info(info, fwd, sfreq)
+    info, data = _fill_measurement_info(info, fwd, sfreq, data)
 
-    evoked = EvokedArray(data, info_out, times[0], nave=1)
+    evoked = EvokedArray(data, info, times[0], nave=1)
 
     evoked.times = times
     evoked._update_first_last()
@@ -1483,7 +1487,7 @@ def apply_forward_raw(fwd, stc, info, start=None, stop=None,
                                  use_cps=use_cps)
 
     sfreq = 1.0 / stc.tstep
-    info = _fill_measurement_info(info, fwd, sfreq)
+    info, data = _fill_measurement_info(info, fwd, sfreq, data)
     info['projs'] = []
     # store sensor data in Raw object using the info
     raw = RawArray(data, info)

mne/simulation/evoked.py

@@ -15,7 +15,7 @@
 
 
 @verbose
-def simulate_evoked(fwd, stc, info, cov, nave=30, iir_filter=None,
+def simulate_evoked(fwd, stc, info, cov=None, nave=30, iir_filter=None,
                     random_state=None, use_cps=True, verbose=None):
     """Generate noisy evoked data.
 
@@ -33,8 +33,8 @@ def simulate_evoked(fwd, stc, info, cov, nave=30, iir_filter=None,
         The source time courses.
     info : dict
         Measurement info to generate the evoked.
-    cov : Covariance object
-        The noise covariance.
+    cov : Covariance object | None
+        The noise covariance. If None, no noise is added.
     nave : int
         Number of averaged epochs (defaults to 30).
 
@@ -70,11 +70,14 @@ def simulate_evoked(fwd, stc, info, cov, nave=30, iir_filter=None,
     .. versionadded:: 0.10.0
     """
     evoked = apply_forward(fwd, stc, info, use_cps=use_cps)
+    if cov is None:
+        return evoked
+
     if nave < np.inf:
         noise = _simulate_noise_evoked(evoked, cov, iir_filter, random_state)
         evoked.data += noise.data / math.sqrt(nave)
         evoked.nave = np.int64(nave)
-    if cov is not None and cov.get('projs', None):
+    if cov.get('projs', None):
         evoked.add_proj(cov['projs']).apply_proj()
     return evoked
 

mne/simulation/tests/test_evoked.py

@@ -158,3 +158,25 @@ def test_rank_deficiency():
     assert cov['names'] == cov_new['names']
     r = np.corrcoef(cov['data'].ravel(), cov_new['data'].ravel())[0, 1]
     assert r > 0.98
+
+
+@testing.requires_testing_data
+def test_order():
+    """Test that order does not matter."""
+    fwd = read_forward_solution(fwd_fname)
+    fwd = convert_forward_solution(fwd, force_fixed=True, use_cps=False)
+    evoked = read_evokeds(ave_fname)[0].pick_types(meg=True, eeg=True)
+    assert 'meg' in evoked
+    assert 'eeg' in evoked
+    meg_picks = pick_types(evoked.info, meg=True)
+    eeg_picks = pick_types(evoked.info, eeg=True)
+    # MEG then EEG
+    assert (eeg_picks > meg_picks.max()).all()
+    times = np.arange(10) / 1000.
+    stc = simulate_sparse_stc(fwd['src'], 1, times=times, random_state=0)
+    evoked_sim = simulate_evoked(fwd, stc, evoked.info, nave=np.inf)
+    reorder = np.concatenate([eeg_picks, meg_picks])
+    evoked.reorder_channels([evoked.ch_names[pick] for pick in reorder])
+    evoked_sim_2 = simulate_evoked(fwd, stc, evoked.info, nave=np.inf)
+    want_data = evoked_sim.data[reorder]
+    assert_allclose(evoked_sim_2.data, want_data)