Fix scaling and interp bugs

doc/changes/dev/14461.bugfix.rst

@@ -0,0 +1,6 @@
+Fixed several issues with readers and interpolations, by `Eric Larson`_:
+
+- :func:`mne.io.read_raw_eximia` now returns EEG and EOG data in V instead of µV (values are now 1e-6 times what they were before)
+- :meth:`mne.io.Raw.interpolate_bads` with fNIRS channels no longer clear all of ``info["bads"]`` even with ``reset_bads=False``
+- ``reset_bads=True`` now removes only the channels that were actually interpolated from ``info["bads"]`` instead of all non-excluded bads
+- :meth:`mne.io.Raw.set_montage` stores EEG digitization points in channel order (with the reference first as ident 0) instead of montage order

mne/channels/channels.py

@@ -1067,7 +1067,9 @@ def interpolate_bads(
         Parameters
         ----------
         reset_bads : bool
-            If True, remove the bads from info.
+            If True, remove the interpolated channels from ``info["bads"]``.
+            Bad channels that were not interpolated (e.g., those in ``exclude``
+            or of a channel type not handled by ``method``) are kept.
         mode : str
             Either ``'accurate'`` or ``'fast'``, determines the quality of the
             Legendre polynomial expansion used for interpolation of channels
@@ -1195,6 +1197,10 @@ def interpolate_bads(
             idx = _picks_to_idx(self.info, ch_type, exclude=(), allow_empty=True)
             if len(pick_info(self.info, idx)["bads"]) == 0:
                 method.pop(ch_type)
+        picks = _picks_to_idx(
+            self.info, list(method), exclude=exclude, allow_empty=True
+        )
+        interpolated = set(self.info["bads"]) & {self.ch_names[pi] for pi in picks}
         logger.info("Interpolating bad channels.")
         needs_origin = [key != "seeg" and val != "nan" for key, val in method.items()]
         if any(needs_origin):
@@ -1231,7 +1237,9 @@ def interpolate_bads(
                     "nan-containing channels can be easily excluded from later "
                     "computations."
                 )
-            self.info["bads"] = [ch for ch in self.info["bads"] if ch in exclude]
+            self.info["bads"] = [
+                ch for ch in self.info["bads"] if ch not in interpolated
+            ]
 
         return self
 

mne/channels/interpolation.py

@@ -303,9 +303,6 @@ def _interpolate_bads_nirs(inst, exclude=(), verbose=None):
         closest_raw_idx = picks_nirs[closest_dist_idx]
         inst._data[..., bad_raw_idx, :] = inst._data[..., closest_raw_idx, :]
 
-    # TODO: this seems like a bug because it does not respect reset_bads
-    inst.info["bads"] = [ch for ch in inst.info["bads"] if ch in exclude]
-
     return inst
 
 

mne/channels/montage.py

@@ -12,6 +12,7 @@
 import numpy as np
 
 from .._fiff._digitization import (
+    DigPoint,
     _coord_frame_const,
     _count_points_by_type,
     _ensure_fiducials_head,
@@ -1596,13 +1597,20 @@ def _backcompat_value(pos, ref_pos):
         ch["coord_frame"] = FIFF.FIFFV_COORD_HEAD
     del ch_pos
 
-    # XXX this is probably wrong as it uses the order from the montage
-    # rather than the order of our info['ch_names'] ...
     digpoints = [
-        mnt_head.dig[ii]
+        mnt_head.dig[ii] if name is None else DigPoint(mnt_head.dig[ii], ident=0)
         for ii, name in enumerate(dig_names_use)
-        if name in (info_names_use + ref_names)
+        if name in ref_names
     ]
+    # renumber so enforce_order (which sorts by ident) keeps info["ch_names"] order
+    dig_idx = {name: ii for ii, name in enumerate(dig_names_use)}
+    ch_names_dig = [
+        name
+        for name in dict.fromkeys(info_names_use)
+        if name in dig_idx and name not in ref_names
+    ]
+    for ident, name in enumerate(ch_names_dig, 1):
+        digpoints.append(DigPoint(mnt_head.dig[dig_idx[name]], ident=ident))
 
     # get a copy of the old dig
     if info["dig"] is not None:

mne/channels/tests/test_interpolation.py

@@ -223,9 +223,10 @@ def test_interpolation_meg():
     raw.crop(0, 0.1).load_data().pick(epochs_meg.ch_names)
     raw.info.normalize_proj()
     raw.set_channel_types({raw.ch_names[0]: "stim"}, on_unit_change="ignore")
-    raw.info["bads"] = [raw.ch_names[1]]
+    raw.info["bads"] = raw.ch_names[:2]
     raw.load_data()
     raw.interpolate_bads(mode="fast")
+    assert raw.info["bads"] == raw.ch_names[:1]  # stim not interpolated, stays bad
     del raw
 
     # check that interpolation works for MEG
@@ -267,6 +268,7 @@ def test_interpolation_meg():
     data2 = evoked.interpolate_bads(origin="auto", exclude=["MEG 0121"]).data[pick]
     assert np.corrcoef(data1[0], data2[0])[0, 1] > thresh
     assert np.all(data2[1] == 1e10)
+    assert evoked.info["bads"] == ["MEG 0121"]
 
 
 def _this_interpol(inst, ref_meg=False):
@@ -383,10 +385,12 @@ def test_interpolation_nirs_reordered_picks():
     """Test NIRS interpolation uses the closest donor in raw channel space."""
     raw = _nirs_raw()
     raw.info["bads"] = ["S2_D2 760", "S2_D2 850"]
-    raw.interpolate_bads(
-        method=dict(fnirs="nearest"), origin=(0.0, 0.0, 0.0), verbose=False
-    )
+    kwargs = dict(method=dict(fnirs="nearest"), origin=(0.0, 0.0, 0.0), verbose=False)
+    raw.interpolate_bads(reset_bads=False, **kwargs)
     _assert_nirs_donor(raw)
+    assert raw.info["bads"] == ["S2_D2 760", "S2_D2 850"]
+    raw.interpolate_bads(**kwargs)
+    assert raw.info["bads"] == []
 
 
 def test_interpolation_nirs_epochs():

mne/channels/tests/test_montage.py

@@ -1044,6 +1044,17 @@ def test_set_dig_montage():
         ],
     )
 
+    # EEG dig follows info["ch_names"] order, with the reference first as ident 0
+    montage_rev = make_dig_montage(
+        ch_pos=dict(**dict(reversed(ch_pos.items())), REF=np.full(3, 42)),
+        coord_frame="head",
+    )
+    info = create_info(ch_names, sfreq=1, ch_types="eeg")
+    info.set_montage(montage_rev)
+    eeg_dig = [d for d in info["dig"] if d["kind"] == FIFF.FIFFV_POINT_EEG]
+    assert [d["ident"] for d in eeg_dig] == [0, 1, 2, 3]
+    assert_allclose([d["r"] for d in eeg_dig], np.r_[[[42] * 3], list(ch_pos.values())])
+
 
 def test_set_dig_montage_with_nan_positions():
     """Test that fiducials are not NaN.

mne/io/edf/tests/test_edf.py

@@ -214,16 +214,12 @@ def test_subject_info(tmp_path):
 
 def test_bdf_data():
     """Test reading raw bdf files."""
-    # XXX BDF data for these is around 0.01 when it should be in the uV range,
-    # probably some bug
-    test_scaling = False
     raw_py = _test_raw_reader(
         read_raw_bdf,
         input_fname=bdf_path,
         eog=eog,
         misc=misc,
         exclude=["M2", "IEOG"],
-        test_scaling=test_scaling,
     )
     assert len(raw_py.ch_names) == 71
     raw_py = _test_raw_reader(
@@ -233,7 +229,6 @@ def test_bdf_data():
         eog=eog,
         misc=misc,
         exclude=["M2", "IEOG"],
-        test_scaling=test_scaling,
     )
     assert len(raw_py.ch_names) == 71
     assert "RawBDF" in repr(raw_py)

mne/io/edf/tests/test_gdf.py

@@ -144,7 +144,6 @@ def test_gdf2_data():
         input_fname=gdf2_path.with_name(gdf2_path.name + ".gdf"),
         eog=None,
         misc=None,
-        test_scaling=False,  # XXX this should be True
     )
 
 

mne/io/egi/tests/test_egi.py

@@ -82,7 +82,7 @@ def test_egi_mff_pause(fname, skip_times, event_times):
             raw = _test_raw_reader(
                 read_raw_egi,
                 input_fname=fname,
-                test_scaling=False,  # XXX probably some bug
+                test_scaling=False,  # large DC offsets plus zero-filled skips
                 test_rank="less",
                 events_as_annotations=False,
             )
@@ -167,7 +167,7 @@ def test_io_egi_mff(events_as_annotations):
         input_fname=egi_mff_fname,
         include=include,
         channel_naming="EEG %03d",
-        test_scaling=False,  # XXX probably some bug
+        test_scaling=False,  # real data with transients up to ~140 mV
         events_as_annotations=events_as_annotations,
     )
     assert raw.info["sfreq"] == 1000.0
@@ -260,7 +260,6 @@ def test_io_egi():
         input_fname=egi_fname,
         include=include,
         test_rank="less",
-        test_scaling=False,  # XXX probably some bug
         events_as_annotations=False,
     )
 
@@ -322,7 +321,6 @@ def test_io_egi_pns_mff(tmp_path):
         channel_naming="EEG %03d",
         verbose="error",
         test_rank="less",
-        test_scaling=False,  # XXX probably some bug
     )
     assert names == pns_names
     mat_names = [
@@ -682,7 +680,7 @@ def test_read_event_keys(fname, expected):
         raw = _test_raw_reader(
             read_raw_egi,
             input_fname=fname,
-            test_scaling=False,
+            test_scaling=False,  # large DC offsets plus zero-filled skips
             test_rank="less",
             events_as_annotations=True,
             event_key="cel#",

mne/io/eximia/eximia.py

@@ -110,7 +110,7 @@ def __init__(self, fname, preload=False, verbose=None):
             0.0015259021896696422,
             0.0015259021896696422,
             0.0015259021896696422,
-            0.3814755474174106,
+            0.3814755474174106e-6,
         ]
         ch_names += (
             "Fp1 Fpz Fp2 AF1 AFz AF2 "
@@ -122,7 +122,7 @@ def __init__(self, fname, preload=False, verbose=None):
             "P10 PO3 POz PO4 O1 Oz O2 Iz".split()
         )
         n_eeg = len(ch_names) - len(cals)
-        cals += [0.07629510948348212] * n_eeg
+        cals += [0.07629510948348212e-6] * n_eeg
         ch_types += ["eeg"] * n_eeg
         assert len(ch_names) == n_chan
         info = create_info(ch_names, sfreq, ch_types)

mne/io/eximia/tests/test_eximia.py

@@ -2,7 +2,8 @@
 # License: BSD-3-Clause
 # Copyright the MNE-Python contributors.
 
-from numpy.testing import assert_array_equal
+import numpy as np
+from numpy.testing import assert_allclose, assert_array_equal
 from scipy import io as sio
 
 from mne.datasets.testing import data_path, requires_testing_data
@@ -18,11 +19,7 @@ def test_eximia_nxe():
     fname = testing_path / "eximia" / "test_eximia.nxe"
     raw = read_raw_eximia(fname, preload=True)
     assert "RawEximia" in repr(raw)
-    _test_raw_reader(
-        read_raw_eximia,
-        fname=fname,
-        test_scaling=False,  # XXX probably a scaling problem
-    )
+    _test_raw_reader(read_raw_eximia, fname=fname)
     fname_mat = testing_path / "eximia" / "test_eximia.mat"
     mc = sio.loadmat(fname_mat)
     m_data = mc["data"]
@@ -47,4 +44,6 @@ def test_eximia_nxe():
     ]
     assert ch_types == m_ch_types
 
-    assert_array_equal(m_data, raw._data)
+    stim = np.array(ch_types) == "stim"
+    assert_array_equal(m_data[stim], raw._data[stim])
+    assert_allclose(m_data[~stim] * 1e-6, raw._data[~stim], rtol=1e-12)

mne/io/fieldtrip/tests/helpers.py

@@ -161,6 +161,9 @@ def _get_raw_data(system, drop_extra_chs):
     if system in ["CNT"]:
         raw_data._data = raw_data._data * 1e6
 
+    if system == "eximia":
+        raw_data._data[mne.pick_types(raw_data.info, eeg=True, eog=True)] *= 1e6
+
     if system in ignore_channels_dict:
         raw_data.drop_channels(ignore_channels_dict[system])
 

mne/io/nihon/tests/test_nihon.py

@@ -24,7 +24,7 @@ def test_nihon_eeg():
     fname = data_path / "NihonKohden" / "MB0400FU.EEG"
     raw = read_raw_nihon(fname.as_posix(), preload=True)
     assert "RawNihon" in repr(raw)
-    _test_raw_reader(read_raw_nihon, fname=fname, test_scaling=False)
+    _test_raw_reader(read_raw_nihon, fname=fname)
     fname_edf = data_path / "NihonKohden" / "MB0400FU.EDF"
     raw_edf = read_raw_edf(fname_edf, preload=True)
     raw_edf.drop_channels(["Events/Markers"])

mne/io/nsx/tests/test_nsx.py

@@ -91,7 +91,7 @@ def test_nsx_ver_31():
             input_fname=nsx_31_fname,
             eog=None,
             misc=None,
-            test_scaling=False,
+            test_scaling=False,  # synthetic data with a +/-5 V analog range
             test_rank=False,
         )
     raw_data, _ = raw[:]
@@ -162,7 +162,7 @@ def test_nsx_ver_22():
         input_fname=nsx_22_fname,
         eog=None,
         misc=None,
-        test_scaling=False,  # XXX this should be True
+        test_scaling=False,  # synthetic data with a +/-5 V analog range
         test_rank=False,
     )
     raw_data, _ = raw[:]

mne/time_frequency/tfr.py

@@ -1547,17 +1547,7 @@ def __init__(
                 for k, v in dict(method=method, freqs=freqs).items()
                 if v is None
             ]
-            # TODO when py3.11 is min version, replace if/elif/else block with
-            # classname = inspect.currentframe().f_back.f_code.co_qualname.split(".")[0]
-            _varnames = inspect.currentframe().f_back.f_code.co_varnames
-            if "BaseRaw" in _varnames:
-                classname = "RawTFR"
-            elif "Evoked" in _varnames:
-                classname = "AverageTFR"
-            else:
-                assert "BaseEpochs" in _varnames and "Evoked" not in _varnames
-                classname = "EpochsTFR"
-            # end TODO
+            classname = inspect.currentframe().f_back.f_code.co_qualname.split(".")[0]
             raise ValueError(
                 f"{classname} got unsupported parameter value{_pl(problem)} "
                 f"{' and '.join(problem)}."

mne/viz/epochs.py

@@ -1092,7 +1092,6 @@ def plot_epochs(
     else:
         # make our own copy of the dict
         event_id = dict() if event_id is True else event_id.copy()  # to dict
-        # TODO: when min py=3.9, change to `epochs.event_id | event_id` (maybe).
         # Passed-in event_id should take precedence, i.e., not replace existing
         # keys *or* repeat existing values. For example, if epochs.event_id has
         # a=1 and passed-in event_id has f=1, the second takes precedence.