ENH: Improve OPM auditory dataset and example (#12539)

mne/datasets/config.py

@@ -92,6 +92,7 @@
     testing="0.152",
     misc="0.27",
     phantom_kit="0.2",
+    ucl_opm_auditory="0.2",
 )
 TESTING_VERSIONED = f'mne-testing-data-{RELEASES["testing"]}'
 MISC_VERSIONED = f'mne-misc-data-{RELEASES["misc"]}'
@@ -149,8 +150,8 @@
 
 MNE_DATASETS["ucl_opm_auditory"] = dict(
     archive_name="auditory_OPM_stationary.zip",
-    hash="md5:9ed0d8d554894542b56f8e7c4c0041fe",
-    url="https://osf.io/download/mwrt3/?version=1",
+    hash="md5:b2d69aa2d656b960bd0c18968dc1a14d",
+    url="https://osf.io/download/tp324/?version=1",  # original is mwrt3
     folder_name="auditory_OPM_stationary",
     config_key="MNE_DATASETS_UCL_OPM_AUDITORY_PATH",
 )

mne/preprocessing/tests/test_ica.py

@@ -1481,7 +1481,7 @@ def test_ica_labels():
 
     # derive reference ICA components and append them to raw
     ica_rf = ICA(n_components=2, max_iter=2, allow_ref_meg=True)
-    with _record_warnings(), pytest.warns(UserWarning, match="did not converge"):
+    with _record_warnings():  # high pass and/or no convergence
         ica_rf.fit(raw.copy().pick("ref_meg"))
     icacomps = ica_rf.get_sources(raw)
     # rename components so they are auto-detected by find_bads_ref

mne/utils/tests/test_numerics.py

@@ -450,7 +450,7 @@ def test_pca(n_components, whiten):
     assert_array_equal(X, X_orig)
     X_mne = pca_mne.fit_transform(X)
     assert_array_equal(X, X_orig)
-    assert_allclose(X_skl, X_mne)
+    assert_allclose(X_skl, X_mne * np.sign(np.sum(X_skl * X_mne, axis=0)))
     assert pca_mne.n_components_ == pca_skl.n_components_
     for key in (
         "mean_",
@@ -459,6 +459,10 @@ def test_pca(n_components, whiten):
         "explained_variance_ratio_",
     ):
         val_skl, val_mne = getattr(pca_skl, key), getattr(pca_mne, key)
+        if key == "components_":
+            val_mne = val_mne * np.sign(
+                np.sum(val_skl * val_mne, axis=1, keepdims=True)
+            )
         assert_allclose(val_skl, val_mne)
     if isinstance(n_components, float):
         assert pca_mne.n_components_ == n_dim - 1

tutorials/preprocessing/80_opm_processing.py

@@ -26,20 +26,20 @@
 # %%
 
 import matplotlib.pyplot as plt
+import nibabel as nib
 import numpy as np
 
 import mne
 
-opm_data_folder = mne.datasets.ucl_opm_auditory.data_path()
+subject = "sub-002"
+data_path = mne.datasets.ucl_opm_auditory.data_path()
 opm_file = (
-    opm_data_folder
-    / "sub-002"
-    / "ses-001"
-    / "meg"
-    / "sub-002_ses-001_task-aef_run-001_meg.bin"
+    data_path / subject / "ses-001" / "meg" / "sub-002_ses-001_task-aef_run-001_meg.bin"
 )
+subjects_dir = data_path / "derivatives" / "freesurfer" / "subjects"
+
 # For now we are going to assume the device and head coordinate frames are
-# identical (even though this is incorrect), so we pass verbose='error' for now
+# identical (even though this is incorrect), so we pass verbose='error'
 raw = mne.io.read_raw_fil(opm_file, verbose="error")
 raw.crop(120, 210).load_data()  # crop for speed
 
@@ -240,7 +240,59 @@
     raw, events, tmin=-0.1, tmax=0.4, baseline=(-0.1, 0.0), verbose="error"
 )
 evoked = epochs.average()
-evoked.plot()
+t_peak = evoked.times[np.argmax(np.std(evoked.copy().pick("meg").data, axis=0))]
+fig = evoked.plot()
+fig.axes[0].axvline(t_peak, color="red", ls="--", lw=1)
+
+# %%
+# Visualizing coregistration
+# --------------------------
+# By design, the sensors in this dataset are already in the scanner RAS coordinate
+# frame. We can thus visualize them in the FreeSurfer MRI coordinate frame by computing
+# the transformation between the FreeSurfer MRI coordinate frame and scanner RAS:
+
+mri = nib.load(subjects_dir / "sub-002" / "mri" / "T1.mgz")
+trans = mri.header.get_vox2ras_tkr() @ np.linalg.inv(mri.affine)
+trans[:3, 3] /= 1000.0  # nibabel uses mm, MNE uses m
+trans = mne.transforms.Transform("head", "mri", trans)
+
+bem = subjects_dir / subject / "bem" / f"{subject}-5120-bem-sol.fif"
+src = subjects_dir / subject / "bem" / f"{subject}-oct-6-src.fif"
+mne.viz.plot_alignment(
+    evoked.info,
+    subjects_dir=subjects_dir,
+    subject=subject,
+    trans=trans,
+    surfaces={"head": 0.1, "inner_skull": 0.2, "white": 1.0},
+    meg=["helmet", "sensors"],
+    verbose="error",
+    bem=bem,
+    src=src,
+)
+
+# %%
+# Plotting the inverse
+# --------------------
+# Now we can compute a forward and inverse:
+
+fwd = mne.make_forward_solution(
+    evoked.info,
+    trans=trans,
+    bem=bem,
+    src=src,
+    verbose=True,
+)
+noise_cov = mne.compute_covariance(epochs, tmax=0)
+inv = mne.minimum_norm.make_inverse_operator(evoked.info, fwd, noise_cov, verbose=True)
+stc = mne.minimum_norm.apply_inverse(
+    evoked, inv, 1.0 / 9.0, method="dSPM", verbose=True
+)
+brain = stc.plot(
+    hemi="split",
+    size=(800, 400),
+    initial_time=t_peak,
+    subjects_dir=subjects_dir,
+)
 
 # %%
 # References