Merge branch 'mne-tools:main' into gsoc_commits

Author: shresth-keshari

.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
 

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`_.

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`_.

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

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)

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
 
 

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],

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

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"

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