Fix BEM solving tolerance (#14071)

Author: larsoner

mne/_fiff/tests/test_reference.py

@@ -413,7 +413,7 @@ def test_set_eeg_reference_rest():
     assert 0.006 < exp_var < 0.008
     evoked.set_eeg_reference("REST", forward=fwd)
     exp_var_old = 1 - np.linalg.norm(evoked.data[:, idx] - old) / norm
-    assert 0.005 < exp_var_old <= 0.009
+    assert 0.005 < exp_var_old <= 0.0095
     exp_var = 1 - np.linalg.norm(evoked.data[:, idx] - want) / norm
     assert 0.995 < exp_var <= 1
 

mne/bem.py

@@ -12,7 +12,7 @@
 import shutil
 from collections import OrderedDict
 from copy import deepcopy
-from functools import partial
+from functools import cache, partial
 from pathlib import Path
 
 import numpy as np
@@ -791,23 +791,44 @@ def _one_step(mu, u):
 
 def _fwd_eeg_fit_berg_scherg(m, nterms, nfit):
     """Fit the Berg-Scherg equivalent spherical model dipole parameters."""
+    # The fit depends only on the relative radii and conductivities of the
+    # layers, not the absolute head radius or origin, so cache on those
+    # Using the exact relative-radius ratio also keeps scipy's COBYLA out of a
+    # near-degenerate regime where it fails to converge
+    rel_rads = tuple(float(layer["rel_rad"]) for layer in m["layers"])
+    sigmas = tuple(float(layer["sigma"]) for layer in m["layers"])
+    mu, lambda_, rv = _fit_berg_scherg_cached(rel_rads, sigmas, nterms, nfit)
+
+    m["mu"] = np.array(mu)
+    # This division takes into account the actual conductivities
+    m["lambda"] = np.array(lambda_) / m["layers"][-1]["sigma"]
+    m["nfit"] = nfit
+    return rv
+
+
+@cache
+def _fit_berg_scherg_cached(rel_rads, sigmas, nterms, nfit):
+    """Fit Berg-Scherg params (pure function of relative radii and sigmas)."""
     assert nfit >= 2
+    # Only rel_rad and sigma are read to compute the coefficients and weighting.
+    m = dict(layers=[dict(rel_rad=r, sigma=s) for r, s in zip(rel_rads, sigmas)])
     u = dict(nfit=nfit, nterms=nterms)
 
     # (1) Calculate the coefficients of the true expansion
     u["fn"] = _fwd_eeg_get_multi_sphere_model_coeffs(m, nterms + 1)
 
-    # (2) Calculate the weighting
-    f = min([layer["rad"] for layer in m["layers"]]) / max(
-        [layer["rad"] for layer in m["layers"]]
-    )
+    # (2) Calculate the weighting from the relative-radius ratio
+    f = min(rel_rads) / max(rel_rads)
 
     # correct weighting
     k = np.arange(1, nterms + 1)
     u["w"] = np.sqrt((2.0 * k + 1) * (3.0 * k + 1.0) / k) * np.power(f, (k - 1.0))
     u["w"][-1] = 0
 
-    # Do the nonlinear minimization, constraining mu to the interval [-1, +1]
+    # rhobeg (initial trust-region radius) is ~half the (-1, 1) variable range
+    # rhoend (final radius) sets the resolution of mu. The dipoles sit at radii
+    # proportional to mu (order 1) while the fit's residual var is ~1e-4, so resolving
+    # below that gains nothing and can prevent convergence
     mu_0 = np.zeros(3)
     fun = partial(_one_step, u=u)
     catol = 1e-6
@@ -816,18 +837,13 @@ def _fwd_eeg_fit_berg_scherg(m, nterms, nfit):
     def cons(x):
         return max_ - np.abs(x)
 
-    mu = fmin_cobyla(fun, mu_0, [cons], rhobeg=0.5, rhoend=1e-5, catol=catol)
+    mu = fmin_cobyla(fun, mu_0, [cons], rhobeg=0.5, rhoend=1e-4, catol=catol)
 
     # (6) Do the final step: calculation of the linear parameters
     rv, lambda_ = _compute_linear_parameters(mu, u)
     order = np.argsort(mu)[::-1]
     mu, lambda_ = mu[order], lambda_[order]  # sort: largest mu first
-
-    m["mu"] = mu
-    # This division takes into account the actual conductivities
-    m["lambda"] = lambda_ / m["layers"][-1]["sigma"]
-    m["nfit"] = nfit
-    return rv
+    return tuple(mu), tuple(lambda_), rv
 
 
 @verbose

mne/chpi.py

@@ -565,6 +565,10 @@ def _fit_magnetic_dipole(B_orig, x0, too_close, whitener, coils, guesses):
         idx = np.argmin(res)
         if res[idx] < res0:
             x0 = guesses["rr"][idx]
+    # x0 is the previous time point's coil position (or a better guess-grid point), so a
+    # rhobeg (initial trust-region radius) of 1 mm covers typical between-window motion.
+    # rhoend (final radius) sets the position resolution; 10 um is finer than cHPI
+    # accuracy but still converges in a few dozen evaluations
     x = fmin_cobyla(objective, x0, (), rhobeg=1e-3, rhoend=1e-5, disp=False)
     gof, moment = objective(x, return_moment=True)
     gof = 1.0 - gof / B2

mne/dipole.py

@@ -1354,10 +1354,14 @@ def _fit_dipole(
     )
 
     # Tested minimizers:
-    #    Simplex, BFGS, CG, COBYLA, L-BFGS-B, Powell, SLSQP, TNC
+    #    Simplex, BFGS, CG, COBYLA, COBYQA, L-BFGS-B, Powell, SLSQP, TNC
     # Several were similar, but COBYLA won for having a handy constraint
     # function we can use to ensure we stay inside the inner skull /
-    # smallest sphere
+    # smallest sphere. COBYQA needs ~10x more constraint evaluations, and the
+    # inner-skull constraint is a surface search rather than a cheap formula
+    # rhobeg (initial trust-region radius) is a few times the 2 cm guess grid so the fit
+    # can leave a poorly chosen grid cell. rhoend (final radius, the ``tol`` argument)
+    # sets the position resolution; loosening it to 1e-4 already shifts fits by ~2 mm
     rd_final = fmin_cobyla(
         fun, x0, (constraint,), consargs=(), rhobeg=5e-2, rhoend=rhoend, disp=False
     )