BUG: Fix bug with simulating head pos and BEM (#13276)

Co-authored-by: autofix-ci[bot] <114827586+autofix-ci[bot]@users.noreply.github.com>

Author: larsoner

doc/changes/devel/13276.bugfix.rst

@@ -0,0 +1 @@
+Fix bug with :func:`mne.simulation.simulate_raw` and :class:`mne.simulation.SourceSimulator` where using different head positions with ``head_pos`` and a BEM would raise an error, by `Eric Larson`_.

mne/chpi.py

@@ -99,8 +99,9 @@ def read_head_pos(fname):
 
     Returns
     -------
-    pos : array, shape (N, 10)
+    quats : array, shape (n_pos, 10)
         The position and quaternion parameters from cHPI fitting.
+        See :func:`mne.chpi.compute_head_pos` for details on the columns.
 
     See Also
     --------
@@ -126,8 +127,9 @@ def write_head_pos(fname, pos):
     ----------
     fname : path-like
         The filename to write.
-    pos : array, shape (N, 10)
+    pos : array, shape (n_pos, 10)
         The position and quaternion parameters from cHPI fitting.
+        See :func:`mne.chpi.compute_head_pos` for details on the columns.
 
     See Also
     --------
@@ -141,7 +143,9 @@ def write_head_pos(fname, pos):
     _check_fname(fname, overwrite=True)
     pos = np.array(pos, np.float64)
     if pos.ndim != 2 or pos.shape[1] != 10:
-        raise ValueError("pos must be a 2D array of shape (N, 10)")
+        raise ValueError(
+            f"pos must be a 2D array of shape (N, 10), got shape {pos.shape}"
+        )
     with open(fname, "wb") as fid:
         fid.write(
             " Time       q1       q2       q3       q4       q5       "
@@ -157,16 +161,17 @@ def head_pos_to_trans_rot_t(quats):
 
     Parameters
     ----------
-    quats : ndarray, shape (N, 10)
+    quats : ndarray, shape (n_pos, 10)
         MaxFilter-formatted position and quaternion parameters.
+        See :func:`mne.chpi.read_head_pos` for details on the columns.
 
     Returns
     -------
-    translation : ndarray, shape (N, 3)
+    translation : ndarray, shape (n_pos, 3)
         Translations at each time point.
-    rotation : ndarray, shape (N, 3, 3)
+    rotation : ndarray, shape (n_pos, 3, 3)
         Rotations at each time point.
-    t : ndarray, shape (N,)
+    t : ndarray, shape (n_pos,)
         The time points.
 
     See Also
@@ -929,7 +934,8 @@ def compute_head_pos(
     Returns
     -------
     quats : ndarray, shape (n_pos, 10)
-        The ``[t, q1, q2, q3, x, y, z, gof, err, v]`` for each fit.
+        MaxFilter-formatted head position parameters. The columns correspond to
+        ``[t, q1, q2, q3, x, y, z, gof, err, v]`` for each time point.
 
     See Also
     --------

mne/forward/_compute_forward.py

@@ -35,6 +35,7 @@ def _dup_coil_set(coils, coord_frame, t):
     if t is not None:
         coord_frame = t["to"]
         for coil in coils:
+            assert isinstance(coil, dict), f"Coil must be a dict, got {type(coil)}"
             for key in ("ex", "ey", "ez"):
                 if key in coil:
                     coil[key] = apply_trans(t["trans"], coil[key], False)
@@ -794,6 +795,7 @@ def _compute_forwards_meeg(rr, *, sensors, fwd_data, n_jobs, silent=False):
     mri_Q, bem_rr, fun = fwd_data["mri_Q"], fwd_data["bem_rr"], fwd_data["fun"]
     solutions = fwd_data["solutions"]
     del fwd_data
+    rr = np.ascontiguousarray(rr)  # usually true but not guaranteed, e.g. in dipole.py
     for coil_type, sens in sensors.items():
         coils = sens["defs"]
         compensator = sens.get("compensator", None)
@@ -835,6 +837,9 @@ def _compute_forwards(rr, *, bem, sensors, n_jobs, verbose=None):
     solver = bem.get("solver", "mne")
     _check_option("solver", solver, ("mne", "openmeeg"))
     if bem["is_sphere"] or solver == "mne":
+        # This modifies "sensors" in place, so let's copy it in case the calling
+        # function needs to reuse it (e.g., in simulate_raw.py)
+        sensors = deepcopy(sensors)
         fwd_data = _prep_field_computation(rr, sensors=sensors, bem=bem, n_jobs=n_jobs)
         Bs = _compute_forwards_meeg(
             rr, sensors=sensors, fwd_data=fwd_data, n_jobs=n_jobs

mne/simulation/tests/test_raw.py

@@ -38,6 +38,7 @@
     compute_head_pos,
     get_chpi_info,
     read_head_pos,
+    write_head_pos,
 )
 from mne.datasets import testing
 from mne.io import RawArray, read_raw_fif
@@ -54,6 +55,7 @@
 from mne.source_space._source_space import _compare_source_spaces
 from mne.surface import _get_ico_surface
 from mne.tests.test_chpi import _assert_quats
+from mne.transforms import _affine_to_quat
 from mne.utils import catch_logging
 
 raw_fname_short = Path(__file__).parents[2] / "io" / "tests" / "data" / "test_raw.fif"
@@ -251,6 +253,11 @@ def test_simulate_raw_sphere(raw_data, tmp_path):
 
     # head pos
     head_pos_sim = _get_head_pos_sim(raw)
+    head_pos_sim_2 = np.zeros((len(head_pos_sim), 10))
+    for ii, (t, mat) in enumerate(head_pos_sim.items()):
+        head_pos_sim_2[ii, :7] = [t] + list(_affine_to_quat(mat))
+    head_pos_sim_3 = tmp_path / "head_pos.txt"
+    write_head_pos(head_pos_sim_3, head_pos_sim_2)
 
     #
     # Test raw simulation with basic parameters
@@ -259,11 +266,9 @@ def test_simulate_raw_sphere(raw_data, tmp_path):
     cov = read_cov(cov_fname)
     cov["projs"] = raw.info["projs"]
     raw.info["bads"] = raw.ch_names[:1]
-    sphere_norad = make_sphere_model("auto", None, raw.info)
     raw_meg = raw.copy().pick("meg")
-    raw_sim = simulate_raw(
-        raw_meg.info, stc, trans, src, sphere_norad, head_pos=head_pos_sim
-    )
+    raw_sim = simulate_raw(raw_meg.info, stc, trans, src, sphere, head_pos=head_pos_sim)
+    raw_data = raw_sim[:][0]
     # Test IO on processed data
     test_outname = tmp_path / "sim_test_raw.fif"
     raw_sim.save(test_outname)
@@ -307,12 +312,14 @@ def test_simulate_raw_sphere(raw_data, tmp_path):
     )
     del raw_sim, raw_sim_2
 
-    # check that different interpolations are similar given small movements
+    # check that different interpolations are similar given small movements,
+    # using different input forms of head_pos
     raw_sim = simulate_raw(
         raw.info, stc, trans, src, sphere, head_pos=head_pos_sim, interp="linear"
     )
+    assert_allclose(raw_sim.get_data("meg"), raw_data, rtol=0.02)
     raw_sim_hann = simulate_raw(
-        raw.info, stc, trans, src, sphere, head_pos=head_pos_sim, interp="hann"
+        raw.info, stc, trans, src, sphere, head_pos=head_pos_sim_3, interp="hann"
     )
     assert_allclose(raw_sim[:][0], raw_sim_hann[:][0], rtol=1e-1, atol=1e-14)
     del raw_sim_hann
@@ -391,7 +398,7 @@ def test_simulate_raw_bem(raw_data):
     cov = make_ad_hoc_cov(raw.info)
     # The tolerance for the BEM is surprisingly high but I get the same
     # result when using MNE-C and Xfit, even when using a proper 5120 BEM :(
-    for use_raw, bem, tol in ((raw_sim_sph, sphere, 4), (raw_sim_bem, bem_fname, 31)):
+    for use_raw, bem, tol in ((raw_sim_sph, sphere, 6), (raw_sim_bem, bem_fname, 31)):
         events = find_events(use_raw, "STI 014")
         assert len(locs) == 6
         evoked = Epochs(use_raw, events, 1, 0, tmax, baseline=None).average()
@@ -425,6 +432,18 @@ def test_simulate_raw_bem(raw_data):
     assert_allclose(amp0 / amp1, wf_sim[0] / wf_sim[1], rtol=1e-5)
     assert amp2 == 0
     assert raw_sim.n_times == ss.n_times
+    # smoke test that different head positions can be used as well
+    head_pos_sim = {1.0: raw.info["dev_head_t"]["trans"]}
+    raw_sim_2 = simulate_raw(
+        raw.info,
+        ss,
+        src=src_ss,
+        bem=bem_fname,
+        first_samp=first_samp,
+        head_pos=head_pos_sim,
+    )
+    data_2 = raw_sim_2.get_data()
+    assert_allclose(data, data_2, rtol=1e-7)
 
 
 @pytest.mark.slowtest  # slow on Windows Azure