Merge branch 'mne-tools:main' into gsoc_commits

shresth-keshari · web-flow · commit 7a401c4accc7 · 2025-04-01T14:06:32.000+05:30
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
@@ -1,7 +1,7 @@
 repos:
   # Ruff mne
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: v0.9.10
+    rev: v0.11.0
     hooks:
       - id: ruff
         name: ruff lint mne
@@ -33,7 +33,7 @@ repos:
 
   # yamllint
   - repo: https://github.com/adrienverge/yamllint.git
-    rev: v1.35.1
+    rev: v1.36.2
     hooks:
       - id: yamllint
         args: [--strict, -c, .yamllint.yml]
@@ -82,7 +82,7 @@ repos:
 
   # zizmor
   - repo: https://github.com/woodruffw/zizmor-pre-commit
-    rev: v1.4.1
+    rev: v1.5.1
     hooks:
       - id: zizmor
 
diff --git a/doc/changes/devel/13174.bugfix.rst b/doc/changes/devel/13174.bugfix.rst
@@ -0,0 +1 @@
+Fix bug where :func:`mne.export.export_raw` might allocate huge intermediate arrays unnecessarily, when padding data blocks during export to EDF format, by `Daniel McCloy`_.
diff --git a/doc/changes/devel/13178.bugfix.rst b/doc/changes/devel/13178.bugfix.rst
@@ -0,0 +1 @@
+Fix bug with least-squares fitting of head origin using digitization points in :func:`mne.preprocessing.maxwell_filter`, by `Eric Larson`_.
diff --git a/doc/development/contributing.rst b/doc/development/contributing.rst
@@ -291,7 +291,7 @@ Now we'll remove the *stable* version of MNE-Python and replace it with the
 the correct environment first (``conda activate mnedev``), and then do::
 
     $ cd $INSTALL_LOCATION/mne-python    # make sure we're in the right folder
-    $ conda remove --force mne  # the --force avoids dependency checking
+    $ conda remove --force mne-base  # the --force avoids dependency checking
     $ pip install -e .
 
 The command ``pip install -e .`` installs a python module into the current
diff --git a/mne/_fiff/tests/test_reference.py b/mne/_fiff/tests/test_reference.py
@@ -312,7 +312,11 @@ def test_set_eeg_reference_rest():
     same = [raw.ch_names.index(raw.info["bads"][0])]
     picks = np.setdiff1d(np.arange(len(raw.ch_names)), same)
     trans = None
-    sphere = make_sphere_model("auto", "auto", raw.info)
+    # Use fixed values from old sphere fit to reduce lines changed with fixed algorithm
+    sphere = make_sphere_model(
+        [-0.00413508, 0.01598787, 0.05175598],
+        0.09100286249131773,
+    )
     src = setup_volume_source_space(pos=20.0, sphere=sphere, exclude=30.0)
     assert src[0]["nuse"] == 223  # low but fast
     fwd = make_forward_solution(raw.info, trans, src, sphere)
diff --git a/mne/bem.py b/mne/bem.py
@@ -1073,10 +1073,10 @@ def get_fitting_dig(info, dig_kinds="auto", exclude_frontal=True, verbose=None):
 
 
 @verbose
-def _fit_sphere_to_headshape(info, dig_kinds, verbose=None):
+def _fit_sphere_to_headshape(info, dig_kinds, *, verbose=None):
     """Fit a sphere to the given head shape."""
     hsp = get_fitting_dig(info, dig_kinds)
-    radius, origin_head = _fit_sphere(np.array(hsp), disp=False)
+    radius, origin_head = _fit_sphere(np.array(hsp))
     # compute origin in device coordinates
     dev_head_t = info["dev_head_t"]
     if dev_head_t is None:
@@ -1105,36 +1105,16 @@ def _fit_sphere_to_headshape(info, dig_kinds, verbose=None):
     return radius, origin_head, origin_device
 
 
-def _fit_sphere(points, disp="auto"):
+def _fit_sphere(points):
     """Fit a sphere to an arbitrary set of points."""
-    if isinstance(disp, str) and disp == "auto":
-        disp = True if logger.level <= 20 else False
-    # initial guess for center and radius
-    radii = (np.max(points, axis=1) - np.min(points, axis=1)) / 2.0
-    radius_init = radii.mean()
-    center_init = np.median(points, axis=0)
-
-    # optimization
-    x0 = np.concatenate([center_init, [radius_init]])
-
-    def cost_fun(center_rad):
-        d = np.linalg.norm(points - center_rad[:3], axis=1) - center_rad[3]
-        d *= d
-        return d.sum()
-
-    def constraint(center_rad):
-        return center_rad[3]  # radius must be >= 0
-
-    x_opt = fmin_cobyla(
-        cost_fun,
-        x0,
-        constraint,
-        rhobeg=radius_init,
-        rhoend=radius_init * 1e-6,
-        disp=disp,
-    )
-
-    origin, radius = x_opt[:3], x_opt[3]
+    # linear least-squares sphere fit, see for example
+    # https://stackoverflow.com/a/78909044
+    # TODO: At some point we should maybe reject outliers first...
+    A = np.c_[2 * points, np.ones((len(points), 1))]
+    b = (points**2).sum(axis=1)
+    x, _, _, _ = np.linalg.lstsq(A, b, rcond=1e-6)
+    origin = x[:3]
+    radius = np.sqrt(x[0] ** 2 + x[1] ** 2 + x[2] ** 2 + x[3])
     return radius, origin
 
 
diff --git a/mne/channels/tests/test_montage.py b/mne/channels/tests/test_montage.py
@@ -1453,7 +1453,7 @@ def renamer(x):
     assert (orig_pos != new_pos).all()
 
     r0 = _fit_sphere(new_pos)[1]
-    assert_allclose(r0, [-0.001021, 0.014554, 0.041404], atol=1e-4)
+    assert_allclose(r0, [-0.001043, 0.01469, 0.041448], atol=1e-4)
     # spot check: Fp1 and Fpz
     assert_allclose(
         new_pos[:2],
diff --git a/mne/dipole.py b/mne/dipole.py
@@ -1573,7 +1573,7 @@ def fit_dipole(
         # Find the best-fitting sphere
         inner_skull = _bem_find_surface(bem, "inner_skull")
         inner_skull = inner_skull.copy()
-        R, r0 = _fit_sphere(inner_skull["rr"], disp=False)
+        R, r0 = _fit_sphere(inner_skull["rr"])
         # r0 back to head frame for logging
         r0 = apply_trans(mri_head_t["trans"], r0[np.newaxis, :])[0]
         inner_skull["r0"] = r0
diff --git a/mne/export/_edf.py b/mne/export/_edf.py
@@ -66,7 +66,18 @@ def _export_raw(fname, raw, physical_range, add_ch_type):
                 f"{pad_width / sfreq:.3g} seconds of edge values were appended to all "
                 "channels when writing the final block."
             )
-            data = np.pad(data, (0, int(pad_width)), "edge")
+            orig_shape = data.shape
+            data = np.pad(
+                data,
+                (
+                    (0, 0),
+                    (0, int(pad_width)),
+                ),
+                "edge",
+            )
+            assert data.shape[0] == orig_shape[0]
+            assert data.shape[1] > orig_shape[1]
+
             annotations.append(
                 EdfAnnotation(
                     raw.times[-1] + 1 / sfreq, pad_width / sfreq, "BAD_ACQ_SKIP"
diff --git a/mne/forward/tests/test_make_forward.py b/mne/forward/tests/test_make_forward.py
@@ -719,7 +719,7 @@ def test_make_forward_dipole(tmp_path):
     # Make sure each coordinate is close to reference
     # NB tolerance should be set relative to snr of simulated evoked!
     assert_allclose(
-        dip_fit.pos, dip_test.pos, rtol=0, atol=1e-2, err_msg="position mismatch"
+        dip_fit.pos, dip_test.pos, rtol=0, atol=1.3e-2, err_msg="position mismatch"
     )
     assert dist < 1e-2  # within 1 cm
     assert corr > 0.985
diff --git a/mne/io/curry/tests/test_curry.py b/mne/io/curry/tests/test_curry.py
@@ -280,7 +280,7 @@ def test_read_raw_curry_rfDC(fname, tol, mock_dev_head_t, tmp_path):
     assert n_other == 15
     pos = np.array(pos)
     nn = np.array(nn)
-    rad, origin = _fit_sphere(pos, disp=False)
+    rad, origin = _fit_sphere(pos)
     assert 0.11 < rad < 0.13
     pos -= origin
     pos /= np.linalg.norm(pos, axis=1, keepdims=True)
diff --git a/mne/preprocessing/tests/test_fine_cal.py b/mne/preprocessing/tests/test_fine_cal.py
@@ -77,7 +77,7 @@ def test_compute_fine_cal(kind):
         angle_limit = 5
         gwoma = [66, 68]
         ggoma = [55, 150]
-        ggwma = [62, 86]
+        ggwma = [61, 86]
         sfs = [26, 27, 61, 63, 61, 63, 68, 70]
         cl3 = [0.6, 0.7]
     else:
@@ -222,13 +222,13 @@ def test_fine_cal_systems(system, tmp_path):
         err_limit = 500
         n_ref = 3
         corrs = (0.58, 0.61, 0.57)
-        sfs = [0.9, 1.1, 2.3, 2.8]
+        sfs = [0.9, 1.1, 2.2, 2.8]
         corr_tol = 0.3
     elif system == "ctf":
         raw = read_raw_ctf(ctf_fname_continuous).crop(0, 1)
         raw.apply_gradient_compensation(0)
         angle_limit = 170
-        err_limit = 6000
+        err_limit = 12600
         n_ref = 28
         corrs = (0.19, 0.41, 0.49)
         sfs = [0.5, 0.7, 0.9, 1.55]
@@ -242,8 +242,8 @@ def test_fine_cal_systems(system, tmp_path):
         err_limit = 15
         int_order = 5
         corrs = (0.13, 0.0, 0.12)
-        sfs = [4, 5, 125, 145]
-        corr_tol = 0.15
+        sfs = [4, 5, 125, 155]
+        corr_tol = 0.34
     else:
         assert system == "triux", f"Unknown system {system}"
         raw = read_raw_fif(tri_fname)
diff --git a/mne/preprocessing/tests/test_maxwell.py b/mne/preprocessing/tests/test_maxwell.py
@@ -173,6 +173,7 @@
 
 def _assert_n_free(raw_sss, lower, upper=None):
     """Check the DOF."""
+    __tracebackhide__ = True
     upper = lower if upper is None else upper
     n_free = raw_sss.info["proc_history"][0]["max_info"]["sss_info"]["nfree"]
     assert lower <= n_free <= upper, f"nfree fail: {lower} <= {n_free} <= {upper}"
@@ -366,8 +367,10 @@ def test_other_systems():
     elp_path = kit_dir / "test_elp.txt"
     hsp_path = kit_dir / "test_hsp.txt"
     raw_kit = read_raw_kit(sqd_path, str(mrk_path), str(elp_path), str(hsp_path))
-    with pytest.warns(RuntimeWarning, match="fit"):
-        pytest.raises(RuntimeError, maxwell_filter, raw_kit)
+    with (
+        pytest.raises(NotImplementedError, match="Cannot create forward solution with"),
+    ):
+        maxwell_filter(raw_kit, verbose=True)
     with catch_logging() as log:
         raw_sss = maxwell_filter(
             raw_kit, origin=(0.0, 0.0, 0.04), ignore_ref=True, verbose=True
@@ -378,24 +381,21 @@ def test_other_systems():
         raw_kit, origin=(0.0, 0.0, 0.04), ignore_ref=True, mag_scale="auto"
     )
     assert_allclose(raw_sss._data, raw_sss_auto._data)
-    # The KIT origin fit is terrible
-    with pytest.warns(RuntimeWarning, match="more than 20 mm"):
-        with catch_logging() as log:
-            pytest.raises(
-                RuntimeError, maxwell_filter, raw_kit, ignore_ref=True, regularize=None
-            )  # bad condition
-            raw_sss = maxwell_filter(
-                raw_kit,
-                origin="auto",
-                ignore_ref=True,
-                bad_condition="info",
-                verbose=True,
-            )
+    with catch_logging() as log:
+        pytest.raises(
+            RuntimeError, maxwell_filter, raw_kit, ignore_ref=True, regularize=None
+        )  # bad condition
+        raw_sss = maxwell_filter(
+            raw_kit,
+            origin="auto",
+            ignore_ref=True,
+            bad_condition="info",
+            verbose=True,
+        )
     log = log.getvalue()
-    assert "badly conditioned" in log
-    assert "more than 20 mm from" in log
-    # fits can differ slightly based on scipy version, so be lenient here
-    _assert_n_free(raw_sss, 28, 34)  # bad origin == brutal reg
+    assert "badly conditioned" not in log
+    assert "more than 20 mm from" not in log
+    _assert_n_free(raw_sss, 67, 67)
     # Let's set the origin
     with catch_logging() as log:
         raw_sss = maxwell_filter(
diff --git a/mne/report/report.py b/mne/report/report.py
@@ -1613,7 +1613,7 @@ def add_trans(
         title : str
             The title to add.
         subject : str | None
-            The name of the FreeSurfer subject the ``trans```` belong to. The
+            The name of the FreeSurfer subject the ``trans`` belongs to. The
             name is not stored with the ``trans`` and therefore needs to be
             specified. If ``None``, will use the value of ``subject`` passed on
             report creation.
diff --git a/mne/simulation/tests/test_raw.py b/mne/simulation/tests/test_raw.py
@@ -220,7 +220,11 @@ def raw_data():
 
     src = read_source_spaces(src_fname)
     trans = read_trans(trans_fname)
-    sphere = make_sphere_model("auto", "auto", raw.info)
+    # Use fixed values from old sphere fit to reduce lines changed with fixed algorithm
+    sphere = make_sphere_model(
+        [-0.00413508, 0.01598787, 0.05175598],
+        0.09100286249131773,
+    )
     stc = _make_stc(raw, src)
     return raw, src, stc, trans, sphere
 
@@ -385,9 +389,9 @@ def test_simulate_raw_bem(raw_data):
     locs = np.concatenate([s["rr"][s["vertno"]] for s in src])
     tmax = (len(locs) - 1) / raw.info["sfreq"]
     cov = make_ad_hoc_cov(raw.info)
-    # The tolerance for the BEM is surprisingly high (28) but I get the same
+    # 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, 2), (raw_sim_bem, bem_fname, 31)):
+    for use_raw, bem, tol in ((raw_sim_sph, sphere, 4), (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()
diff --git a/mne/surface.py b/mne/surface.py
@@ -220,7 +220,7 @@ def get_meg_helmet_surf(info, trans=None, *, verbose=None):
         )
         hull = ConvexHull(rr)
         rr = rr[np.unique(hull.simplices)]
-        R, center = _fit_sphere(rr, disp=False)
+        R, center = _fit_sphere(rr)
         sph = _cart_to_sph(rr - center)[:, 1:]
         # add a point at the front of the helmet (where the face should be):
         # 90 deg az and maximal el (down from Z/up axis)
diff --git a/mne/tests/test_bem.py b/mne/tests/test_bem.py
@@ -481,17 +481,16 @@ def test_fit_sphere_to_headshape():
     assert_allclose(r, r2, atol=1e-7)
     assert_allclose(oh, oh2, atol=1e-7)
     assert_allclose(od, od2, atol=1e-7)
-    # this one should pass, 1 EXTRA point and 3 EEG (but the fit is terrible)
+    # this one should pass, 1 EXTRA point and 3 EEG
     info = Info(dig=dig[:7], dev_head_t=dev_head_t)
-    with (
-        _record_warnings(),
-        pytest.warns(RuntimeWarning, match="Estimated head radius"),
-    ):
+    with pytest.warns(RuntimeWarning, match="fitting may be inaccurate"):
         r, oh, od = fit_sphere_to_headshape(info, units="m")
-    # this one should fail, 1 EXTRA point and 3 EEG (but the fit is terrible)
+    # this one should fail
     info = Info(dig=dig[:6], dev_head_t=dev_head_t)
-    pytest.raises(ValueError, fit_sphere_to_headshape, info, units="m")
-    pytest.raises(TypeError, fit_sphere_to_headshape, 1, units="m")
+    with pytest.raises(ValueError, match="at least 4"):
+        fit_sphere_to_headshape(info, units="m")
+    with pytest.raises(TypeError, match="Info"):
+        fit_sphere_to_headshape(1, units="m")
 
 
 @pytest.mark.slowtest  # ~2 min on Azure Windows
diff --git a/mne/transforms.py b/mne/transforms.py
@@ -1134,10 +1134,10 @@ def fit(
         if center:
             logger.info("    Centering data")
             hsp = np.array([p for p in source if not (p[2] < -1e-6 and p[1] > 1e-6)])
-            src_center = _fit_sphere(hsp, disp=False)[1]
+            src_center = _fit_sphere(hsp)[1]
             source = source - src_center
             hsp = np.array([p for p in destination if not (p[2] < 0 and p[1] > 0)])
-            dest_center = _fit_sphere(hsp, disp=False)[1]
+            dest_center = _fit_sphere(hsp)[1]
             destination = destination - dest_center
             logger.info(
                 "    Using centers {np.array_str(src_center, None, 3)} -> "
diff --git a/tools/install_pre_requirements.sh b/tools/install_pre_requirements.sh
@@ -15,18 +15,19 @@ echo "PyQt6 and scientific-python-nightly-wheels dependencies"
 python -m pip install $STD_ARGS pip setuptools packaging \
 	threadpoolctl cycler fonttools kiwisolver pyparsing pillow python-dateutil \
 	patsy pytz tzdata nibabel tqdm trx-python joblib numexpr "$QT_BINDING" \
-	py-cpuinfo blosc2 hatchling
+	py-cpuinfo blosc2 hatchling "h5py>=3.12.1"  # TODO move h5py back to nightlies when possible
 echo "NumPy/SciPy/pandas etc."
 python -m pip uninstall -yq numpy
 python -m pip install --upgrade matplotlib  # TODO: Until https://github.com/matplotlib/matplotlib/pull/29427 lands
 python -m pip install $STD_ARGS --only-binary ":all:" --default-timeout=60 \
 	--index-url "https://pypi.anaconda.org/scientific-python-nightly-wheels/simple" \
 	"numpy>=2.1.0.dev0" "scikit-learn>=1.6.dev0" "scipy>=1.15.0.dev0" \
 	"pandas>=3.0.0.dev0" \
-	"h5py>=3.12.1" "dipy>=1.10.0.dev0"
-# TODO: should have above: "pyarrow>=20.0.0.dev0" "tables>=3.10.2.dev0"
-# https://github.com/apache/arrow/issues/40216#issuecomment-2678899168
+	"dipy>=1.10.0.dev0" \
+    "pyarrow>=20.0.0.dev0"
+# TODO: should have above: "tables>=3.10.2.dev0" "h5py>=3.12.1"
 # https://github.com/PyTables/PyTables/issues/1115
+# https://github.com/h5py/h5py/issues/2563
 
 # statsmodels requires formulaic@main so we need to use --extra-index-url
 echo "statsmodels"

Original file line number	Diff line number	Diff line change
`@@ -0,0 +1 @@`
	`1`	+Fix bug where :func:`mne.export.export_raw` might allocate huge intermediate arrays unnecessarily, when padding data blocks during export to EDF format, by `Daniel McCloy`_.
Original file line number	Diff line number	Diff line change
`@@ -0,0 +1 @@`
	`1`	+Fix bug with least-squares fitting of head origin using digitization points in :func:`mne.preprocessing.maxwell_filter`, by `Eric Larson`_.
Original file line number	Diff line number	Diff line change
`@@ -719,7 +719,7 @@ def test_make_forward_dipole(tmp_path):`
`719`	`719`	`# Make sure each coordinate is close to reference`
`720`	`720`	`# NB tolerance should be set relative to snr of simulated evoked!`
`721`	`721`	`assert_allclose(`
`722`		`- dip_fit.pos, dip_test.pos, rtol=0, atol=1e-2, err_msg="position mismatch"`
	`722`	`+ dip_fit.pos, dip_test.pos, rtol=0, atol=1.3e-2, err_msg="position mismatch"`
`723`	`723`	`)`
`724`	`724`	`assert dist < 1e-2 # within 1 cm`
`725`	`725`	`assert corr > 0.985`