ENH: Add in_argmax mode to better match older MaxFilter (#14059)

Author: larsoner

doc/changes/dev/14059.bugfix.rst

@@ -0,0 +1 @@
+Fixed two bugs in how SSS information is stored and read: :func:`mne.preprocessing.maxwell_filter` recorded which multipole moments it had kept in the wrong order (the filtered data itself was unaffected), and :func:`mne.compute_rank` reported the full MEG rank rather than the reduced SSS rank for data processed with MaxFilter™-3.0, by `Eric Larson`_.

doc/changes/dev/14059.newfeature.rst

@@ -0,0 +1 @@
+Added ``regularize="in_argmax"`` to :func:`mne.preprocessing.maxwell_filter` to allow component selection more similar to that of MaxFilter™-2.2, by `Eric Larson`_.

doc/references.bib

@@ -1025,6 +1025,24 @@ @article{KnuutilaEtAl1993
   year = {1993}
 }
 
+@inproceedings{NenonenEtAl2004,
+  author = {Nenonen, Jukka and Kajola, Matti and Simola, Juha and Ahonen, Antti},
+  booktitle = {Proceedings of the 14th International Conference on Biomagnetism (BIOMAG2004)},
+  pages = {630-631},
+  title = {Total Information of Multichannel {{MEG}} Sensor Arrays},
+  year = {2004}
+}
+
+@article{NenonenEtAl2007,
+  author = {Nenonen, Jukka and Taulu, Samu and Kajola, Matti and Ahonen, Antti},
+  doi = {10.1016/j.ics.2007.01.058},
+  journal = {International Congress Series},
+  pages = {245-248},
+  title = {Total Information Extracted from {{MEG}} Measurements},
+  volume = {1300},
+  year = {2007}
+}
+
 @article{Koles1991,
   author = {Koles, Zoltan J.},
   doi = {10.1016/0013-4694(91)90163-X},

mne/datasets/config.py

@@ -87,7 +87,7 @@
 # update the checksum in the MNE_DATASETS dict below, and change version
 # here: ↓↓↓↓↓↓↓↓
 RELEASES = dict(
-    testing="0.175",
+    testing="0.176",
     misc="0.27",
     phantom_kit="0.2",
     ucl_opm_auditory="0.2",
@@ -115,7 +115,7 @@
 # Testing and misc are at the top as they're updated most often
 MNE_DATASETS["testing"] = dict(
     archive_name=f"{TESTING_VERSIONED}.tar.gz",
-    hash="md5:0934f219e59aa727988b95352cfc1fcb",
+    hash="md5:4e5f515c71b6897766934170d1d94a7f",
     url=(
         "https://codeload.github.com/mne-tools/mne-testing-data/"
         f"tar.gz/{RELEASES['testing']}"

mne/preprocessing/maxwell.py

@@ -1423,9 +1423,14 @@ def _regularize(
     n_in = _get_n_moments(int_order)
     n_out = S_decomp.shape[1] - n_in
     t_str = f"{t:8.3f}"
-    if regularize is not None:  # regularize='in'
+    if regularize is not None:  # regularize='in' or 'in_argmax'
         in_removes, out_removes = _regularize_in(
-            int_order, ext_order, S_decomp, mag_or_fine, extended_remove
+            int_order,
+            ext_order,
+            S_decomp,
+            mag_or_fine,
+            extended_remove,
+            argmax=regularize == "in_argmax",
         )
     else:
         in_removes = []
@@ -1510,10 +1515,9 @@ def _get_mf_picks_fix_mags(info, int_order, ext_order, ignore_ref=False, verbose
 
 def _check_regularize(regularize):
     """Ensure regularize is valid."""
-    if not (
-        regularize is None or (isinstance(regularize, str) and regularize in ("in",))
-    ):
-        raise ValueError('regularize must be None or "in"')
+    _validate_type(regularize, (str, None), "regularize")
+    if regularize is not None:
+        _check_option("regularize", regularize, ("in", "in_argmax"), extra="when str")
 
 
 def _check_usable(inst, ignore_ref):
@@ -1735,7 +1739,6 @@ def _sss_basis(exp, all_coils):
     sin_az = rmags[:, 1] / r_xy  # sin(phi)
     sin_az[z_only] = 0.0
     # Appropriate vector spherical harmonics terms
-    #  JNE 2012-02-08: modified alm -> 2*alm, blm -> -2*blm
     r_nn2 = r_n.copy()
     r_nn1 = 1.0 / (r_n * r_n)
     S_tot = np.empty((n_coils, n_in + n_out), np.float64)
@@ -1789,7 +1792,13 @@ def _sss_basis(exp, all_coils):
             sin_order = np.sin(ord_phi)
             cos_order = np.cos(ord_phi)
             mult /= np.sqrt((degree - order + 1) * (degree + order))
-            factor = mult * np.sqrt(2)  # equivalence fix (MF uses 2.)
+            # √2 keeps the real basis orthonormal. MaxFilter's real
+            # ("even-odd") coefficients instead absorb a factor of 2 relative to
+            # the complex ones, so our order != 0 columns are its columns
+            # divided by √2. Being a per-column scaling, this cancels in
+            # S_in @ pinv(S_tot); it only matters where column norms do, i.e.
+            # when regularizing (see _regularize_in).
+            factor = mult * np.sqrt(2)
 
             # Real
             idx = _deg_ord_idx(degree, order)
@@ -1952,6 +1961,25 @@ def _get_degrees_orders(order):
     return degrees, orders
 
 
+def _mne_ord_to_mf_idx(order):
+    """Get indices reordering a moment array from our order into MaxFilter's.
+
+    We index moments by signed order ascending within each degree (see
+    ``_deg_ord_idx``, i.e. -degree, ..., 0, ..., +degree). MaxFilter instead
+    walks its even-odd formulation, m = 0, ..., 2 * degree, where m = 0 is the
+    zero order, odd m is the real part of order (m + 1) // 2, and even m > 0 is
+    the imaginary part (our negative order).
+    """
+    return np.array(
+        [
+            _deg_ord_idx(degree, order_)
+            for degree in range(1, order + 1)
+            for order_ in [0, *(s * m for m in range(1, degree + 1) for s in (1, -1))]
+        ],
+        int,
+    )
+
+
 def _alegendre_deriv(order, degree, val):
     """Compute the derivative of the associated Legendre polynomial at a value.
 
@@ -1996,8 +2024,11 @@ def _bases_complex_to_real(complex_tot, int_order, ext_order):
                 idx_neg = _deg_ord_idx(deg, -order)
                 real[:, idx_pos] = _sh_complex_to_real(comp[:, idx_pos], order)
                 if order != 0:
-                    # This extra mult factor baffles me a bit, but it works
-                    # in round-trip testing, so we'll keep it :(
+                    # idx_neg holds (-1)**order * conj(comp[:, idx_pos]) (the
+                    # conjugation property), and _sh_complex_to_real takes
+                    # √2·imag for order < 0, flipping the sign again. This
+                    # mult == -(-1)**order undoes both, leaving
+                    # real[:, idx_neg] == √2 · imag(comp[:, idx_pos]).
                     mult = -1 if order % 2 == 0 else 1
                     real[:, idx_neg] = mult * _sh_complex_to_real(
                         comp[:, idx_neg], -order
@@ -2097,6 +2128,14 @@ def _update_sss_info(
         raw.info["maxshield"] = False
     components = np.zeros(n_in + n_out + len(extended_proj)).astype("int32")
     components[reg_moments] = 1
+    # MaxFilter lays the moments out in its even-odd order rather than ours, so
+    # reindex both expansions before writing them out (nfree is a count, so it
+    # is unaffected). Any extended (eSSS) entries are appended after the two
+    # expansions and have no degree/order structure, so they stay put.
+    components[:n_in] = components[:n_in][_mne_ord_to_mf_idx(int_order)]
+    components[n_in : n_in + n_out] = components[n_in : n_in + n_out][
+        _mne_ord_to_mf_idx(ext_order)
+    ]
     sss_info_dict = dict(
         in_order=int_order,
         out_order=ext_order,
@@ -2400,7 +2439,9 @@ def _regularize_out(int_order, ext_order, mag_or_fine, extended_remove):
     return list(range(n_in, n_in + 3 * remove_homog)) + extended_remove
 
 
-def _regularize_in(int_order, ext_order, S_decomp, mag_or_fine, extended_remove):
+def _regularize_in(
+    int_order, ext_order, S_decomp, mag_or_fine, extended_remove, *, argmax=False
+):
     """Regularize basis set using idealized SNR measure."""
     n_in, n_out = _get_n_moments([int_order, ext_order])
 
@@ -2424,22 +2465,7 @@ def _regularize_in(int_order, ext_order, S_decomp, mag_or_fine, extended_remove)
     # if plot:
     #     import matplotlib.pyplot as plt
     #     fig, axs = plt.subplots(3, figsize=[6, 12])
-    #     plot_ord = np.empty(n_in, int)
-    #     plot_ord.fill(-1)
-    #     count = 0
-    #     # Reorder plot to match MF
-    #     for degree in range(1, int_order + 1):
-    #         for order in range(0, degree + 1):
-    #             assert plot_ord[count] == -1
-    #             plot_ord[count] = _deg_ord_idx(degree, order)
-    #             count += 1
-    #             if order > 0:
-    #                 assert plot_ord[count] == -1
-    #                 plot_ord[count] = _deg_ord_idx(degree, -order)
-    #                 count += 1
-    #     assert count == n_in
-    #     assert (plot_ord >= 0).all()
-    #     assert len(np.unique(plot_ord)) == n_in
+    #     plot_ord = _mne_ord_to_mf_idx(int_order)  # reorder plot to match MF
     noise_lev = 5e-13  # noise level in T/m
     noise_lev *= noise_lev  # effectively what would happen by earlier multiply
     for ii in range(n_in):
@@ -2455,8 +2481,20 @@ def _regularize_in(int_order, ext_order, S_decomp, mag_or_fine, extended_remove)
         eta_lm_sq = eta_lm_sq.sum(axis=1)
         eta_lm_sq *= noise_lev
 
-        # Mysterious scale factors to match MF, likely due to differences
-        # in the basis normalizations...
+        # Together these scale snr by 400 (order != 0) / 200 (order == 0). Only
+        # the ratio of 2 is derivable: our order != 0 columns are MaxFilter's
+        # divided by √2 (see _sss_basis), making eta_lm_sq (a squared pinv)
+        # 2x its value there, which doubling order == 0 evens out.
+        #
+        # The shared factor of 100 is empirical, and puts I_tots in the ~480
+        # bits/sample range reported for this array in Nenonen et al. 2007 (Int
+        # Congr Ser 1300:245). It cannot affect *which* components are dropped,
+        # since remove_order below uses argmin(snr), but it does move the cut:
+        # every snr here is >> 1, so I_tots ~ 0.5*sum(log2(snr)) and a global
+        # scale adds 0.5*n_keepers*log2(scale) -- a ramp, not an offset, as
+        # n_keepers shrinks each iteration. It is nonetheless well enough
+        # calibrated that we reproduce MaxFilter 3.0's selection; see
+        # test_regularization_mf3.
         eta_lm_sq[orders[in_keepers] == 0] *= 2
         eta_lm_sq *= 0.0025
         snr = a_lm_sq[in_keepers] / eta_lm_sq
@@ -2474,12 +2512,17 @@ def _regularize_in(int_order, ext_order, S_decomp, mag_or_fine, extended_remove)
     #     axs[1].set(ylabel='Information', xlabel='Iteration')
     #     axs[2].plot(eigs[:, 0] / eigs[:, 1])
     #     axs[2].set(ylabel='Condition', xlabel='Iteration')
-    # Pick the components that give at least 98% of max info
-    # This is done because the curves can be quite flat, and we err on the
-    # side of including rather than excluding components
+    # argmax is what MaxFilter 2.2 does, and what Nenonen et al. 2007 (Int Congr
+    # Ser 1300:245) describes: "Iteration is stopped when the maximum Itot is
+    # found". Over 27 expansion origins it reproduces 2.2's nfree exactly 16/27
+    # times against 4/27 for the 98% rule, which is 3.0's and which we track
+    # closely instead (see test_regularization_mf3).
     if n_in:
         max_info = np.max(I_tots)
-        lim_idx = np.where(I_tots >= 0.98 * max_info)[0][0]
+        if argmax:
+            lim_idx = np.argmax(I_tots)
+        else:
+            lim_idx = np.where(I_tots >= 0.98 * max_info)[0][0]
         in_removes = remove_order[:lim_idx]
         for ii, ri in enumerate(in_removes):
             eig = eigs[ii]
@@ -2514,10 +2557,29 @@ def _compute_sphere_activation_in(degrees):
     rho_i : float
         The current density.
 
+    Notes
+    -----
+    The model is that of :footcite:`NenonenEtAl2004`, which is also where the
+    constants below and ``mag_scale=100`` come from.
+
     References
     ----------
     .. footbibliography::
     """
+    # Per NenonenEtAl2004: rho_i is Knuutila's 0.6 µA/m²/√Hz random brain
+    # current density (the 100 fT below inverts to 0.591), r_in is its
+    # "integration radius of 80 mm", and mag_scale=100 is its noise ratio,
+    # 3 fT/√Hz mag vs 4.8 fT/√Hz grad over a 16 mm baseline (3e-15 T
+    # vs 3.0e-13 T/m). rho_i and the noise are both spectral densities, so this
+    # is bandwidth free -- the 100 fT is only the RMS it gives over 0.1-30 Hz,
+    # and applying √bandwidth to one but not the other measurably hurts.
+    #
+    # We use the *surface* rather than the volume variant, though that paper
+    # describes a volume current density: the two differ by an l-dependent
+    # R_in²/(l+3)², and surface reproduces MaxFilter's component selection far
+    # better. r_in=0.080 is likewise a sharp empirical optimum across expansion
+    # origins, and making the sphere aware of the origin's offset from the head
+    # centre does not help.
     r_in = 0.080  # radius of the randomly-activated sphere
 
     # set the observation point r=r_s, az=el=0, so we can just look at m=0 term