From 74ec8481ffd49c05f3e8700e9af30a1d73a16c73 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 15 Jul 2026 14:42:05 -0400 Subject: [PATCH 1/7] ENH: Add in_argmax mode to better match older MaxFilter --- doc/changes/dev/14058.bugfix.rst | 1 + doc/changes/dev/14058.newfeature.rst | 1 + doc/references.bib | 18 +++ mne/datasets/config.py | 4 +- mne/preprocessing/maxwell.py | 128 +++++++++++++++----- mne/preprocessing/tests/test_maxwell.py | 153 ++++++++++++++++++++++-- mne/rank.py | 22 ++-- mne/utils/docs.py | 18 ++- 8 files changed, 287 insertions(+), 58 deletions(-) create mode 100644 doc/changes/dev/14058.bugfix.rst create mode 100644 doc/changes/dev/14058.newfeature.rst diff --git a/doc/changes/dev/14058.bugfix.rst b/doc/changes/dev/14058.bugfix.rst new file mode 100644 index 00000000000..ca75d7760f3 --- /dev/null +++ b/doc/changes/dev/14058.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`_. diff --git a/doc/changes/dev/14058.newfeature.rst b/doc/changes/dev/14058.newfeature.rst new file mode 100644 index 00000000000..61a40d5fcfc --- /dev/null +++ b/doc/changes/dev/14058.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`_. diff --git a/doc/references.bib b/doc/references.bib index 1020008ee70..9a489ee6a63 100644 --- a/doc/references.bib +++ b/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}, diff --git a/mne/datasets/config.py b/mne/datasets/config.py index 0ea542fdca0..5599c20e5d4 100644 --- a/mne/datasets/config.py +++ b/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.173", + 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:0973c25cd0e3ca43a75ea6953ace1daa", + hash="md5:4e5f515c71b6897766934170d1d94a7f", url=( "https://codeload.github.com/mne-tools/mne-testing-data/" f"tar.gz/{RELEASES['testing']}" diff --git a/mne/preprocessing/maxwell.py b/mne/preprocessing/maxwell.py index 7e649acc3cf..6348c6b92a5 100644 --- a/mne/preprocessing/maxwell.py +++ b/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 diff --git a/mne/preprocessing/tests/test_maxwell.py b/mne/preprocessing/tests/test_maxwell.py index d5988825436..768dec3b147 100644 --- a/mne/preprocessing/tests/test_maxwell.py +++ b/mne/preprocessing/tests/test_maxwell.py @@ -42,7 +42,10 @@ from mne.preprocessing.maxwell import ( _bases_complex_to_real, _bases_real_to_complex, + _compute_sphere_activation_in, + _get_degrees_orders, _get_n_moments, + _mne_ord_to_mf_idx, _prep_mf_coils, _sh_complex_to_real, _sh_negate, @@ -70,6 +73,10 @@ sss_nonstd_fname = sss_path / "test_move_anon_nonStdOrigin_raw_sss.fif" sss_bad_recon_fname = sss_path / "test_move_anon_badRecon_raw_sss.fif" sss_reg_in_fname = sss_path / "test_move_anon_regIn_raw_sss.fif" +sss_mf3_reg_in_fname = sss_path / "test_move_anon_mf3_regIn_raw_sss.fif" +sss_mf3_reg_in_nonstd_fname = ( + sss_path / "test_move_anon_mf3_regIn_nonStdOrigin_raw_sss.fif" +) sss_fine_cal_fname = sss_path / "test_move_anon_fineCal_raw_sss.fif" sss_ctc_fname = sss_path / "test_move_anon_crossTalk_raw_sss.fif" sss_trans_default_fname = sss_path / "test_move_anon_transDefault_raw_sss.fif" @@ -552,13 +559,15 @@ def test_multipolar_bases(): S_tot_fast = _trans_sss_basis( exp, all_coils=_prep_mf_coils(info), trans=info["dev_head_t"] ) - # there are some sign differences for columns (order/degrees) - # in here, likely due to Condon-Shortley. Here we use a - # Magnetometer channel to figure out the flips because the - # gradiometer channels have effectively zero values for first three - # external components (i.e., S_tot[grad_picks, 80:83]) - flips = np.sign(S_tot_fast[2]) != np.sign(S_tot[2]) - flips = 1 - 2 * flips + # The two implementations differ in sign for negative orders only: + # _sss_basis_basic carries the Condon-Shortley (-1)**order there (via + # _sh_negate), whereas _sss_basis follows MaxFilter and builds the + # -order column directly without it. The flips are therefore known + # analytically rather than detected empirically. + orders = np.concatenate( + [_get_degrees_orders(order)[1] for order in (int_order, ext_order)] + ) + flips = np.where(orders < 0, (-1.0) ** orders, 1.0) assert_allclose(S_tot, S_tot_fast * flips, atol=1e-16) @@ -947,24 +956,138 @@ def test_regularization(): raw_fnames = (raw_fname, erm_fname, sample_fname) sss_fnames = (sss_reg_in_fname, sss_erm_reg_in_fname, sss_samp_reg_in_fname) comp_tols = [0, 1, 4] + # These references are from MaxFilter 2.2, whose cut "in_argmax" implements. + # It is exact for the empty-room file but looser for the two head-frame + # cases, since the cut is not the only thing that changed between 2.2 and + # 3.0; test_regularization_mf3 matches 3.0 exactly. + argmax_comp_tols = [3, 0, 8] for ii, rf in enumerate(raw_fnames): raw = read_crop(rf, (0.0, 1.0)) sss_reg_in = read_crop(sss_fnames[ii]) # Test "in" regularization - raw_sss = maxwell_filter(raw, coord_frame=coord_frames[ii], origin=origins[ii]) + with catch_logging() as log: + raw_sss = maxwell_filter( + raw, coord_frame=coord_frames[ii], origin=origins[ii], verbose="debug" + ) assert_meg_snr(raw_sss, sss_reg_in, min_tols[ii], med_tols[ii], msg=rf) # check components match _check_reg_match(raw_sss, sss_reg_in, comp_tols[ii]) + if rf is erm_fname: + # NenonenEtAl2007 report a peak total information of ~480 + # bits/sample for this array. Our snr carries an empirical scale + # factor (see _regularize_in) and this is the only check on it, so + # keep the range loose: the peak also grows with int_order and + # varies with the expansion origin. + peak = float(re.search(r"of peak ([0-9.]+)", log.getvalue()).group(1)) + assert 400 < peak < 560, peak + + # Test "in_argmax" (MaxFilter 2.2-like) regularization + raw_sss = maxwell_filter( + raw, + coord_frame=coord_frames[ii], + origin=origins[ii], + regularize="in_argmax", + ) + _check_reg_match(raw_sss, sss_reg_in, argmax_comp_tols[ii]) + + with pytest.raises(ValueError, match="Invalid value for the 'regularize'"): + maxwell_filter(raw, regularize="both") + with pytest.raises(TypeError, match="must be an instance of str or None"): + maxwell_filter(raw, regularize=1) + + +@pytest.mark.slowtest +@testing.requires_testing_data +@pytest.mark.parametrize( + "origin, fname, nfree", + [ + pytest.param(mf_head_origin, sss_mf3_reg_in_fname, 67, id="std origin"), + pytest.param( + (0.0, 0.02, 0.02), sss_mf3_reg_in_nonstd_fname, 54, id="non-std origin" + ), + ], +) +def test_regularization_mf3(origin, fname, nfree): + """Test regularize="in" against MaxFilter 3.0.""" + # "in" implements 3.0's cut, so unlike the 2.2 references in + # test_regularization these match component-for-component rather than to a + # tolerance. The two MaxFilter versions differ from each other by more than + # we differ from 3.0: at the non-standard origin 2.2 keeps 65, 3.0 keeps 54. + raw = read_crop(raw_fname, (0.0, 1.0)) + sss_mf3 = read_crop(fname) + raw_sss = maxwell_filter(raw, origin=origin, regularize="in") + assert_meg_snr(raw_sss, sss_mf3, 20.0, 100.0, msg=fname.name) + info_py, info_mf3 = _sss_info(raw_sss), _sss_info(sss_mf3) + assert info_mf3["nfree"] == nfree + assert_array_equal(info_py["components"], info_mf3["components"]) + # ...so the SSS rank read back from the file agrees too + assert _get_rank_sss(sss_mf3) == nfree + + +def test_sphere_activation_model(): + """Test the random brain current model against Nenonen et al. 2004.""" + degrees = np.arange(int_order + 1) + a_power, rho_i = _compute_sphere_activation_in(degrees) + # Knuutila's random brain current density estimate, 0.6 µA/m²/√Hz + assert_allclose(rho_i, 0.6e-6, rtol=0.02) + # rho_i is stored precomputed, so make sure it still matches the derivation + # documented beside it (100 fT radial field at r_s=0.13 m, r_in=0.08 m) + ll = np.arange(1, 100) + in_sum = ( + 4 + * np.pi + * np.sum(ll * (ll + 1.0) / (2.0 * ll + 1.0) * (0.08 / 0.13) ** (2 * ll + 2)) + ) + assert_allclose(rho_i, 100e-15 * 1e7 / np.sqrt(in_sum), rtol=1e-6) + # Degree 0 carries no signal; the rest falls off steeply and monotonically + assert a_power[0] == 0 + assert (np.diff(a_power[1:]) < 0).all() + assert a_power[int_order] / a_power[1] < 1e-15 + + +def test_mf_component_order(): + """Test that used moments are encoded in MaxFilter's even-odd order.""" + idx = _mne_ord_to_mf_idx(int_order) + n_in = _get_n_moments(int_order) + assert_array_equal(np.sort(idx), np.arange(n_in)) # a permutation + # NenonenEtAl2007 Fig 2 sweeps L=7..11, i.e. N=63..143 virtual channels + assert_array_equal( + [_get_n_moments(L) for L in range(7, 12)], [63, 80, 99, 120, 143] + ) + degrees, orders = _get_degrees_orders(int_order) + # MaxFilter gives each degree 2 * degree + 1 slots, so these are the flat + # indices of each degree's first moment. + starts = np.cumsum([0] + [2 * d + 1 for d in range(1, int_order)]) + assert_array_equal(starts, [0, 3, 8, 15, 24, 35, 48, 63]) + assert_array_equal(degrees[idx][starts], np.arange(1, int_order + 1)) + assert_array_equal(orders[idx][starts], 0) # every degree starts at m=0 + # and within a degree the orders run 0, +1, -1, +2, -2, ... + assert_array_equal(degrees[idx][24:35], 5) + assert_array_equal(orders[idx][24:35], [0, 1, -1, 2, -2, 3, -3, 4, -4, 5, -5]) + + +def _sss_info(inst): + """Get the SSS record from proc_history. + + MaxFilter 3.0 writes an empty sss_info ahead of the populated one, so the + first proc_history entry is not necessarily the one we want. + """ + for pp in inst.info["proc_history"]: + sss_info = pp.get("max_info", {}).get("sss_info", {}) + if sss_info.get("in_order") is not None: + return sss_info + raise RuntimeError(f"No SSS info found in {inst}") + def _check_reg_match(sss_py, sss_mf, comp_tol): """Check regularization.""" - info_py = sss_py.info["proc_history"][0]["max_info"]["sss_info"] + info_py = _sss_info(sss_py) assert info_py is not None assert len(info_py) > 0 - info_mf = sss_mf.info["proc_history"][0]["max_info"]["sss_info"] + info_mf = _sss_info(sss_mf) n_in = None for inf in (info_py, info_mf): if n_in is None: @@ -975,6 +1098,16 @@ def _check_reg_match(sss_py, sss_mf, comp_tol): assert_allclose( info_py["nfree"], info_mf["nfree"], atol=comp_tol, err_msg=sss_py.filenames[0] ) + # nfree is a sum, so it says nothing about *which* components were dropped. + # Both masks use MaxFilter's order (see _mne_ord_to_mf_idx) and we follow its + # removal sequence, so whichever dropped fewer must be a subset of the other. + drop_mf = set(np.nonzero(np.asarray(info_mf["components"][:n_in]) == 0)[0]) + drop_py = set(np.nonzero(np.asarray(info_py["components"][:n_in]) == 0)[0]) + smaller, larger = sorted((drop_py, drop_mf), key=len) + assert smaller <= larger, ( + f"{sss_py.filenames[0]}: removal sequences diverge, " + f"MNE-only={sorted(drop_py - drop_mf)}, MF-only={sorted(drop_mf - drop_py)}" + ) @pytest.mark.slowtest diff --git a/mne/rank.py b/mne/rank.py index 49708583ebf..65833a5363f 100644 --- a/mne/rank.py +++ b/mne/rank.py @@ -316,19 +316,21 @@ def _get_rank_sss( info = inst if isinstance(inst, Info) else inst.info del inst - proc_info = info.get("proc_history", []) - if len(proc_info) > 1: - logger.info("Found multiple SSS records. Using the first.") - if ( - len(proc_info) == 0 - or "max_info" not in proc_info[0] - or "in_order" not in proc_info[0]["max_info"]["sss_info"] - ): + # Take the first record that actually holds an SSS expansion. MaxFilter 2.2 + # writes it first, but 3.0 emits an empty sss_info ahead of the populated + # one, so we cannot just use proc_history[0]. + proc_info = [ + pp + for pp in info.get("proc_history", []) + if "in_order" in pp.get("max_info", {}).get("sss_info", {}) + ] + if len(proc_info) == 0: raise ValueError( f'Could not find Maxfilter information in info["proc_history"]. {msg}' ) - proc_info = proc_info[0] - max_info = proc_info["max_info"] + if len(proc_info) > 1: + logger.info("Found multiple SSS records. Using the first.") + max_info = proc_info[0]["max_info"] inside = max_info["sss_info"]["in_order"] nfree = (inside + 1) ** 2 - 1 nfree -= ( diff --git a/mne/utils/docs.py b/mne/utils/docs.py index b8377103c6a..ae697fc97ac 100644 --- a/mne/utils/docs.py +++ b/mne/utils/docs.py @@ -3872,9 +3872,21 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75): docdict["regularize_maxwell_reg"] = """ regularize : str | None - Basis regularization type, must be ``"in"`` or None. - ``"in"`` is the same algorithm as the ``-regularize in`` option in - MaxFilter™. + Basis regularization type, must be ``"in"``, ``"in_argmax"``, or None. + Both ``"in"`` options use the same information-theoretic component ordering + as the ``-regularize in`` option in MaxFilter™, and differ only in where + the total-information curve is cut: + + ``"in"`` (default) + Keeps the components giving at least 98% of the peak total information. The curve + can be quite flat, so this errs on the side of including rather than excluding + components. This is the criterion MaxFilter™ 3.0 uses. + ``"in_argmax"`` + Keeps the components at the peak itself, which is what MaxFilter™ 2.2 does. + Use this to match MaxFilter™ 2.2 output more closely; it generally excludes more + components than ``"in"``. + + .. versionadded:: 1.13 """ From 16048f428b71938fc6b37cb30332d28a36550bfa Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 15 Jul 2026 17:51:51 -0400 Subject: [PATCH 2/7] FIX: Whoops --- doc/changes/dev/{14058.bugfix.rst => 14059.bugfix.rst} | 0 doc/changes/dev/{14058.newfeature.rst => 14059.newfeature.rst} | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename doc/changes/dev/{14058.bugfix.rst => 14059.bugfix.rst} (100%) rename doc/changes/dev/{14058.newfeature.rst => 14059.newfeature.rst} (100%) diff --git a/doc/changes/dev/14058.bugfix.rst b/doc/changes/dev/14059.bugfix.rst similarity index 100% rename from doc/changes/dev/14058.bugfix.rst rename to doc/changes/dev/14059.bugfix.rst diff --git a/doc/changes/dev/14058.newfeature.rst b/doc/changes/dev/14059.newfeature.rst similarity index 100% rename from doc/changes/dev/14058.newfeature.rst rename to doc/changes/dev/14059.newfeature.rst From af236ab8770873e8070aed7f22c4e41a40341e69 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 15 Jul 2026 19:02:48 -0400 Subject: [PATCH 3/7] FIX: Read-only --- mne/tests/test_morph.py | 4 ++-- mne/tests/test_source_estimate.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/mne/tests/test_morph.py b/mne/tests/test_morph.py index c363e05ed24..78b8cbed5ae 100644 --- a/mne/tests/test_morph.py +++ b/mne/tests/test_morph.py @@ -263,7 +263,7 @@ def test_surface_source_morph_round_trip(smooth, lower, upper, n_warn, dtype): @testing.requires_testing_data -def test_surface_source_morph_shortcut(): +def test_surface_source_morph_shortcut(subjects_dir_tmp): """Test that our shortcut for smooth=0 works.""" stc = mne.read_source_estimate(fname_smorph) morph_identity = compute_source_morph( @@ -272,7 +272,7 @@ def test_surface_source_morph_shortcut(): "sample", spacing=stc.vertices, smooth=0, - subjects_dir=subjects_dir, + subjects_dir=subjects_dir_tmp, ) stc_back = morph_identity.apply(stc) assert_allclose(stc_back.data, stc.data, rtol=1e-4) diff --git a/mne/tests/test_source_estimate.py b/mne/tests/test_source_estimate.py index 2dfd3bb04d4..6044455586d 100644 --- a/mne/tests/test_source_estimate.py +++ b/mne/tests/test_source_estimate.py @@ -144,7 +144,7 @@ def test_stc_baseline_correction(): @testing.requires_testing_data -def test_spatial_inter_hemi_adjacency(): +def test_spatial_inter_hemi_adjacency(subjects_dir_temp): """Test spatial adjacency between hemispheres.""" # trivial cases conn = spatial_inter_hemi_adjacency(fname_src_3, 5e-6) @@ -169,7 +169,7 @@ def test_spatial_inter_hemi_adjacency(): for hi, hemi in enumerate(("lh", "rh")): has_neighbors = src[hi]["vertno"][np.where(np.any(upper_right, axis=1 - hi))[0]] labels = read_labels_from_annot( - "sample", "aparc.a2009s", hemi, subjects_dir=subjects_dir + "sample", "aparc.a2009s", hemi, subjects_dir=subjects_dir_temp ) use_labels = [ label.name[:-3] From ee50bc9b778e09285eb4e4ba5f195a4ac2f4f261 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 15 Jul 2026 21:14:17 -0400 Subject: [PATCH 4/7] FIX: RW --- mne/tests/test_source_estimate.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/mne/tests/test_source_estimate.py b/mne/tests/test_source_estimate.py index 6044455586d..51f38f432cf 100644 --- a/mne/tests/test_source_estimate.py +++ b/mne/tests/test_source_estimate.py @@ -77,6 +77,7 @@ from mne.source_space._source_space import _get_src_nn from mne.transforms import apply_trans, invert_transform from mne.utils import ( + _chmod_rw_R, _record_warnings, catch_logging, ) @@ -144,7 +145,7 @@ def test_stc_baseline_correction(): @testing.requires_testing_data -def test_spatial_inter_hemi_adjacency(subjects_dir_temp): +def test_spatial_inter_hemi_adjacency(subjects_dir_tmp): """Test spatial adjacency between hemispheres.""" # trivial cases conn = spatial_inter_hemi_adjacency(fname_src_3, 5e-6) @@ -169,7 +170,7 @@ def test_spatial_inter_hemi_adjacency(subjects_dir_temp): for hi, hemi in enumerate(("lh", "rh")): has_neighbors = src[hi]["vertno"][np.where(np.any(upper_right, axis=1 - hi))[0]] labels = read_labels_from_annot( - "sample", "aparc.a2009s", hemi, subjects_dir=subjects_dir_temp + "sample", "aparc.a2009s", hemi, subjects_dir=subjects_dir_tmp ) use_labels = [ label.name[:-3] @@ -1686,6 +1687,7 @@ def test_stc_near_sensors(tmp_path): subjects_dir / "sample" / "surf" / f"{hemi}.white", tmp_path / "sample" / "surf" / f"{hemi}.pial", ) + _chmod_rw_R(tmp_path / "sample" / "surf") # for read-only FSes # here we use a distance is smaller than the inter-sensor distance kwargs = dict( subject="sample", From d88245266d565fbb27ee5c9441769ee5af1378aa Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 15 Jul 2026 21:55:17 -0400 Subject: [PATCH 5/7] FIX: It caught another --- mne/tests/test_source_estimate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mne/tests/test_source_estimate.py b/mne/tests/test_source_estimate.py index 51f38f432cf..1420cd802fa 100644 --- a/mne/tests/test_source_estimate.py +++ b/mne/tests/test_source_estimate.py @@ -1722,7 +1722,7 @@ def test_stc_near_sensors(tmp_path): "sample", smooth=5, spacing=None, - subjects_dir=subjects_dir, + subjects_dir=tmp_path, ).apply(stc_src) lh_idx = np.searchsorted(stc_src_full.vertices[0], stc.vertices[0]) rh_idx = np.searchsorted(stc_src_full.vertices[1], stc.vertices[1]) From 5fe612a674409a1530425722efff42613d4e6894 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Wed, 15 Jul 2026 23:52:37 -0400 Subject: [PATCH 6/7] FIX: Okay --- mne/tests/test_source_estimate.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/mne/tests/test_source_estimate.py b/mne/tests/test_source_estimate.py index 1420cd802fa..fbca8b3daf7 100644 --- a/mne/tests/test_source_estimate.py +++ b/mne/tests/test_source_estimate.py @@ -1715,15 +1715,14 @@ def test_stc_near_sensors(tmp_path): assert src[0]["coord_frame"] == FIFF.FIFFV_COORD_HEAD stc_src = stc_near_sensors(evoked, src=src, **kwargs) assert len(stc_src.data) == 7928 - with pytest.warns(RuntimeWarning, match="not included"): # some removed - stc_src_full = compute_source_morph( - stc_src, - "sample", - "sample", - smooth=5, - spacing=None, - subjects_dir=tmp_path, - ).apply(stc_src) + stc_src_full = compute_source_morph( + stc_src, + "sample", + "sample", + smooth=5, + spacing=None, + subjects_dir=tmp_path, + ).apply(stc_src) lh_idx = np.searchsorted(stc_src_full.vertices[0], stc.vertices[0]) rh_idx = np.searchsorted(stc_src_full.vertices[1], stc.vertices[1]) rh_idx += len(stc_src_full.vertices[0]) From 2a39aedb21c1f69453c5395e8dec065656ce1f6d Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Thu, 16 Jul 2026 08:22:44 -0400 Subject: [PATCH 7/7] FIX: Passing locally --- mne/tests/test_source_estimate.py | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/mne/tests/test_source_estimate.py b/mne/tests/test_source_estimate.py index fbca8b3daf7..1ed48bd8243 100644 --- a/mne/tests/test_source_estimate.py +++ b/mne/tests/test_source_estimate.py @@ -1687,6 +1687,10 @@ def test_stc_near_sensors(tmp_path): subjects_dir / "sample" / "surf" / f"{hemi}.white", tmp_path / "sample" / "surf" / f"{hemi}.pial", ) + copyfile( + subjects_dir / "sample" / "surf" / f"{hemi}.sphere.reg", + tmp_path / "sample" / "surf" / f"{hemi}.sphere.reg", + ) _chmod_rw_R(tmp_path / "sample" / "surf") # for read-only FSes # here we use a distance is smaller than the inter-sensor distance kwargs = dict( @@ -1715,14 +1719,15 @@ def test_stc_near_sensors(tmp_path): assert src[0]["coord_frame"] == FIFF.FIFFV_COORD_HEAD stc_src = stc_near_sensors(evoked, src=src, **kwargs) assert len(stc_src.data) == 7928 - stc_src_full = compute_source_morph( - stc_src, - "sample", - "sample", - smooth=5, - spacing=None, - subjects_dir=tmp_path, - ).apply(stc_src) + with pytest.warns(RuntimeWarning, match="not included"): # some removed + stc_src_full = compute_source_morph( + stc_src, + "sample", + "sample", + smooth=5, + spacing=None, + subjects_dir=tmp_path, + ).apply(stc_src) lh_idx = np.searchsorted(stc_src_full.vertices[0], stc.vertices[0]) rh_idx = np.searchsorted(stc_src_full.vertices[1], stc.vertices[1]) rh_idx += len(stc_src_full.vertices[0])