MAINT: Delegate SSP-SIR reconstruction to MNE

docs/api.rst

@@ -128,8 +128,6 @@ SSP-SIR
    :nosignatures:
 
    mne_denoise.sspsir.SSPSIR
-   mne_denoise.sspsir.compute_sspsir
-   mne_denoise.sspsir.compute_sir
 
 Overcorrection metrics
 ----------------------

docs/changes/devel/90.feature.rst

@@ -0,0 +1,2 @@
+SSP-SIR now delegates Forward-based source-informed projection reconstruction
+to MNE-Python.

examples/sspsir/plot_01_sspsir_basics.py

@@ -12,7 +12,8 @@
 space whose channels no longer correspond to electrodes. SSP-SIR repairs both
 problems: it reconstructs the projected data through a forward model, which
 restores what the head model says must have been there and returns the signal
-to interpretable sensor space.
+to interpretable sensor space. mne-denoise estimates the artifact subspace and
+crossfade, while MNE-Python performs the Forward-based reconstruction.
 
 This example shows the two design choices that matter most in practice: how
 many artifact components to remove, and the fact that the projection is
@@ -34,9 +35,10 @@
 # %%
 # A montage and a forward model
 # -----------------------------
-# As with SOUND, the lead field is built from the montage when no ``forward``
-# is supplied. We build it explicitly here too, so the simulated brain
-# response is something the head model can actually account for.
+# As with SOUND, a spherical Forward is built from the montage when no
+# ``forward`` is supplied. We build its gain explicitly here too, so the
+# simulated brain response is something the head model can actually account
+# for; SSP-SIR delegates its reconstruction to MNE-Python.
 
 montage = mne.channels.make_standard_montage("standard_1020")
 ch_names = [ch for ch in montage.ch_names if ch not in ("A1", "A2")][:32]

mne_denoise/_leadfield.py

@@ -1,31 +1,25 @@
-"""Lead-field construction shared by the forward-model denoisers.
-
-SOUND and SSP-SIR are both forward-model methods: they need a lead-field matrix
-``L`` whose columns are the scalp topographies that cortical current sources
-produce, and they use it as a prior for what a plausible brain signal looks
-like.  Whatever a sensor records that ``L`` cannot explain is treated as noise
-or artifact, so the lead field is what separates "brain" from "everything
-else" in both algorithms.
-
-An individualised forward model computed from the participant's anatomy is
-always the better input.  When none is available, both methods fall back to a
-three-layer spherical head model derived from the electrode montage alone
-(Mutanen et al., 2016, 2018).  That is workable because neither algorithm uses
-``L`` directly — only ``L @ L.T``, the lead-field covariance describing the
-typical cross-correlations between channels — so a head model needs to capture
-those correlations, not the anatomy that produced them.
-
-This module centralises the construction so the two estimators resolve their
-lead field identically:
+"""Forward and lead-field construction shared by the forward-model denoisers.
+
+SOUND consumes an average-referenced lead-field matrix ``L`` whose columns are
+the scalp topographies that cortical current sources produce. SSP-SIR keeps the
+MNE ``Forward`` used to produce those topographies and delegates its
+projection-reconstruction mapping to MNE-Python. Both methods share the same
+spherical ``Forward`` construction when an anatomy-specific model is
+unavailable.
+
+This module centralises the construction so SOUND and SSP-SIR use the same
+forward-model geometry:
 
 - :class:`SphericalHeadModel` — the shell radii and conductivities of the
   fallback head model, with the published values as :data:`REFERENCE_HEAD`.
 - :func:`fibonacci_sphere` — deterministic, quasi-uniform directions on a
   sphere, used to place the source dipoles.
+- :func:`resolve_forward` — align a user-supplied ``Forward`` or build the
+  spherical fallback.
 - :func:`make_spherical_leadfield` — build the fallback lead field from a
   montage.
 - :func:`resolve_leadfield` — choose between a user-supplied forward model and
-  that fallback; this is the entry point the estimators call.
+  that fallback as an ndarray for SOUND.
 
 Every lead field returned here is average referenced, the reference both
 algorithms operate in.
@@ -63,6 +57,7 @@
     "SphericalHeadModel",
     "fibonacci_sphere",
     "make_spherical_leadfield",
+    "resolve_forward",
     "resolve_leadfield",
 ]
 
@@ -173,31 +168,6 @@ def _forward_gain(forward: mne.Forward) -> np.ndarray:
     )
 
 
-def _leadfield_from_forward(forward: mne.Forward, info: mne.Info) -> np.ndarray:
-    """Extract an average-referenced lead field from a user forward solution.
-
-    The forward's rows are reordered to match ``info``'s channel order, so a
-    forward computed elsewhere (with its own channel ordering) lines up with
-    the data being cleaned.
-    """
-    gain = _forward_gain(forward)
-    row_names = list(forward["sol"]["row_names"])
-    if len(row_names) != gain.shape[0]:
-        raise ValueError(
-            "The supplied forward has a different number of row names and "
-            "gain-matrix rows."
-        )
-    wanted = list(info["ch_names"])
-    missing = [ch for ch in wanted if ch not in row_names]
-    if missing:
-        raise ValueError(
-            "The supplied forward model is missing channels present in the data: "
-            f"{missing[:5]}{'...' if len(missing) > 5 else ''}."
-        )
-    idx = [row_names.index(ch) for ch in wanted]
-    return _average_reference(gain[idx])
-
-
 def fibonacci_sphere(n_points: int) -> np.ndarray:
     """Generate unit vectors quasi-uniformly covering the sphere.
 
@@ -237,6 +207,73 @@ def fibonacci_sphere(n_points: int) -> np.ndarray:
     )
 
 
+def _make_spherical_forward(
+    info: mne.Info,
+    *,
+    n_dipoles: int = 5000,
+    head_model: SphericalHeadModel = REFERENCE_HEAD,
+    verbose: bool = False,
+) -> mne.Forward:
+    """Build the spherical fallback ``Forward`` from an EEG montage."""
+    if (
+        isinstance(n_dipoles, (bool, np.bool_))
+        or not isinstance(n_dipoles, Integral)
+        or n_dipoles < 1
+    ):
+        raise ValueError(f"n_dipoles must be a positive integer, got {n_dipoles!r}.")
+    _mne.require_mne("automatic spherical lead-field construction")
+    eeg_picks = _mne.mne.pick_types(info, meg=False, eeg=True, exclude=())
+    if len(eeg_picks) != len(info["ch_names"]):
+        raise ValueError(
+            "Automatic spherical lead-field construction supports EEG channels "
+            "only; provide an explicit forward model for MEG or mixed channel types."
+        )
+
+    with warnings.catch_warnings():
+        # The best-fit sphere centre can sit >20 mm from the head-frame origin
+        # for partial or idealised montages; harmless for spanning the
+        # topography subspace that SOUND and SSP-SIR rely on.
+        warnings.filterwarnings("ignore", message=".*from head frame origin.*")
+        sphere = _mne.mne.make_sphere_model(
+            r0="auto",
+            head_radius="auto",
+            info=info,
+            relative_radii=head_model.relative_radii,
+            sigmas=head_model.conductivities,
+            verbose=verbose,
+        )
+        # ``sphere["layers"][-1]["rad"]`` is the fitted scalp radius in
+        # metres; the source shell tracks it so the geometry stays
+        # proportional.
+        head_radius = float(sphere["layers"][-1]["rad"])
+        directions = fibonacci_sphere(n_dipoles)
+        positions = (
+            directions * (head_model.dipole_relative_radius * head_radius)
+            + sphere["r0"]
+        )
+        src = _mne.mne.setup_volume_source_space(
+            pos={"rr": positions, "nn": directions}, sphere_units="m", verbose=verbose
+        )
+        forward = _mne.mne.make_forward_solution(
+            info,
+            trans=None,
+            src=src,
+            bem=sphere,
+            eeg=True,
+            meg=False,
+            verbose=verbose,
+        )
+        forward = _mne.mne.convert_forward_solution(
+            forward, force_fixed=True, use_cps=False, verbose=verbose
+        )
+    # MNE's forward solver stores gains as float32. Keep the canonical
+    # Forward in double precision so its public gain view and MNE's Gram-based
+    # reconstruction use the same values.
+    forward["sol"]["data"] = np.asarray(forward["sol"]["data"], dtype=float)
+    forward["_orig_sol"] = np.asarray(forward["_orig_sol"], dtype=float)
+    return forward
+
+
 def make_spherical_leadfield(
     info: mne.Info,
     *,
@@ -275,12 +312,12 @@ def make_spherical_leadfield(
 
     Notes
     -----
-    Only ``leadfield @ leadfield.T`` enters SOUND and SSP-SIR, where it acts as
-    the source-covariance prior: it sets the minimum-norm weighting in SOUND
-    and the truncation scale of the source-informed reconstruction in SSP-SIR.
-    Its *spectrum* therefore matters, not merely its span, which is why the
+    SOUND uses ``leadfield @ leadfield.T`` as its source-covariance prior. SSP-SIR
+    retains the Forward and passes its gain, projections, whitening, and
+    source-informed reconstruction to MNE-Python. The *spectrum* of the gain
+    still matters for both methods, not merely its span, which is why the
     published shell-of-radial-dipoles geometry is followed rather than a volume
-    grid.  A volume grid with free orientations spans a comparable subspace
+    grid. A volume grid with free orientations spans a comparable subspace
     (mean principal-angle cosine ~0.94 over the leading topographies for a
     32-channel montage) but has a visibly faster-decaying spectrum, changing
     the effective regularisation.
@@ -291,51 +328,67 @@ def make_spherical_leadfield(
     from such a draw — inside that sampling spread, while being exactly
     repeatable.
     """
-    if (
-        isinstance(n_dipoles, (bool, np.bool_))
-        or not isinstance(n_dipoles, Integral)
-        or n_dipoles < 1
-    ):
-        raise ValueError(f"n_dipoles must be a positive integer, got {n_dipoles!r}.")
-    _mne.require_mne("automatic spherical lead-field construction")
-    eeg_picks = _mne.mne.pick_types(info, meg=False, eeg=True, exclude=())
-    if len(eeg_picks) != len(info["ch_names"]):
-        raise ValueError(
-            "Automatic spherical lead-field construction supports EEG channels "
-            "only; provide an explicit forward model for MEG or mixed channel types."
-        )
+    forward = _make_spherical_forward(
+        info, n_dipoles=n_dipoles, head_model=head_model, verbose=verbose
+    )
+    return _average_reference(_forward_gain(forward))
 
-    with warnings.catch_warnings():
-        # The best-fit sphere centre can sit >20 mm from the head-frame origin
-        # for partial or idealised montages; harmless for spanning the
-        # topography subspace that SOUND and SSP-SIR rely on.
-        warnings.filterwarnings("ignore", message=".*from head frame origin.*")
-        sphere = _mne.mne.make_sphere_model(
-            r0="auto",
-            head_radius="auto",
-            info=info,
-            relative_radii=head_model.relative_radii,
-            sigmas=head_model.conductivities,
-            verbose=verbose,
-        )
-        # ``sphere["layers"][-1]["rad"]`` is the fitted scalp radius in metres;
-        # the source shell tracks it so the geometry stays proportional.
-        head_radius = float(sphere["layers"][-1]["rad"])
-        directions = fibonacci_sphere(n_dipoles)
-        positions = (
-            directions * (head_model.dipole_relative_radius * head_radius)
-            + sphere["r0"]
-        )
-        src = _mne.mne.setup_volume_source_space(
-            pos={"rr": positions, "nn": directions}, sphere_units="m", verbose=verbose
-        )
-        fwd = _mne.mne.make_forward_solution(
-            info, trans=None, src=src, bem=sphere, eeg=True, meg=False, verbose=verbose
+
+def resolve_forward(
+    *,
+    inst: mne.io.BaseRaw | mne.BaseEpochs | mne.Evoked | None,
+    ch_names: list[str] | None,
+    n_channels: int,
+    method: str,
+    forward: mne.Forward | None = None,
+    n_dipoles: int = 5000,
+    head_model: SphericalHeadModel = REFERENCE_HEAD,
+) -> mne.Forward:
+    """Resolve and align the ``Forward`` used by a forward-model denoiser.
+
+    MNE-Python performs the channel-name alignment for explicit forwards. For
+    array input, the forward's own row order is the channel contract because
+    the array carries no names.
+    """
+    _mne.require_mne(f"{method} Forward resolution")
+    if inst is not None:
+        info = inst.copy().pick(ch_names).info
+        if forward is None:
+            return _make_spherical_forward(
+                info, n_dipoles=n_dipoles, head_model=head_model
+            )
+        _forward_gain(forward)
+        resolved = _mne.mne.pick_channels_forward(
+            forward, include=info["ch_names"], ordered=True
         )
-        fwd = _mne.mne.convert_forward_solution(
-            fwd, force_fixed=True, use_cps=False, verbose=verbose
+        resolved["sol"]["data"] = np.asarray(resolved["sol"]["data"], dtype=float)
+        resolved["_orig_sol"] = np.asarray(resolved["_orig_sol"], dtype=float)
+        return resolved
+
+    if forward is not None:
+        gain = _forward_gain(forward)
+        if gain.shape[0] != n_channels:
+            raise ValueError(
+                "For array input, the forward must have the same number of "
+                f"channels as the data ({gain.shape[0]} vs {n_channels})."
+            )
+        row_names = list(forward["sol"]["row_names"])
+        if len(row_names) != gain.shape[0]:
+            raise ValueError(
+                "The supplied forward has a different number of row names and "
+                "gain-matrix rows."
+            )
+        resolved = _mne.mne.pick_channels_forward(
+            forward, include=row_names, ordered=True
         )
-    return _average_reference(np.asarray(fwd["sol"]["data"], dtype=float))
+        resolved["sol"]["data"] = np.asarray(resolved["sol"]["data"], dtype=float)
+        resolved["_orig_sol"] = np.asarray(resolved["_orig_sol"], dtype=float)
+        return resolved
+
+    raise ValueError(
+        f"{method} needs channel positions: pass an MNE object with a montage, "
+        "or provide a `forward` for array input."
+    )
 
 
 def resolve_leadfield(
@@ -350,9 +403,9 @@ def resolve_leadfield(
 ) -> np.ndarray:
     """Resolve the lead field an estimator should use.
 
-    Both :class:`~mne_denoise.sound.SOUND` and
-    :class:`~mne_denoise.sspsir.SSPSIR` accept the same three kinds of input,
-    and this function is where that choice is made once for both:
+    SOUND's ndarray lead-field input can be resolved from the same forward
+    model used by SSP-SIR. This function handles the three supported input
+    cases:
 
     - **MNE object, no forward** — build the spherical fallback from its
       montage.
@@ -407,23 +460,13 @@ def resolve_leadfield(
     --------
     make_spherical_leadfield : The fallback this dispatches to.
     """
-    if inst is not None:
-        _mne.require_mne("MNE lead-field resolution")
-        info = inst.copy().pick(ch_names).info
-        if forward is not None:
-            return _leadfield_from_forward(forward, info)
-        return make_spherical_leadfield(
-            info, n_dipoles=n_dipoles, head_model=head_model
-        )
-    if forward is not None:
-        gain = _forward_gain(forward)
-        if gain.shape[0] != n_channels:
-            raise ValueError(
-                "For array input, the forward must have the same number of "
-                f"channels as the data ({gain.shape[0]} vs {n_channels})."
-            )
-        return _average_reference(gain)
-    raise ValueError(
-        f"{method} needs channel positions: pass an MNE object with a montage, "
-        "or provide a `forward` for array input."
+    forward = resolve_forward(
+        inst=inst,
+        ch_names=ch_names,
+        n_channels=n_channels,
+        method=method,
+        forward=forward,
+        n_dipoles=n_dipoles,
+        head_model=head_model,
     )
+    return _average_reference(_forward_gain(forward))