"""SSP-SIR source-informed reconstruction."""
from __future__ import annotations
import warnings
from numbers import Integral, Real
import numpy as np
from scipy.ndimage import uniform_filter1d
from scipy.signal import butter, filtfilt
from scipy.special import expit
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.utils.validation import check_is_fitted
from . import _mne
from ._data import extract_data_from_mne, reconstruct_mne_object
from ._leadfield import (
_average_reference,
_make_spherical_forward,
_validate_leadfield,
)
from ._logging import logger, verbose
from ._validation import (
check_channel_layout,
check_matching_sfreq,
check_option,
check_positive_real,
)
__all__ = ["SSPSIR"]
#: 10-90% transition width (s) of the crossfade around a user artifact window.
_SMOOTH_LENGTH = 0.010
def _artifact_subspace(
svd_input: np.ndarray, n_components
) -> tuple[np.ndarray, int, np.ndarray]:
"""Estimate the high-frequency artifact subspace and component count."""
svd_input = np.asarray(svd_input, dtype=float)
if svd_input.ndim != 2 or 0 in svd_input.shape:
raise ValueError(
f"svd_input must be a non-empty 2D array, got shape {svd_input.shape}."
)
if not np.isfinite(svd_input).all():
raise ValueError("svd_input must contain only finite values.")
if isinstance(n_components, (bool, np.bool_)) or not isinstance(n_components, Real):
raise ValueError(
"n_components must be a positive integer or a variance fraction in (0, 1)."
)
is_count = isinstance(n_components, Integral)
if is_count:
n_pc = int(n_components)
if n_pc < 1:
raise ValueError(
f"n_components must be a positive integer, got {n_components}."
)
elif not 0.0 < float(n_components) < 1.0:
raise ValueError(
"A floating-point n_components must be a variance fraction in "
f"(0, 1), got {n_components}."
)
u, s, _ = np.linalg.svd(svd_input, full_matrices=False)
if s[0] == 0.0:
raise ValueError("Cannot estimate an artifact subspace from all-zero data.")
if is_count:
if n_pc > u.shape[1]:
raise ValueError(
f"n_components={n_pc} exceeds the {u.shape[1]} available "
"artifact dimensions."
)
else:
power = (s / s[0]) ** 2
n_pc = int(
np.searchsorted(np.cumsum(power), float(n_components) * power.sum()) + 1
)
if n_pc >= svd_input.shape[0]:
raise ValueError(
"n_components must leave at least one channel dimension outside "
"the artifact subspace."
)
return u[:, :n_pc], n_pc, s
def _as_mne_projection(topography: np.ndarray, ch_names: list[str], desc: str):
"""Wrap one channel topography as an inactive MNE projection."""
_mne.require_mne("SSP-SIR projections")
topography = np.asarray(topography, dtype=float)
if topography.ndim == 1:
topography = topography[:, np.newaxis]
return _mne.mne.Projection(
data={
"data": topography.T,
"col_names": list(ch_names),
"row_names": None,
"ncol": len(ch_names),
"nrow": 1,
},
active=False,
desc=desc,
kind=1, # FIFFV_PROJ_ITEM_FIELD
explained_var=None,
)
def _as_mne_projections(topographies: np.ndarray, ch_names: list[str]) -> list:
"""Wrap artifact topographies as MNE projection objects."""
return [
_as_mne_projection(
topographies[:, i], ch_names, desc=f"SSP-SIR artifact {i + 1}"
)
for i in range(topographies.shape[1])
]
[docs]
class SSPSIR(BaseEstimator, TransformerMixin):
"""Source-informed signal-space projection for TMS-evoked muscle artifact removal.
SSP-SIR projects out an artifact subspace, reconstructs through a lead field,
and blends the projected and unprojected reconstructions around the artifact
window.
Parameters
----------
n_components : int or float
Number of artifact components, or high-frequency variance fraction in (0, 1).
forward : mne.Forward or None, default=None
Optional explicit forward solution. For compatible EEG MNE input with a
montage, None uses a spherical fallback; MEG or mixed-channel MNE input
and NumPy input require an explicit forward.
art_window : tuple of float or None, default=None
Artifact interval (tmin, tmax) in seconds; None derives a high-frequency
envelope from the epoch.
blend : {"auto", "constant"}, default="auto"
Crossfade rule.
high_pass : float, default=100.0
High-pass cutoff in Hz for artifact-subspace estimation.
M : int or None, default=None
Source-informed reconstruction rank; None uses data rank minus artifact rank.
smooth_length : float, default=0.010
Crossfade transition width in seconds.
sfreq : float or None, default=None
Sampling frequency for NumPy input.
n_dipoles : int, default=5000
Dipoles for a generated spherical EEG lead field.
verbose : bool, str, int, or None, default=None
Logging level.
Attributes
----------
leadfield_ : ndarray
Lead field used during fitting.
artifact_topographies_ : ndarray
Fitted artifact subspace.
operator_ : ndarray
Projected reconstruction operator.
operator_orig_ : ndarray
Unprojected reconstruction operator.
kernel_ : ndarray
Temporal crossfade weights.
n_components_ : int
Effective artifact-component count.
M_ : int
Effective reconstruction rank.
singular_values_ : ndarray
Singular values used for artifact-subspace selection.
projs_ : list of mne.Projection
Artifact directions when fitted on named MNE data.
See Also
--------
mne_denoise.sound.SOUND
Another forward-model-based denoising method with a different noise model.
Notes
-----
NumPy input uses (n_channels, n_times) or (n_epochs, n_channels, n_times).
MNE Raw, Epochs, and Evoked inputs are supported and returned without mutation.
The artifact-window reconstruction is time-locked to the fitted data
:footcite:p:`mutanen2016_sspsir,mutanen2022_source_artifact,mutanen2024_sspsir_simulation,hernandez_pavon2022_tms_review`.
References
----------
.. footbibliography::
Examples
--------
A preloaded MNE ``Epochs`` object with a compatible EEG montage and a
sampling rate whose Nyquist frequency is above the configured high-pass
cutoff can use the spherical fallback:
.. code-block:: python
from mne_denoise.sspsir import SSPSIR
model = SSPSIR(
n_components=3,
art_window=(0.005, 0.050),
)
clean = model.fit_transform(epochs)
"""
def __init__(
self,
*,
n_components=None,
forward=None,
art_window=None,
blend: str = "auto",
high_pass: float = 100.0,
M=None,
smooth_length: float = _SMOOTH_LENGTH,
sfreq=None,
n_dipoles: int = 5000,
verbose: bool | str | int | None = None,
):
self.n_components = n_components
self.forward = forward
self.art_window = art_window
self.blend = blend
self.high_pass = high_pass
self.M = M
self.smooth_length = smooth_length
self.sfreq = sfreq
self.n_dipoles = n_dipoles
self.verbose = verbose
def _resolve_sfreq(self, sfreq):
sfreq = sfreq if sfreq is not None else self.sfreq
if sfreq is None:
raise ValueError(
"SSP-SIR needs a sampling frequency: pass an MNE object or set sfreq."
)
if isinstance(sfreq, (bool, np.bool_)) or not isinstance(sfreq, Real):
raise ValueError(f"sfreq must be a positive finite number, got {sfreq!r}.")
sfreq = float(sfreq)
if not np.isfinite(sfreq) or sfreq <= 0.0:
raise ValueError(f"sfreq must be a positive finite number, got {sfreq!r}.")
if not 0.0 < float(self.high_pass) < sfreq / 2.0:
raise ValueError(
f"high_pass ({self.high_pass} Hz) must be between 0 and the "
f"Nyquist frequency ({sfreq / 2.0} Hz) for sfreq={sfreq} Hz."
)
return sfreq
def _svd_input(self, evoked, sfreq, times):
"""Filter the evoked data and construct the artifact-subspace input and crossfade."""
b, a = butter(2, self.high_pass / (sfreq / 2.0), btype="high")
data_high = filtfilt(b, a, evoked, axis=1)
if self.art_window is not None:
tmin, tmax = self.art_window
mask = (times >= tmin) & (times <= tmax)
if not mask.any():
raise ValueError("art_window does not overlap the data time range.")
# Smooth step function around the artifact window.
slope = 4.0 / float(self.smooth_length)
kernel = expit(slope * (times - tmin)) - expit(slope * (times - tmax))
return data_high[:, mask], kernel
# Weight by a 50 ms sliding RMS of high-frequency power.
win = max(1, int(round(sfreq / 1000.0 * 50.0)))
power = uniform_filter1d(data_high**2, size=win, axis=1, mode="nearest")
kernel = power.mean(axis=0)
kernel = np.sqrt(kernel / kernel.max()) if kernel.max() > 0 else kernel
return kernel[None, :] * data_high, kernel
[docs]
@verbose
def fit(
self,
X,
y=None,
*,
verbose: bool | str | int | None = None,
):
"""Fit the SSP-SIR artifact subspace and reconstruction operators.
Parameters
----------
X : ndarray, Raw, Epochs, or Evoked
EEG data used for fitting. NumPy input is channel-first.
y : None, default=None
Ignored for scikit-learn compatibility.
verbose : bool, str, int, or None, default=None
Logging level.
Returns
-------
SSPSIR
The fitted estimator.
"""
if self.n_components is None:
raise ValueError(
"n_components must be set (number of artifact PCs to remove, or a "
"variance fraction in (0, 1))."
)
check_option(self.blend, name="blend", allowed=("auto", "constant"))
check_positive_real(self.high_pass, name="high_pass")
check_positive_real(self.smooth_length, name="smooth_length")
if isinstance(self.n_dipoles, (bool, np.bool_)) or not isinstance(
self.n_dipoles, Integral
):
raise ValueError(
f"n_dipoles must be a positive integer, got {self.n_dipoles!r}."
)
if self.n_dipoles < 1:
raise ValueError(
f"n_dipoles must be a positive integer, got {self.n_dipoles!r}."
)
if self.M is not None and (
isinstance(self.M, (bool, np.bool_))
or not isinstance(self.M, Integral)
or self.M < 1
):
raise ValueError(f"M must be a positive integer or None, got {self.M!r}.")
if self.art_window is not None:
if not isinstance(self.art_window, tuple) or len(self.art_window) != 2:
raise ValueError("art_window must be a (tmin, tmax) tuple or None.")
tmin, tmax = self.art_window
if any(
isinstance(value, (bool, np.bool_))
or not isinstance(value, Real)
or not np.isfinite(value)
for value in (tmin, tmax)
):
raise ValueError("art_window values must be finite numbers.")
if tmin >= tmax:
raise ValueError(
f"art_window must satisfy tmin < tmax, got {self.art_window}."
)
data, sfreq, _, orig_inst, _, ch_names = extract_data_from_mne(X)
sfreq = self._resolve_sfreq(sfreq)
times = getattr(orig_inst, "times", None)
# Average reference and an evoked (trial-averaged) view for the subspace.
evoked = data.mean(axis=0) if data.ndim == 3 else np.asarray(data, float)
evoked = evoked - evoked.mean(axis=0, keepdims=True)
n_channels = evoked.shape[0]
if times is None:
times = np.arange(evoked.shape[1], dtype=float) / sfreq
else:
times = np.asarray(times, dtype=float)
svd_input, kernel = self._svd_input(evoked, sfreq, times)
(
self.artifact_topographies_,
self.n_components_,
self.singular_values_,
) = _artifact_subspace(svd_input, self.n_components)
_mne.require_mne("SSP-SIR Forward resolution")
if orig_inst is not None:
fitted_ch_names = list(ch_names)
if self.forward is None:
info = orig_inst.copy().pick(fitted_ch_names).info
resolved_forward = _make_spherical_forward(
info,
n_dipoles=self.n_dipoles,
)
else:
_validate_leadfield(
self.forward["sol"]["data"],
what="The supplied forward gain matrix",
)
resolved_forward = _mne.mne.pick_channels_forward(
self.forward,
include=fitted_ch_names,
ordered=True,
copy=True,
)
else:
if self.forward is None:
raise ValueError(
"SSP-SIR needs channel positions: pass an MNE object with a "
"montage, or provide a `forward` for array input."
)
gain = _validate_leadfield(
self.forward["sol"]["data"],
what="The supplied forward gain matrix",
)
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(self.forward["sol"]["row_names"])
if len(row_names) != n_channels:
raise ValueError(
"The supplied forward has a different number of row names and "
"gain-matrix rows."
)
fitted_ch_names = row_names
resolved_forward = _mne.mne.pick_channels_forward(
self.forward,
include=fitted_ch_names,
ordered=True,
copy=True,
)
gain = _validate_leadfield(
resolved_forward["sol"]["data"],
what="The supplied forward gain matrix",
)
self.leadfield_ = _average_reference(gain)
if orig_inst is not None:
source_info = orig_inst.info
else:
source_info = resolved_forward["info"]
source_picks = _mne.mne.pick_channels(
source_info["ch_names"], include=fitted_ch_names, ordered=True
)
fitted_ch_types = source_info.get_channel_types(picks=source_picks)
clean_info = _mne.mne.create_info(fitted_ch_names, sfreq, fitted_ch_types)
artifact_projs = _as_mne_projections(
self.artifact_topographies_, fitted_ch_names
)
self.projs_ = artifact_projs if orig_inst is not None else []
channel_types = set(fitted_ch_types)
has_eeg = "eeg" in channel_types
has_meg = bool(channel_types.intersection(("mag", "grad", "meg")))
if has_eeg and has_meg:
raise ValueError(
"SSP-SIR requires one homogeneous sensor family for a single M; "
"mixed EEG and MEG Forward rows are not supported."
)
if has_eeg:
rank_key = "eeg"
common_desc = "Average EEG reference"
elif has_meg:
rank_key = "meg"
common_desc = "SSP-SIR common mode"
else:
raise ValueError("SSP-SIR requires EEG or MEG data channels.")
common_mode = np.ones(n_channels) / np.sqrt(n_channels)
common_mode_proj = _as_mne_projection(common_mode, fitted_ch_names, common_desc)
data_rank = int(np.linalg.matrix_rank(evoked))
requested_M = (
int(self.M) if self.M is not None else data_rank - self.n_components_
)
if requested_M < 1:
raise ValueError(f"M must be a positive integer, got {requested_M!r}.")
projected_leadfield = self.leadfield_ - self.artifact_topographies_ @ (
self.artifact_topographies_.T @ self.leadfield_
)
available_rank = int(np.linalg.matrix_rank(projected_leadfield))
if available_rank == 0:
raise ValueError(
"Cannot reconstruct from the projected lead field: its numerical "
"rank is zero."
)
M = requested_M
if available_rank < M:
warnings.warn(
f"M={M} exceeds the numerical rank ({available_rank}) of the "
"projected lead field; using "
f"M={available_rank} instead.",
RuntimeWarning,
stacklevel=2,
)
M = available_rank
identity = _mne.mne.EvokedArray(
np.eye(n_channels),
clean_info,
tmin=0.0,
verbose=False,
)
all_projs = [common_mode_proj, *artifact_projs]
identity.add_proj(all_projs, verbose=False)
projected = identity.copy()
self.operator_ = projected.reconstruct_proj(
projs=projected.info["projs"],
forward=resolved_forward,
rank={rank_key: M},
verbose=False,
).data
self.M_ = int(np.linalg.matrix_rank(self.operator_))
unprojected = identity.copy()
self.operator_orig_ = unprojected.reconstruct_proj(
projs=[unprojected.info["projs"][0]],
forward=resolved_forward,
rank={rank_key: self.M_},
verbose=False,
).data
self.kernel_ = (
np.ones(evoked.shape[1])
if self.blend == "constant"
else np.clip(kernel, 0.0, 1.0)
)
self.sfreq_ = sfreq
self.times_ = times.copy()
self._mne_ch_names_ = fitted_ch_names
logger.info(
"SSP-SIR: channels=%d, removed %d artifact component(s), "
"SIR truncation M=%d (data rank %d), blend=%s",
n_channels,
self.n_components_,
self.M_,
data_rank,
self.blend,
)
return self
