feat(eeg): add remove_artifacts_with_ica with EOG correlation rejection

Implements ICA-based EOG artifact removal using measure="correlation",
threshold=0.9 to reliably flag components on small (4-channel) fixtures
where the default z-score threshold is algebraically unreachable.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

src/pipelines/eeg_pipeline.py

@@ -14,6 +14,7 @@ from __future__ import annotations
 
 import mne
 import numpy as np
+from mne.preprocessing import ICA
 
 from src.core.logger import get_logger
 
@@ -73,3 +74,69 @@ def bandpass_filter(
     out.filter(l_freq=l_freq, h_freq=h_freq, picks="all", verbose="ERROR")
     logger.info("Bandpass filter applied: %.1f-%.1f Hz", l_freq, h_freq)
     return out
+
+
+def remove_artifacts_with_ica(
+    raw: mne.io.BaseRaw,
+    eog_ch_name: str | None = None,
+    n_components: int = 15,
+    random_state: int = 97,
+) -> mne.io.BaseRaw:
+    """Remove EOG-like artifacts using MNE's ICA + EOG correlation.
+
+    Fits an ICA decomposition on `raw`, finds components whose time courses
+    correlate (Pearson) with the named EOG channel via `find_bads_eog` using
+    `measure="correlation"`, marks them as "bad" and reconstructs the signal
+    without them. Returns a copy; the input `raw` is unchanged.
+
+    If `eog_ch_name` is None or no bad components are found, returns a copy of
+    `raw` unchanged. This keeps the function safe to call on recordings
+    without an EOG reference.
+
+    Args:
+        raw: Loaded, ideally bandpass-filtered, `mne.io.BaseRaw`.
+        eog_ch_name: Name of the EOG channel for correlation-based detection.
+            None disables auto-rejection.
+        n_components: Cap on ICA components. For small recordings, MNE will
+            silently cap this at the rank of the data.
+        random_state: Seed for ICA's underlying solver. Required for §4
+            Determinism.
+
+    Returns:
+        A copy of `raw` with EOG-correlated ICA components removed.
+    """
+    out = raw.copy()
+    if eog_ch_name is None or eog_ch_name not in out.ch_names:
+        logger.info("ICA skipped: no EOG channel reference provided")
+        return out
+
+    # Cap n_components at the rank of the data to avoid solver complaints
+    # on small synthetic fixtures.
+    n_eeg = len(mne.pick_types(out.info, eeg=True, meg=False))
+    safe_n = min(n_components, max(n_eeg - 1, 1))
+
+    ica = ICA(
+        n_components=safe_n,
+        random_state=random_state,
+        max_iter="auto",
+        method="fastica",
+        verbose="ERROR",
+    )
+    ica.fit(out, picks="eeg", verbose="ERROR")
+    # Use raw correlation (not z-score) so we can reliably flag artifact
+    # components on small recordings where n_components < 10 makes the
+    # default z-score threshold algebraically unreachable.
+    bad_idx, _ = ica.find_bads_eog(
+        out,
+        ch_name=eog_ch_name,
+        measure="correlation",
+        threshold=0.9,
+        verbose="ERROR",
+    )
+    ica.exclude = list(bad_idx)
+    logger.info(
+        "ICA fit: n_components=%d, EOG-correlated rejected=%d",
+        safe_n, len(ica.exclude),
+    )
+    ica.apply(out, verbose="ERROR")
+    return out

tests/pipelines/test_eeg_pipeline.py

@@ -10,6 +10,7 @@ import pytest
 from src.pipelines.eeg_pipeline import (
     bandpass_filter,
     is_valid_epoch,
+    remove_artifacts_with_ica,
 )
 
 
@@ -89,3 +90,61 @@ class TestBandpassFilter:
             bandpass_filter(raw, l_freq=40.0, h_freq=1.0)
         with pytest.raises(ValueError, match="must be strictly less than"):
             bandpass_filter(raw, l_freq=10.0, h_freq=10.0)
+
+
+class TestRemoveArtifactsWithIca:
+    def _load(self) -> mne.io.BaseRaw:
+        return mne.io.read_raw_fif(FIXTURE, preload=True, verbose="ERROR")
+
+    def test_returns_raw_instance(self) -> None:
+        raw = bandpass_filter(self._load(), l_freq=1.0, h_freq=40.0)
+        out = remove_artifacts_with_ica(
+            raw, eog_ch_name="EOG061", n_components=4, random_state=97,
+        )
+        assert isinstance(out, mne.io.BaseRaw)
+
+    def test_preserves_shape(self) -> None:
+        raw = bandpass_filter(self._load(), l_freq=1.0, h_freq=40.0)
+        before = raw.get_data().shape
+        out = remove_artifacts_with_ica(
+            raw, eog_ch_name="EOG061", n_components=4, random_state=97,
+        )
+        assert out.get_data().shape == before
+
+    def test_reduces_eog_correlation_on_frontal_channel(self) -> None:
+        """ICA must reduce correlation between EOG and Cz (the bleed channel)."""
+        raw = bandpass_filter(self._load(), l_freq=1.0, h_freq=40.0)
+        before = raw.get_data()
+        cz_idx = raw.ch_names.index("Cz")
+        eog_idx = raw.ch_names.index("EOG061")
+        corr_before = abs(np.corrcoef(before[cz_idx], before[eog_idx])[0, 1])
+
+        out = remove_artifacts_with_ica(
+            raw, eog_ch_name="EOG061", n_components=4, random_state=97,
+        )
+        after = out.get_data()
+        corr_after = abs(np.corrcoef(after[cz_idx], after[eog_idx])[0, 1])
+        # Allow for noise — but the dominant EOG bleed must be reduced.
+        assert corr_after < corr_before
+
+    def test_no_eog_channel_is_a_noop(self) -> None:
+        """Without an EOG reference, ICA can't auto-reject — should pass through."""
+        raw = bandpass_filter(self._load(), l_freq=1.0, h_freq=40.0)
+        out = remove_artifacts_with_ica(
+            raw, eog_ch_name=None, n_components=4, random_state=97,
+        )
+        # Identical shape; data approximately equal (no rejection happened).
+        assert out.get_data().shape == raw.get_data().shape
+        np.testing.assert_allclose(
+            out.get_data(), raw.get_data(), rtol=1e-6, atol=1e-12
+        )
+
+    def test_is_deterministic_with_seed(self) -> None:
+        raw = bandpass_filter(self._load(), l_freq=1.0, h_freq=40.0)
+        a = remove_artifacts_with_ica(
+            raw, eog_ch_name="EOG061", n_components=4, random_state=97,
+        )
+        b = remove_artifacts_with_ica(
+            raw, eog_ch_name="EOG061", n_components=4, random_state=97,
+        )
+        np.testing.assert_allclose(a.get_data(), b.get_data(), rtol=1e-12, atol=1e-15)