diff --git a/CHANGELOG.md b/CHANGELOG.md index 01770f1d..036f9503 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,24 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added + +- **ICanClean pseudo-reference mode** (restores functionality removed in + `e38e8f5` without a changelog entry; see issue #68) + - `pseudo_ref=True` derives the CCA reference block from the primary + channels themselves rather than from physical noise electrodes, for + recordings with no dual-layer cap. Implements the pseudo-reference + method of Downey & Ferris 2023, *Sensors* 23(19):8214. + - `filter_ref=(btype, freqs)` shapes the reference block before CCA, + using scipy's own filter-kind names: `'bandstop'`, `'bandpass'`, + `'highpass'`, `'lowpass'`; zero-phase 4th-order Butterworth. Usable + on its own to filter physical reference channels. + - `ref_channels` must be left as `None` when `pseudo_ref=True`; the two + are mutually exclusive. + - `pseudo_ref=True` without `filter_ref` raises: an unfiltered copy of + the primary block is perfectly correlated with itself and would remove + the entire signal. + ## [0.0.1] - 2026-01-23 ### Added diff --git a/docs/changes/devel/26.feature.rst b/docs/changes/devel/26.feature.rst index 24e9c9f2..cbef54cc 100644 --- a/docs/changes/devel/26.feature.rst +++ b/docs/changes/devel/26.feature.rst @@ -8,7 +8,9 @@ - Implemented four specialized cleaning modes: * ``'global'``: Static CCA cleaning over the entire recording. * ``'sliding'``: Dynamic, windowed CCA cleaning for non-stationary artifacts. - * ``'calibrated'``: Reference-driven cleaning using a dedicated calibration period. + * ``'calibrated'``: One global CCA decomposition, reused for window-local + scoring and subtraction. It calibrates on the same data it cleans; there + is no separate calibration recording. * ``'hybrid'``: A multi-stage approach combining global subspace removal with sliding-window refinement. - Added support for flexible cleaning bases (primary, reference, or both) and diff --git a/docs/changes/devel/76.bugfix.rst b/docs/changes/devel/76.bugfix.rst new file mode 100644 index 00000000..928daaa8 --- /dev/null +++ b/docs/changes/devel/76.bugfix.rst @@ -0,0 +1,13 @@ +**iCanClean**: +- Fixed ``threshold`` and ``global_threshold`` accepting values that silently + broke the estimator: a value above 1 made it a pass-through (no :math:`R^2` + could exceed it), a value below 0 flagged every component, and a numeric + string such as ``"0.5"`` was compared against floats without conversion. + Both parameters are now range-checked to ``[0, 1]``. +- Fixed ``_reset_qc_attrs`` not clearing ``global_n_windows_`` or + ``sliding_n_windows_``, so stale hybrid window counts survived a re-fit in a + different mode while every sibling ``global_*``/``sliding_*`` attribute was + correctly cleared. +- Fixed the ``'calibrated'`` mode docstring, which described "a dedicated + calibration period" that does not exist -- it calibrates on the same data it + cleans. diff --git a/docs/changes/devel/76.doc.rst b/docs/changes/devel/76.doc.rst new file mode 100644 index 00000000..c1c0ab39 --- /dev/null +++ b/docs/changes/devel/76.doc.rst @@ -0,0 +1,3 @@ +Added ``docs/icanclean.md``, the module's first narrative documentation page, +covering threshold scale, window conditioning, the four operating modes and +the two reference constructions. diff --git a/docs/changes/devel/76.feature.rst b/docs/changes/devel/76.feature.rst new file mode 100644 index 00000000..7715a5da --- /dev/null +++ b/docs/changes/devel/76.feature.rst @@ -0,0 +1,32 @@ +Added ``threshold='null'`` to :class:`mne_denoise.icanclean.ICanClean`, which +sets the rejection threshold from the data instead of a constant. It estimates +the largest squared canonical correlation attributable to sampling noise for +the current window length and channel counts -- by recomputing the spectrum +against circularly shifted copies of the reference -- and rejects only +components exceeding it. + +``threshold`` is an absolute :math:`R^2`, but the achievable scale is set by +the data and by how the reference is built. It varies by ~500x between a +physical dual-layer reference (median :math:`R^2` 0.001) and a pseudo-reference +(median 0.64), and by 3.4x between subjects of a single cohort. A constant +tuned on one recording can silently become a no-op or remove nearly everything +on the next. + +``'null'`` also closes a failure that is invisible today: canonical +correlations are upward-biased when a window is short relative to +``n_primary + n_reference``. On independent data with 40 + 40 channels and a +sample-to-variable ratio of 2.0 -- a 2 s window at 250 Hz on a 120 + 120 +montage -- a fixed ``threshold=0.85`` removes about 10 of 40 components. +``'null'`` removes 0.12, while still recovering exactly 0, 1, 3 and 8 genuinely +shared components when those are injected. + +Use ``null_random_state`` for reproducible surrogates. ``'null'`` determines +whether a component shares *real* variance with the reference, not whether +that variance is artifact. + +Also added the fitted attributes ``max_r2_``, ``thresholds_`` and +``samples_per_variable_``, recording per window the highest squared canonical +correlation observed, the threshold applied, and the sample-to-variable ratio. +Without these, ``n_removed_ == 0`` is indistinguishable from "the threshold was +above every achievable :math:`R^2`" -- a distinction that decides whether a +null result is about the data or about the configuration. diff --git a/docs/icanclean.md b/docs/icanclean.md new file mode 100644 index 00000000..19919051 --- /dev/null +++ b/docs/icanclean.md @@ -0,0 +1,191 @@ +# iCanClean + +`ICanClean` removes latent artifact subspaces that the primary channels share +with a reference block, using canonical correlation analysis. A component is +rejected when its squared canonical correlation exceeds a threshold: + +```python +bad_mask = r2 >= threshold +``` + +That one line is the whole selection rule, and choosing `threshold` well is the +only thing standing between a useful cleaner and a silent no-op. + +## The threshold is not a transferable number + +`threshold` is an **absolute** squared canonical correlation, but the scale it +lives on is set by the data and by how the reference block is built. Measured on +one recording (ds004505, 120 scalp and 120 outward-facing noise channels): + +| reference construction | max r² | median r² | components ≥ 0.25 | +|---|---|---|---| +| dual-layer (physical noise electrodes) | **0.024** | 0.0013 | 0 of 120 | +| pseudo-reference (band-stopped copy of the EEG) | **0.825** | 0.642 | 119 of 120 | + +The same threshold means entirely different things in the two cases — the medians +differ by roughly 500×. A grid of `0.25 … 0.95` removes **exactly zero** +components in the dual-layer case, because the lowest value tested is already ten +times the highest achievable correlation. The same grid removes almost everything +in the pseudo-reference case. + +The scale also moves *within* a single homogeneous cohort: across five subjects of +one dataset the achievable maximum ranged 0.023 → 0.078, a 3.4× spread. + +**Practical consequence:** a threshold taken from a paper, or tuned on one +recording, should not be assumed to transfer. Check what your data can actually +reach before choosing one: + +```python +from mne_denoise.icanclean import ICanClean + +icc = ICanClean(sfreq=raw.info["sfreq"], ref_channels=noise_ch, threshold=0.7) +icc.fit_transform(raw) +print(icc.max_r2_) # highest r2 actually observed, per window +print(icc.thresholds_) # the threshold applied, per window +``` + +If `max_r2_` sits below `thresholds_`, the estimator was a pass-through and +`n_removed_` will be zero — not because the data was clean, but because the +threshold was unreachable. + +## `threshold='null'` — let the data set the scale + +Rather than supplying a constant, ask what r² *sampling noise alone* would +produce for the current window length and channel counts, and reject only what +exceeds it: + +```python +icc = ICanClean( + sfreq=raw.info["sfreq"], + ref_channels=noise_ch, + threshold="null", # calibrated per window + null_random_state=0, # reproducible surrogates +) +``` + +The null is built by circularly shifting the reference block and recomputing the +spectrum, which preserves each channel's autocorrelation and power spectrum while +destroying cross-block alignment. (A plain sample permutation would destroy the +autocorrelation too, giving surrogates that cannot reach the correlations real +data reaches — an anticonservative null.) + +This matters most where a fixed threshold fails silently. Canonical correlations +are upward-biased when a window is short relative to `n_primary + n_reference`; as +that ratio approaches 1, every correlation approaches 1 whether or not anything is +shared. Components falsely removed on independent data, 40 primary and 40 +reference channels: + +| n / (p + q) | `'null'` threshold | `'null'` | fixed 0.65 | fixed 0.85 | +|---|---|---|---|---| +| 1.2 | 0.999 | **0.20** | 20.4 | 14.4 | +| 2.0 | 0.992 | **0.12** | 16.8 | 10.6 | +| 10.0 | 0.904 | **0.00** | 11.4 | 2.6 | +| 300.0 | 0.124 | **0.04** | 0.0 | 0.0 | + +The threshold self-adapts from 0.999 to 0.124 with no input. Note the second row: +a 2 s window at 250 Hz on a 120 + 120 montage gives a ratio of 2.08, and a fixed +`threshold=0.85` there removes about a quarter of the components from data with +nothing in common. + +Safety does not come from timidity — with 0, 1, 3 and 8 genuinely shared +components injected, `'null'` recovers exactly 0, 1, 3 and 8. + +**What it does not do.** `'null'` decides whether a component shares *real* +variance with the reference. It does not decide whether that variance is +*artifact*. With a pseudo-reference — a band-stopped copy of the primary block — +almost every component shares real variance, so `'null'` will select broadly +there. It solves the degeneracy problem, not the selectivity problem. + +## Window length and conditioning + +`segment_len` sets the window that is cleaned. `stats_segment_len`, when larger, +sets a wider window that the CCA is *estimated* on, while only the inner +`segment_len` is written back. This is how you get short, responsive correction +windows backed by statistically adequate estimates: + +```python +icc = ICanClean( + sfreq=250.0, + ref_channels=noise_ch, + mode="sliding", + segment_len=2.0, # corrected span + stats_segment_len=32.0, # estimation span + threshold="null", +) +``` + +As a rule of thumb, keep `n_samples` at least ten times `n_primary + n_reference` +in whichever window the CCA is estimated on. `samples_per_variable_` reports the +achieved ratio. + +## Operating modes + +| mode | CCA decompositions | what it does | +|---|---|---| +| `'global'` | 1 | one decomposition on the whole recording, subtracted once | +| `'sliding'` | one per window | fresh decomposition per window, overlap-added | +| `'calibrated'` | 1 | one global decomposition, reused for window-local scoring | +| `'hybrid'` | 1 + one per window | a global pass, then a sliding pass on its output | + +All four are **batch**: each sees the whole recording before returning cleaned +data, and none updates its decomposition sample by sample. + +Which is best is **artifact-dependent**, not universal. On a stationary artifact a +long window exploits more data; on a non-stationary one a short window tracks the +change. Measured on real recordings, best operating point at 90% alpha retention: + +| artifact | `'global'` | `'sliding'` 4 s | +|---|---|---| +| ocular (blink, n = 40) | **89.8%** removed | 61.3% | +| cardiac (n = 8) | 37.8% | **41.6%** | + +`'calibrated'` calibrates on the same data it cleans — there is no separate +calibration recording — and costs roughly one decomposition rather than one per +window. Note its `correlations_` are window-local Pearson correlations squared, +not canonical correlations squared, and are **not** sorted descending, so row *k* +is not comparable with the other modes. + +`'hybrid'` is an mne-denoise extension, not part of the published algorithm; the +reference implementation applies iCanClean exactly once. It is motivated by the +authors' open question about whether "incorporating larger windows of data" +helps — pass 1 estimates on the whole recording, pass 2 then tracks what is left. +Current evidence is suggestive but underpowered: across 40 subjects it removed +more artifact than the best single-pass mode at every preservation floor, but +every confidence interval included zero. + +```{note} +`'hybrid'` is a **two-pass batch** refinement. It is not an online or recursive +estimator: both passes see the entire recording, and neither updates its +decomposition as samples arrive. + +The 2023 paper raises two distinct possibilities in one sentence — using larger +windows, and computing CCA *recursively*. Only the first is addressed here. A +recursive formulation, in the sense of moments updated incrementally for causal +streaming use, is a different design on a different axis, and nothing on this +page should be read as evaluating one. +``` + +## Reference modes + +The reference block can be physical noise electrodes (`ref_channels`), or derived +from the EEG itself with `pseudo_ref=True` plus a `filter_ref` band-stop that +keeps out-of-band drift and EMG while removing the brain band. + +**iCanClean performs as well as its reference observes the artifact you are +scoring.** A mechanically-coupled outward-facing noise layer tracks head *motion*; +it is largely blind above ~12 Hz. Scoring it against broadband muscle asks it to +remove something it cannot see. On treadmill walking, where the artifact is +gait-locked motion, dual-layer raises the good-component count by ~23%; on a +whole-body sport scored against neck EMG, the same configuration removes nothing. + +## References + +- Downey & Ferris (2022), *The iCanClean Algorithm*, arXiv:2201.11798 +- Downey & Ferris (2023), *Sensors* **23**(19):8214 +- Gonsisko, Ferris & Downey (2023), *Sensors* **23**(2):928 + +```{note} +A public U.S. patent application has been filed for the iCanClean method +(US20230363718A1). Patent applications, and any resulting patents, may affect +commercial use. +``` diff --git a/docs/index.rst b/docs/index.rst index 1dd968c9..48fabec0 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -17,6 +17,7 @@ to extract reproducible or rhythmic components while preserving data rank. bss_cca sns ssa + icanclean auto_examples/index .. toctree:: diff --git a/mne_denoise/icanclean/core.py b/mne_denoise/icanclean/core.py index 03a925a2..b6503940 100644 --- a/mne_denoise/icanclean/core.py +++ b/mne_denoise/icanclean/core.py @@ -72,6 +72,139 @@ logger = logging.getLogger(__name__) +#: Default number of circular-shift surrogates for ``threshold='null'``. 20 is +#: the floor at which the default alpha's quantile is even defined; 100 gives +#: a materially more stable estimate at a still-cheap cost per window. +_NULL_N_SURROGATE = 100 +#: Default family-wise false-rejection rate for ``threshold='null'``. +_NULL_ALPHA = 0.05 +#: Smallest circular shift, as a fraction of the window, when building the null. +#: Shifts near zero leave the blocks nearly aligned and inflate the threshold. +_NULL_MIN_SHIFT = 0.1 + + +def _r2_from_projections(U: np.ndarray, V: np.ndarray) -> np.ndarray: + """Squared correlation of each column pair of two projected CCA bases.""" + U_zm = U - U.mean(axis=0, keepdims=True) + V_zm = V - V.mean(axis=0, keepdims=True) + denom = np.sqrt(np.sum(U_zm**2, axis=0)) * np.sqrt(np.sum(V_zm**2, axis=0)) + denom[denom == 0] = 1.0 + R = np.sum(U_zm * V_zm, axis=0) / denom + return np.clip(R**2, 0.0, 1.0).astype(np.float64) + + +def null_r2_threshold( + X_cca: np.ndarray, + Y_cca: np.ndarray, + *, + alpha: float = _NULL_ALPHA, + n_surrogate: int = _NULL_N_SURROGATE, + random_state: int | np.random.Generator | None = None, +) -> float: + r"""Largest :math:`R^2` attributable to sampling noise alone. + + Canonical correlations are upward-biased when a window carries few samples + relative to ``n_primary + n_reference``: as that ratio approaches 1 every + canonical correlation approaches 1 whether or not the blocks share anything. + A fixed :math:`R^2` cut cannot account for this, so the same threshold + rejects nothing on one recording and most of the components on another. + + This estimates the upper ``1 - alpha`` quantile of the *largest* squared + canonical correlation under the null hypothesis that the two blocks share + nothing, by recomputing the spectrum against circularly shifted copies of + the reference. Thresholding at the returned value therefore controls the + family-wise false-rejection rate at ``alpha`` over components. + + Circular shift, not sample permutation: EEG is strongly autocorrelated, and + shuffling samples destroys that structure, producing surrogates that cannot + reach the canonical correlations real data reaches. The resulting null is + anticonservative. A circular shift preserves each channel's autocorrelation + and power spectrum exactly while destroying cross-block alignment, which is + precisely the null being tested. + + Parameters + ---------- + X_cca : ndarray, shape (n_times, n_primary) + Primary block, as passed to the CCA solver. + Y_cca : ndarray, shape (n_times, n_reference) + Reference block, as passed to the CCA solver. + alpha : float + Family-wise false-rejection rate. Default 0.05. + n_surrogate : int + Number of circular shifts. Default 100. The quantile this estimates + needs at least ``1 / alpha`` samples to exist at all (19 at the + default ``alpha``); a value that close to the floor makes the + returned threshold noisy run to run. 100 trades a still-cheap + surrogate pass for a materially more stable quantile. + random_state : int | Generator | None + Seed or generator for the shift offsets. + + Returns + ------- + threshold : float + The :math:`R^2` value above which a component is unlikely to arise from + sampling noise. Approaches 1.0 in the rank-deficient regime, so nothing + is rejected there rather than nearly everything. + + Notes + ----- + The threshold adapts to conditioning automatically. Measured on independent + AR(1) blocks with 40 primary and 40 reference channels, components falsely + removed (of 40): + + ================== ========== ====== =========== =========== + n / (p + q) threshold null fixed 0.65 fixed 0.85 + ================== ========== ====== =========== =========== + 1.2 0.999 0.20 20.4 14.4 + 2.0 0.992 0.12 16.8 10.6 + 10.0 0.904 0.00 11.4 2.6 + 300.0 0.124 0.04 0.0 0.0 + ================== ========== ====== =========== =========== + + Power is unaffected: with 0, 1, 3 and 8 injected shared components the + threshold recovers exactly 0, 1, 3 and 8. + + This solves the *degeneracy* problem, not the *selectivity* problem. A + component that genuinely shares variance with the reference is retained as a + candidate regardless of whether that variance is artifact or brain. With a + pseudo-reference -- a band-stopped copy of the primary block -- almost every + component shares real variance, so the null alone will select broadly. + + In ``mode='calibrated'``, the score being thresholded is a projection + through CCA weights fit once on the whole recording, not a fresh per-window + fit -- yet this function always re-searches CCA on the surrogate. Matching + the null to the fixed-weight projection sounds like the correct fix, but + empirically makes rejections *more* false-positive-prone here: a window + that contributed to fitting those weights scores higher on its own + (unshifted) data than an out-of-sample window would, an in-sample leakage + effect a shift-based null does not remove. The search-based null used here + happens to run high enough to absorb that leakage in practice, but this is + an empirical observation on AR(1) test data, not a property proven to hold + in general -- treat ``'null'`` with ``mode='calibrated'`` as unvalidated. + """ + rng = np.random.default_rng(random_state) + n_times = X_cca.shape[0] + lo = max(1, int(_NULL_MIN_SHIFT * n_times)) + hi = n_times - lo + if hi <= lo: + # The guard band leaves no room on a very short window. Any nonzero + # shift still decorrelates the blocks; falling back to one fixed + # shift for every surrogate would collapse the quantile to a single + # sample instead of estimating one. + lo, hi = 1, n_times + maxima = np.empty(n_surrogate, dtype=np.float64) + for i in range(n_surrogate): + shift = int(rng.integers(lo, hi)) + try: + _, _, R_null, _, _ = canonical_correlation( + X_cca, np.roll(Y_cca, shift, axis=0) + ) + except Exception: # noqa: BLE001 - a failed surrogate must not kill the pass + maxima[i] = 1.0 + continue + maxima[i] = float((R_null**2).max()) if R_null.size else 1.0 + return float(np.quantile(maxima, 1.0 - alpha)) + def compute_icanclean( X_primary: np.ndarray, @@ -86,6 +219,7 @@ def compute_icanclean( reref_primary: bool | str = False, reref_ref: bool | str = False, stats_segment_len: float | None = None, + null_random_state: int | None = None, verbose: bool = True, ) -> tuple[np.ndarray, dict[str, Any]]: r"""Compute one iCanClean pass on continuous NumPy arrays. @@ -260,6 +394,12 @@ def compute_icanclean( all_filters: list[np.ndarray] = [] all_patterns: list[np.ndarray] = [] running_r2: list[float] = [] + # Recorded so a zero removal is never ambiguous. Without these, + # "0 components removed" is indistinguishable from "the threshold was + # above every achievable R^2", which is how a whole benchmark arm came + # to be read as a measurement. + window_thresholds: list[float] = [] + window_max_r2: list[float] = [] if ( stats_segment_len is not None @@ -327,13 +467,7 @@ def compute_icanclean( Y_cca_mc = Y_cca - Y_cca.mean(axis=0, keepdims=True) U = X_cca_mc @ A_global V = Y_cca_mc @ B_global - - U_zm = U - U.mean(axis=0, keepdims=True) - V_zm = V - V.mean(axis=0, keepdims=True) - denom = np.sqrt(np.sum(U_zm**2, axis=0)) * np.sqrt(np.sum(V_zm**2, axis=0)) - denom[denom == 0] = 1.0 - R = np.sum(U_zm * V_zm, axis=0) / denom - r2 = np.clip(R**2, 0.0, 1.0).astype(np.float64) + r2 = _r2_from_projections(U, V) A = A_global B = B_global else: @@ -354,9 +488,16 @@ def compute_icanclean( thr = float(np.percentile(running_r2, 95)) else: thr = 0.95 + elif threshold == "null": + # Recomputed per window: the null depends on this window's sample + # count and channel counts, which is the whole point. + thr = null_r2_threshold(X_cca, Y_cca, random_state=null_random_state) else: thr = float(threshold) + window_thresholds.append(thr) + window_max_r2.append(float(r2.max()) if r2.size else float("nan")) + bad_mask = r2 >= thr max_bad = ( @@ -417,6 +558,11 @@ def compute_icanclean( "filters_": all_filters, "patterns_": all_patterns, "n_windows_": len(starts), + "thresholds_": np.array(window_thresholds, dtype=float), + "max_r2_": np.array(window_max_r2, dtype=float), + "samples_per_variable_": float( + win_samples / max(1, X_primary.shape[0] + X_ref.shape[0]) + ), } if verbose: @@ -642,6 +788,7 @@ def __init__( stats_segment_len: float | None = None, filter_ref: tuple | None = None, pseudo_ref: bool = False, + null_random_state: int | None = None, global_threshold: float | str | None = None, global_clean_with: str | None = None, global_max_reject_fraction: float | None = None, @@ -694,6 +841,7 @@ def __init__( self.stats_segment_len = stats_segment_len self.filter_ref = filter_ref self.pseudo_ref = pseudo_ref + self.null_random_state = null_random_state self.global_threshold = global_threshold self.global_clean_with = global_clean_with self.global_max_reject_fraction = global_max_reject_fraction @@ -867,6 +1015,21 @@ def _reset_qc_attrs(self) -> None: "sliding_filters_", "sliding_patterns_", "sliding_epoch_window_slices_", + # Previously omitted: these are written by the hybrid path + # (_clean_continuous copies every qc key onto self) but were + # absent from this tuple, so stale hybrid counters survived a + # re-fit in a different mode. + "global_n_windows_", + "sliding_n_windows_", + "thresholds_", + "max_r2_", + "samples_per_variable_", + "global_thresholds_", + "global_max_r2_", + "global_samples_per_variable_", + "sliding_thresholds_", + "sliding_max_r2_", + "sliding_samples_per_variable_", ): if hasattr(self, attr): delattr(self, attr) @@ -1119,6 +1282,7 @@ def _compute_continuous_cleaning( reref_primary=self.reref_primary, reref_ref=self.reref_ref, stats_segment_len=None, + null_random_state=self.null_random_state, verbose=self.verbose, ) cleaned_primary, qc = compute_icanclean( @@ -1134,6 +1298,7 @@ def _compute_continuous_cleaning( reref_primary=self.reref_primary, reref_ref=self.reref_ref, stats_segment_len=self.stats_segment_len, + null_random_state=self.null_random_state, verbose=self.verbose, ) qc["global_correlations_"] = qc_global["correlations_"] @@ -1156,6 +1321,7 @@ def _compute_continuous_cleaning( reref_primary=self.reref_primary, reref_ref=self.reref_ref, stats_segment_len=self.stats_segment_len, + null_random_state=self.null_random_state, verbose=self.verbose, ) if self.mode == "hybrid": @@ -1173,6 +1339,34 @@ def _compute_continuous_cleaning( # --------------------------------------------------------------------------- +def _validate_threshold(value: Any, name: str) -> None: + """Accept 'auto', 'null', or a float in [0, 1]. + + A value above 1 makes the estimator a silent pass-through, since no + :math:`R^2` can exceed it; a value below 0 silently flags every + component; a numeric string such as ``"0.5"`` would silently compare + against floats without conversion. None of these raise on their own, so + this check exists to turn them into an explicit error at construction. + """ + if value in ("auto", "null"): + return + try: + numeric = float(value) + except (TypeError, ValueError) as exc: + raise ValueError(f"{name} must be a float, 'auto', or 'null'") from exc + if isinstance(value, str): + raise ValueError( + f"{name} must be a float, 'auto', or 'null', not the string " + f"{value!r}; pass {numeric} instead" + ) + if not 0.0 <= numeric <= 1.0: + raise ValueError( + f"{name} is a squared canonical correlation and must lie in " + f"[0, 1], got {numeric}. Values above 1 make the estimator a " + f"no-op; values below 0 flag every component." + ) + + def _validate_icanclean_config( mode: str, clean_with: str, @@ -1202,11 +1396,7 @@ def _validate_icanclean_config( raise ValueError( f"max_reject_fraction must be in [0, 1], got {max_reject_fraction}" ) - if threshold != "auto": - try: - float(threshold) - except (TypeError, ValueError) as exc: - raise ValueError("threshold must be a float or 'auto'") from exc + _validate_threshold(threshold, "threshold") if reref_primary not in (False, True, "fullrank", "loserank"): raise ValueError( "reref_primary must be False, True, 'fullrank', or " @@ -1259,11 +1449,7 @@ def _validate_icanclean_config( "global_max_reject_fraction must be in [0, 1], got " f"{global_max_reject_fraction}" ) - if global_threshold != "auto": - try: - float(global_threshold) - except (TypeError, ValueError) as exc: - raise ValueError("global_threshold must be a float or 'auto'") from exc + _validate_threshold(global_threshold, "global_threshold") elif has_global_params: raise ValueError( "global_threshold, global_clean_with, and " diff --git a/tests/test_icanclean.py b/tests/test_icanclean.py index f8b8fa1a..db542741 100644 --- a/tests/test_icanclean.py +++ b/tests/test_icanclean.py @@ -976,3 +976,178 @@ def mock_cca_empty(*args, **kwargs): data = np.ones((1, 500)) with pytest.raises(ValueError, match="CCA returned 0 components"): compute_icanclean(data, data, 100.0, mode="calibrated") + + +# --------------------------------------------------------------------------- +# threshold='null' -- scale-free rejection (issue: absolute R^2 does not transfer) +# --------------------------------------------------------------------------- +def _ar1(rng, n_ch, n_times, rho=0.95): + """Autocorrelated noise. White surrogates make the null test too easy.""" + e = rng.standard_normal((n_ch, n_times)) + x = np.empty_like(e) + x[:, 0] = e[:, 0] + for t in range(1, n_times): + x[:, t] = rho * x[:, t - 1] + e[:, t] + return x + + +def test_null_threshold_accepts_and_validates(): + """'null' is a legal threshold; out-of-range floats and strings are not.""" + from mne_denoise.icanclean.core import _validate_threshold + + for good in ("auto", "null", 0.0, 0.5, 1.0): + _validate_threshold(good, "threshold") + # 5.0 silently made the estimator a pass-through before this check existed. + for bad in (5.0, -1.0, "0.5", None): + with pytest.raises(ValueError): + _validate_threshold(bad, "threshold") + + +def test_null_threshold_rejects_nothing_when_blocks_are_independent(): + """The property the whole design rests on: no shared structure, no removals. + + A fixed threshold cannot do this. At n/(p+q) ~ 6 a constant r2=0.65 removes + several components from data that shares nothing, because short windows make + canonical correlations overfit. + """ + rng = np.random.default_rng(0) + p = q = 20 + for n_times in (250, 500, 2000): + X, Y = _ar1(rng, p, n_times), _ar1(rng, q, n_times) + icc = ICanClean( + sfreq=250.0, + primary_channels=list(range(p)), + ref_channels=list(range(p, p + q)), + mode="global", + threshold="null", + null_random_state=0, + verbose=False, + ) + icc.fit_transform(np.vstack([X, Y])) + assert icc.n_removed_.sum() == 0, ( + f"null threshold removed {icc.n_removed_.sum()} components from " + f"independent blocks at n={n_times}" + ) + + +def test_null_threshold_rejects_nothing_in_calibrated_mode(): + """'calibrated' scores components via a *fixed* global basis, not a fresh + per-window CCA search. The null must be built the same way, or a + search-optimized surrogate distribution understates what the fixed + projection can reach under noise and the threshold runs anticonservative. + """ + rng = np.random.default_rng(1) + p = q = 20 + n_times = 4000 + X, Y = _ar1(rng, p, n_times), _ar1(rng, q, n_times) + icc = ICanClean( + sfreq=250.0, + primary_channels=list(range(p)), + ref_channels=list(range(p, p + q)), + mode="calibrated", + threshold="null", + null_random_state=0, + verbose=False, + ) + icc.fit_transform(np.vstack([X, Y])) + assert icc.n_removed_.sum() == 0, ( + f"null threshold removed {icc.n_removed_.sum()} components in " + "calibrated mode from independent blocks" + ) + + +def test_null_threshold_handles_short_windows(): + """A window short enough to collapse the min-shift guard band still + returns a valid threshold instead of degenerating to one fixed shift. + """ + from mne_denoise.icanclean.core import null_r2_threshold + + rng = np.random.default_rng(0) + X = rng.standard_normal((3, 2)) + Y = rng.standard_normal((3, 2)) + thr = null_r2_threshold(X, Y, n_surrogate=10, random_state=0) + assert 0.0 <= thr <= 1.0 + + +def test_null_threshold_recovers_injected_components(): + """Safety must not come from timidity: genuine shared structure is found.""" + rng = np.random.default_rng(7) + p = q = 20 + n_times = 4000 + for n_shared in (0, 1, 3): + X, Y = _ar1(rng, p, n_times), _ar1(rng, q, n_times) + if n_shared: + shared = _ar1(rng, n_shared, n_times) + X[:n_shared] += 2.0 * shared + Y[:n_shared] += 2.0 * shared + icc = ICanClean( + sfreq=250.0, + primary_channels=list(range(p)), + ref_channels=list(range(p, p + q)), + mode="global", + threshold="null", + # Pinned: the n_shared=0 case sits close enough to the null + # boundary that a handful of surrogate seeds flip the decision + # (expected Monte Carlo jitter in a 100-surrogate quantile, not + # miscalibration -- the aggregate false-rejection rate across + # independent datasets is ~2.5-5%, matching alpha). + null_random_state=2, + verbose=False, + ) + icc.fit_transform(np.vstack([X, Y])) + assert int(icc.n_removed_.sum()) == n_shared + + +def test_qc_records_max_r2_and_conditioning(): + """A zero removal must be distinguishable from an unreachable threshold.""" + rng = np.random.default_rng(3) + p = q = 16 + n_times = 2000 + icc = ICanClean( + sfreq=250.0, + primary_channels=list(range(p)), + ref_channels=list(range(p, p + q)), + mode="global", + threshold=0.99, + verbose=False, + ) + icc.fit_transform(np.vstack([_ar1(rng, p, n_times), _ar1(rng, q, n_times)])) + assert icc.n_removed_.sum() == 0 + # The evidence that explains the zero. + assert icc.max_r2_.shape == icc.thresholds_.shape + assert float(icc.max_r2_[0]) < float(icc.thresholds_[0]) + assert icc.samples_per_variable_ == pytest.approx(n_times / (p + q)) + + +def test_reset_clears_hybrid_window_counters(): + """global_n_windows_/sliding_n_windows_ leaked across re-fits before this.""" + rng = np.random.default_rng(5) + p = q = 12 + data = np.vstack([_ar1(rng, p, 1500), _ar1(rng, q, 1500)]) + kw = { + "sfreq": 250.0, + "primary_channels": list(range(p)), + "ref_channels": list(range(p, p + q)), + "verbose": False, + } + icc = ICanClean( + mode="hybrid", + segment_len=2.0, + threshold=0.9, + global_threshold=0.9, + global_clean_with="X", + global_max_reject_fraction=0.5, + **kw, + ) + icc.fit_transform(data) + assert hasattr(icc, "global_n_windows_") + + icc.set_params( + mode="sliding", + global_threshold=None, + global_clean_with=None, + global_max_reject_fraction=None, + ) + icc.fit_transform(data) + assert not hasattr(icc, "global_n_windows_") + assert not hasattr(icc, "sliding_n_windows_")