From 8ef0ae5cfd1aa70478579e6a76a65567fd94d765 Mon Sep 17 00:00:00 2001 From: Seyed Yahya Shirazi Date: Sat, 18 Jul 2026 07:02:39 -0700 Subject: [PATCH 1/7] Add AMICAICA MNE wrapper (single-model) --- pamica/mne_compat/__init__.py | 15 ++ pamica/mne_compat/core.py | 287 ++++++++++++++++++++++++++++++++++ pyproject.toml | 10 +- uv.lock | 6 +- 4 files changed, 316 insertions(+), 2 deletions(-) create mode 100644 pamica/mne_compat/__init__.py create mode 100644 pamica/mne_compat/core.py diff --git a/pamica/mne_compat/__init__.py b/pamica/mne_compat/__init__.py new file mode 100644 index 0000000..26a5ad9 --- /dev/null +++ b/pamica/mne_compat/__init__.py @@ -0,0 +1,15 @@ +"""MNE-Python compatibility layer for pamica (issue #139). + +This subpackage is an *optional* entry point: it imports ``mne`` at module +import time, so it is intentionally NOT imported by :mod:`pamica.__init__`. +``import pamica`` therefore never requires MNE; reach the wrapper explicitly:: + + from pamica.mne_compat import AMICAICA + +Install the dependency with ``pip install pamica[mne]`` (same lazy-extra +pattern as the optional ``mlx`` backend). +""" + +from .core import AMICAICA + +__all__ = ["AMICAICA"] diff --git a/pamica/mne_compat/core.py b/pamica/mne_compat/core.py new file mode 100644 index 0000000..db808ce --- /dev/null +++ b/pamica/mne_compat/core.py @@ -0,0 +1,287 @@ +"""MNE-Python-facing wrapper over pamica's AMICA backend (issue #139, phase 1). + +:class:`AMICAICA` fits AMICA directly from an :class:`mne.io.Raw` / +:class:`mne.Epochs` and exposes the standard MNE ICA consumer surface +(``get_sources``, ``apply``, ``get_components``, ``plot_components``). Its core +value-add is :meth:`AMICAICA.to_mne_ica`, which maps pamica's fitted +mean/sphere/unmixing into a fully-populated :class:`mne.preprocessing.ICA`; the +other methods delegate to that object, so MNE performs the export/plot/back- +projection machinery unchanged. + +This phase covers the single-model case (``n_models == 1``). Multi-model +exposure is issue #141. +""" + +from typing import Optional, Union + +import numpy as np +import torch + +# MNE is an optional extra (`pip install pamica[mne]`); this subpackage is never +# imported by ``pamica.__init__``, so ``import pamica`` does not require it. The +# base type-check/CI env has no mne, hence the scoped ignore (issue #139). +import mne # ty: ignore[unresolved-import] +from mne.preprocessing import ICA as _MNEICA # ty: ignore[unresolved-import] + +from ..amica import AMICA + + +class AMICAICA: + """Fit AMICA from MNE objects and interoperate with ``mne.preprocessing.ICA``. + + The wrapper fits pamica's natural-gradient AMICA backend on the data of an + MNE :class:`~mne.io.Raw` or :class:`~mne.Epochs` and lets MNE consume the + result: :meth:`get_sources`, :meth:`apply`, :meth:`get_components` and + :meth:`plot_components` all delegate to a real + :class:`mne.preprocessing.ICA` built by :meth:`to_mne_ica`. + + Parameters + ---------- + n_mix : int, default=3 + Number of mixture components per source (AMICA ``n_mix``). + random_state : int or None, default=None + Seed for the AMICA fit (passed through as the backend ``seed``) and + stored on the exported :class:`~mne.preprocessing.ICA`. + device : str or torch.device, optional + Torch device for the fit (``None`` = auto; the float64 parity backend + falls back to CPU when auto-selection lands on MPS). See :class:`AMICA`. + verbose : bool, default=True + Whether the underlying :class:`AMICA` prints fit progress. + + Attributes + ---------- + amica_ : AMICA + The fitted single-model pamica model. + info_ : mne.Info + The picked measurement info the fit was run on (channel subset only). + ch_names_ : list of str + Names of the fitted channels, in order. + n_components_ : int + Number of ICA components (equals the number of fitted channels; AMICA + keeps ``n_sources == n_channels``). + + Notes + ----- + The single-model AMICA transform is + ``S = W_fort @ sphere @ (X - mean)`` (the per-model data-space center ``c`` + is exactly zero for one model). MNE computes sources as + ``S = unmixing_matrix_ @ pca_components_ @ (X - pca_mean_)`` (with a unit + pre-whitener). Writing the symmetric-ZCA ``sphere`` as + ``V @ diag(1/sqrt(e)) @ V.T`` with ``V`` orthonormal, the exported ICA uses + ``pca_components_ = V.T``, ``unmixing_matrix_ = W_fort @ sphere @ V`` and + ``pca_mean_ = mean``. Keeping ``pca_components_`` orthonormal is what makes + MNE's ``get_components`` (scalp maps ``inv(sphere) @ inv(W_fort)``) come out + right, since MNE assumes orthonormal PCA rows. The mapping is pinned by a + round-trip test (``to_mne_ica().get_sources(raw)`` equals + ``amica_.transform(X)``). + """ + + def __init__( + self, + n_mix: int = 3, + random_state: Optional[int] = None, + device: Optional[Union[str, torch.device]] = None, + verbose: bool = True, + ): + self.n_mix = n_mix + self.random_state = random_state + self.device = device + self.verbose = verbose + + self.amica_: Optional[AMICA] = None + self.info_ = None + self.ch_names_: Optional[list] = None + self.n_components_: Optional[int] = None + self._fit_kind: Optional[str] = None + self._mne_ica: Optional[_MNEICA] = None + + # ------------------------------------------------------------------ + # Fit + # ------------------------------------------------------------------ + def fit( + self, + inst, + picks=None, + start: Optional[int] = None, + stop: Optional[int] = None, + **fit_kwargs, + ) -> "AMICAICA": + """Fit AMICA to the data of an MNE ``Raw`` or ``Epochs``. + + Parameters + ---------- + inst : mne.io.BaseRaw | mne.BaseEpochs + The data to decompose. ``Epochs`` are concatenated along time + (``np.hstack``), matching how MNE's own ICA fits epoched data. + picks : str | list | slice | None, default=None + Channels to fit, in any form MNE accepts (e.g. ``"eeg"``, + ``"data"``, a name list). ``None`` selects all good data channels + (bads excluded), matching MNE's ICA default. + start, stop : int | None + Sample range (``Raw`` only) passed to ``get_data``; ``None`` uses + the full recording. + **fit_kwargs + Forwarded to :meth:`AMICA.fit` (e.g. ``max_iter``, ``lrate``, + ``do_newton``) and the backend constructor (e.g. ``block_size``). + + Returns + ------- + self : AMICAICA + """ + if not isinstance(inst, (mne.io.BaseRaw, mne.BaseEpochs)): + raise TypeError( + "AMICAICA.fit expects an mne.io.Raw or mne.Epochs, got " + f"{type(inst).__name__}." + ) + + picked = inst.copy().pick("data" if picks is None else picks, exclude="bads") + self.info_ = picked.info + self.ch_names_ = list(picked.ch_names) + + if isinstance(inst, mne.io.BaseRaw): + range_kw = {} + if start is not None: + range_kw["start"] = start + if stop is not None: + range_kw["stop"] = stop + X = picked.get_data(**range_kw) + self._fit_kind = "raw" + else: + # Concatenate epochs along time, as MNE's ICA does for fitting. + X = np.hstack(picked.get_data()) + self._fit_kind = "epochs" + + X = np.ascontiguousarray(X, dtype=np.float64) + self.n_components_ = X.shape[0] + + if isinstance(self.random_state, (int, np.integer)): + fit_kwargs.setdefault("seed", int(self.random_state)) + + amica = AMICA( + n_models=1, n_mix=self.n_mix, device=self.device, verbose=self.verbose + ) + amica.fit(X, **fit_kwargs) + + self.amica_ = amica + self._mne_ica = None # invalidate any cached export + return self + + # ------------------------------------------------------------------ + # Export + # ------------------------------------------------------------------ + def to_mne_ica(self) -> _MNEICA: + """Build (and cache) a fully-populated :class:`mne.preprocessing.ICA`. + + The returned object is a genuine MNE ICA: ``get_sources``, ``apply``, + ``get_components`` and the ``mne.viz`` plotters operate on it natively. + See the class :class:`Notes ` for the mean/sphere/unmixing to + ``pca_mean_``/``pca_components_``/``unmixing_matrix_`` mapping. + + Returns + ------- + ica : mne.preprocessing.ICA + """ + self._check_fitted() + if self._mne_ica is not None: + return self._mne_ica + + assert ( + self.amica_ is not None + and self.amica_.model_ is not None + and self.ch_names_ is not None + ) + backend = self.amica_.model_ + assert ( + backend.mean is not None + and backend.sphere is not None + and backend.c is not None + ) + mean = backend.mean.cpu().numpy().ravel() + sphere = backend.sphere.cpu().numpy() + w_fort = self.amica_.get_unmixing_matrix(model_idx=0) + c = backend.c.cpu().numpy()[:, 0] + n_ch = sphere.shape[0] + + # Orthonormal eigenbasis of the symmetric-ZCA sphere + # (sphere = V diag(1/sqrt(cov_eval)) V.T). eigh gives ascending + # sphere-eigenvalues (= 1/sqrt(cov_eval)); reorder to descending + # explained variance so pca_components_ matches MNE's PCA convention. + sphere_evals, evecs = np.linalg.eigh(sphere) + cov_evals = 1.0 / sphere_evals**2 + order = np.argsort(cov_evals)[::-1] + v = evecs[:, order] + cov_evals = cov_evals[order] + + pca_components = v.T + unmixing = w_fort @ sphere @ v + # c is the per-model center in sphered space; fold it into pca_mean via + # the data-space offset inv(sphere) @ c. For a single model c is exactly + # zero, so this leaves pca_mean == mean bit-for-bit (no eigh residual). + if np.any(c): + pca_mean = mean + np.linalg.solve(sphere, c) + else: + pca_mean = mean + + ica = _MNEICA( + n_components=n_ch, + method="infomax", + max_iter="auto", + random_state=self.random_state, + ) + ica.info = self.info_ + ica.ch_names = list(self.ch_names_) + ica.n_components_ = n_ch + ica.pca_mean_ = pca_mean + ica.pca_components_ = pca_components + ica.pca_explained_variance_ = cov_evals + ica.unmixing_matrix_ = unmixing + ica.pre_whitener_ = np.ones((n_ch, 1)) + ica.n_iter_ = max(int(getattr(backend, "iteration", 0)), 1) + ica._update_mixing_matrix() + ica._update_ica_names() + ica.current_fit = self._fit_kind + + self._mne_ica = ica + return ica + + # ------------------------------------------------------------------ + # MNE consumer surface (delegates to the exported ICA) + # ------------------------------------------------------------------ + def get_sources(self, inst): + """Sources for ``inst`` as an MNE object (see ``ICA.get_sources``).""" + return self.to_mne_ica().get_sources(inst) + + def apply(self, inst, **kwargs): + """Remove selected components and back-project (see ``ICA.apply``). + + Pass ``exclude=[...]`` (or set it on the exported ICA) to drop + components; with no exclusions this reconstructs the input. + """ + return self.to_mne_ica().apply(inst, **kwargs) + + def get_components(self) -> np.ndarray: + """Scalp maps, shape ``(n_channels, n_components)`` (``ICA.get_components``).""" + return self.to_mne_ica().get_components() + + def plot_components(self, *args, **kwargs): + """Plot component topographies (see ``ICA.plot_components``).""" + return self.to_mne_ica().plot_components(*args, **kwargs) + + def plot_sources(self, inst, *args, **kwargs): + """Plot component time courses (see ``ICA.plot_sources``).""" + return self.to_mne_ica().plot_sources(inst, *args, **kwargs) + + # ------------------------------------------------------------------ + # Helpers + # ------------------------------------------------------------------ + def _check_fitted(self) -> None: + if self.amica_ is None: + raise RuntimeError("AMICAICA must be fitted before this call; run fit().") + + def __repr__(self) -> str: + if self.amica_ is None: + return f"" + return ( + f"" + ) diff --git a/pyproject.toml b/pyproject.toml index db086d0..528437d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -49,6 +49,10 @@ docs = [ # Visualization for the benchmark sweeps: MNE builds the scalp-map montage for # the IC topomaps and (issue #91) the distributed channel-subset positions. viz = ["mne>=1.6"] +# MNE-Python compatibility layer (AMICAICA, issue #139). MNE is not a core +# dependency, so install with `pip install pamica[mne]`; `import pamica` never +# requires it (pamica.mne_compat is imported lazily, same pattern as mlx). +mne = ["mne>=1.6"] [project.urls] Homepage = "https://github.com/sccn/pAMICA" @@ -57,7 +61,7 @@ Repository = "https://github.com/sccn/pAMICA" [tool.setuptools] # List subpackages explicitly so the wheel is deterministic across platforms # (a bare ["pamica"] dropped torch_impl/ on some setuptools versions). -packages = ["pamica", "pamica.numpy_impl", "pamica.torch_impl", "pamica.mlx_impl", "pamica.metrics", "pamica.native"] +packages = ["pamica", "pamica.numpy_impl", "pamica.torch_impl", "pamica.mlx_impl", "pamica.metrics", "pamica.native", "pamica.mne_compat"] # numpy_impl/params.json is the NumPy backend's default parameter file, loaded at # runtime via importlib/__file__, so it must ship in the wheel. package-data = {"pamica" = ["data/*"], "pamica.numpy_impl" = ["params.json"]} @@ -83,6 +87,10 @@ omit = [ # GPU, so CI (ubuntu, no mlx) cannot exercise it and always measures 0%. It is # covered locally by tests/mlx_tests/ on Apple hardware (issue #76). "pamica/mlx_impl/*", + # Optional MNE compatibility layer (issue #139): imports mne, absent from the + # base test env, so the base coverage run always measures 0%. Exercised by + # the dedicated `test-mne` CI job (uv sync --extra mne) on real sample EEG. + "pamica/mne_compat/*", ] [tool.coverage.report] diff --git a/uv.lock b/uv.lock index 7fd617f..ed80b15 100644 --- a/uv.lock +++ b/uv.lock @@ -1212,6 +1212,9 @@ docs = [ mlx = [ { name = "mlx" }, ] +mne = [ + { name = "mne" }, +] viz = [ { name = "mne" }, ] @@ -1234,6 +1237,7 @@ requires-dist = [ { name = "mkdocs-material", marker = "extra == 'docs'", specifier = ">=9.5" }, { name = "mkdocstrings", extras = ["python"], marker = "extra == 'docs'", specifier = ">=0.26" }, { name = "mlx", marker = "extra == 'mlx'", specifier = ">=0.32" }, + { name = "mne", marker = "extra == 'mne'", specifier = ">=1.6" }, { name = "mne", marker = "extra == 'viz'", specifier = ">=1.6" }, { name = "numpy" }, { name = "scipy" }, @@ -1241,7 +1245,7 @@ requires-dist = [ { name = "torch", specifier = ">=2.12.1" }, { name = "tqdm" }, ] -provides-extras = ["mlx", "docs", "viz"] +provides-extras = ["mlx", "docs", "viz", "mne"] [package.metadata.requires-dev] dev = [ From 0b804e48a0a7ad776b063a56333dd53f31930a31 Mon Sep 17 00:00:00 2001 From: Seyed Yahya Shirazi Date: Sat, 18 Jul 2026 07:02:40 -0700 Subject: [PATCH 2/7] Add MNE wrapper tests and CI job --- .github/workflows/ci.yml | 22 +++ pamica/tests/mne_tests/__init__.py | 0 pamica/tests/mne_tests/test_mne_wrapper.py | 161 +++++++++++++++++++++ 3 files changed, 183 insertions(+) create mode 100644 pamica/tests/mne_tests/__init__.py create mode 100644 pamica/tests/mne_tests/test_mne_wrapper.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e00ed75..a50793d 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -113,6 +113,28 @@ jobs: files: coverage.xml fail_ci_if_error: false + test-mne: + name: Test MNE wrapper (optional [mne] extra) + needs: [lint, typecheck] + runs-on: ubuntu-latest + env: + PYTORCH_ENABLE_MPS_FALLBACK: "1" + steps: + - uses: actions/checkout@v4 + - uses: astral-sh/setup-uv@v5 + with: + python-version: "3.12" + enable-cache: true + # mne IS pip-installable on Linux (unlike the Apple-only mlx extra), so the + # AMICAICA wrapper (issue #139) is exercised for real in CI here. Its tests + # live in pamica/tests/mne_tests and self-skip in the base `test` job + # (importorskip on mne); this job installs the [mne] extra so they run on + # the real sample EEG. pamica/mne_compat is omitted from coverage in + # pyproject (the base env has no mne), so no fail-under gate applies here. + - run: uv sync --extra mne + - name: pytest (mne wrapper, real sample EEG) + run: uv run pytest pamica/tests/mne_tests + build: name: Build + import (Python ${{ matrix.python-version }}) needs: [lint, typecheck] diff --git a/pamica/tests/mne_tests/__init__.py b/pamica/tests/mne_tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/pamica/tests/mne_tests/test_mne_wrapper.py b/pamica/tests/mne_tests/test_mne_wrapper.py new file mode 100644 index 0000000..c223b05 --- /dev/null +++ b/pamica/tests/mne_tests/test_mne_wrapper.py @@ -0,0 +1,161 @@ +"""AMICAICA MNE-wrapper tests (issue #139 phase 1, single-model). + +Real sample EEG only (no synthetic/mock): the bundled EEGLAB ``eeglab_data.set`` +(32 channels, 30504 samples) is loaded through ``mne.io.read_raw_eeglab`` -- the +same entry point a real MNE user would use. The whole module self-skips when the +optional ``mne`` extra is absent, so the base CI env (no mne) skips it; the +dedicated ``test-mne`` CI job installs ``pamica[mne]`` and runs it for real. + +The load-bearing check is the round trip: the ICA that ``to_mne_ica`` builds must +reproduce ``AMICA.transform`` bit-for-bit (up to float64 eigh residual), which +pins the mean/sphere/unmixing -> pca_mean_/pca_components_/unmixing_matrix_ map. +""" + +from pathlib import Path + +import numpy as np +import pytest + +mne = pytest.importorskip("mne") + +from pamica.mne_compat import AMICAICA # noqa: E402 (after importorskip) + +mne.set_log_level("ERROR") + +SAMPLE_DIR = Path(__file__).resolve().parents[2] / "sample_data" +SET_FILE = SAMPLE_DIR / "eeglab_data.set" +SEED = 42 +MAX_ITER = 12 # enough to move off the init; parity is orientation, not convergence + +pytestmark = pytest.mark.skipif( + not SET_FILE.exists(), reason="sample eeglab_data.set missing" +) + + +@pytest.fixture(scope="module") +def raw(): + """Real continuous EEG as an MNE Raw (32 EEG channels, 128 Hz).""" + return mne.io.read_raw_eeglab(str(SET_FILE), preload=True) + + +@pytest.fixture(scope="module") +def fitted(raw): + """A single AMICAICA fit reused across the read-only assertions.""" + return AMICAICA(n_mix=3, random_state=SEED, device="cpu", verbose=False).fit( + raw, max_iter=MAX_ITER + ) + + +def _picked_data(raw): + """The exact array AMICAICA fits on (all good data channels, float64).""" + return raw.copy().pick("data", exclude="bads").get_data().astype(np.float64) + + +# --- the mapping crux ------------------------------------------------------- +def test_get_sources_matches_amica_transform(raw, fitted): + """to_mne_ica().get_sources == AMICA.transform on the same data.""" + s_mne = fitted.to_mne_ica().get_sources(raw).get_data() + s_amica = fitted.amica_.transform(_picked_data(raw)) + assert s_mne.shape == s_amica.shape == (fitted.n_components_, raw.n_times) + np.testing.assert_allclose(s_mne, s_amica, rtol=1e-6, atol=1e-9) + + +def test_get_components_are_channel_space_maps(fitted): + """get_components() equals the channel-space mixing inv(sphere) @ inv(W).""" + comps = fitted.get_components() + assert comps.shape == (fitted.n_components_, fitted.n_components_) + sphere = fitted.amica_.model_.sphere.cpu().numpy() + w_fort = fitted.amica_.get_unmixing_matrix(0) + maps_ref = np.linalg.inv(sphere) @ np.linalg.inv(w_fort) + np.testing.assert_allclose(comps, maps_ref, rtol=1e-6, atol=1e-10) + + +def test_apply_without_exclude_reconstructs(raw, fitted): + """apply() with no excluded components returns the input unchanged.""" + recon = fitted.apply(raw.copy()) + np.testing.assert_allclose(recon.get_data(), raw.get_data(), rtol=1e-6, atol=1e-12) + + +def test_apply_excludes_component(raw, fitted): + """Excluding a component changes the data and zeroes that source.""" + ica = fitted.to_mne_ica() + cleaned = ica.apply(raw.copy(), exclude=[0]) + # The reconstruction must differ from the input (a component was removed). + assert not np.allclose(cleaned.get_data(), raw.get_data()) + # Re-deriving sources from the cleaned data: component 0 is gone. + src_after = ica.get_sources(cleaned).get_data() + assert np.abs(src_after[0]).max() < 1e-6 * np.abs(src_after).max() + + +# --- object validity -------------------------------------------------------- +def test_to_mne_ica_returns_valid_fitted_ica(raw, fitted): + ica = fitted.to_mne_ica() + assert isinstance(ica, mne.preprocessing.ICA) + assert ica.current_fit != "unfitted" + assert ica.n_components_ == fitted.n_components_ + assert ica.ch_names == fitted.ch_names_ + # A genuine MNE ICA exposes the full component interface. + assert ica.get_components().shape == (fitted.n_components_, fitted.n_components_) + + +def test_export_is_cached_and_invalidated_on_refit(raw): + ica = AMICAICA(random_state=SEED, device="cpu", verbose=False).fit( + raw, max_iter=MAX_ITER + ) + first = ica.to_mne_ica() + assert ica.to_mne_ica() is first # cached + ica.fit(raw, max_iter=MAX_ITER) # refit invalidates the cache + assert ica.to_mne_ica() is not first + + +# --- picks / epochs / plotting ---------------------------------------------- +def test_fit_with_channel_subset(raw): + picks = raw.ch_names[:10] + ica = AMICAICA(random_state=SEED, device="cpu", verbose=False).fit( + raw, picks=picks, max_iter=MAX_ITER + ) + assert ica.n_components_ == 10 + assert ica.ch_names_ == picks + s_mne = ica.to_mne_ica().get_sources(raw).get_data() + x = raw.copy().pick(picks).get_data().astype(np.float64) + assert ica.amica_ is not None + np.testing.assert_allclose(s_mne, ica.amica_.transform(x), rtol=1e-6, atol=1e-9) + + +def test_fit_from_epochs(raw): + epochs = mne.make_fixed_length_epochs(raw, duration=2.0, preload=True) + ica = AMICAICA(random_state=SEED, device="cpu", verbose=False).fit( + epochs, max_iter=MAX_ITER + ) + assert ica._fit_kind == "epochs" + src = ica.to_mne_ica().get_sources(epochs).get_data() # (n_ep, n_comp, n_time) + x = np.hstack(epochs.copy().pick("data", exclude="bads").get_data()) + assert ica.amica_ is not None + np.testing.assert_allclose( + np.hstack(src), ica.amica_.transform(x), rtol=1e-6, atol=1e-9 + ) + + +def test_plot_components_returns_figure(fitted): + import matplotlib + + matplotlib.use("Agg") + import matplotlib.figure + import matplotlib.pyplot as plt + + out = fitted.plot_components(picks=[0, 1], show=False) + figs = out if isinstance(out, list) else [out] + assert figs and all(isinstance(f, matplotlib.figure.Figure) for f in figs) + plt.close("all") + + +# --- guards ----------------------------------------------------------------- +def test_unfitted_calls_raise(): + ica = AMICAICA() + with pytest.raises(RuntimeError, match="must be fitted"): + ica.to_mne_ica() + + +def test_fit_rejects_non_mne_input(): + with pytest.raises(TypeError, match="mne.io.Raw or mne.Epochs"): + AMICAICA().fit(np.random.randn(8, 500)) From 3dc6cb0e19689a29b8ba369255e2a1230f58c481 Mon Sep 17 00:00:00 2001 From: Seyed Yahya Shirazi Date: Sat, 18 Jul 2026 07:02:41 -0700 Subject: [PATCH 3/7] Document AMICAICA MNE wrapper --- docs/api/index.md | 10 +++++++ docs/api/mne-compat.md | 63 ++++++++++++++++++++++++++++++++++++++++++ docs/changelog.md | 20 ++++++++++++++ mkdocs.yml | 1 + 4 files changed, 94 insertions(+) create mode 100644 docs/api/mne-compat.md diff --git a/docs/api/index.md b/docs/api/index.md index e9b7e45..3be3169 100644 --- a/docs/api/index.md +++ b/docs/api/index.md @@ -27,3 +27,13 @@ from pamica.mlx_impl import AMICAMLXNG # requires the `mlx` extra - **[`AMICAMLXNG`](mlx-backend.md)** — the optional Apple-GPU (MLX) backend; the fastest option on Apple Silicon (float32). + +The optional MNE-Python wrapper is likewise imported explicitly: + +```python +from pamica.mne_compat import AMICAICA # requires the `mne` extra +``` + +- **[`AMICAICA`](mne-compat.md)** — fit AMICA from an MNE `Raw`/`Epochs` and + interoperate with `mne.preprocessing.ICA` (`get_sources`, `apply`, + `plot_components`, `to_mne_ica`). diff --git a/docs/api/mne-compat.md b/docs/api/mne-compat.md new file mode 100644 index 0000000..e44c37b --- /dev/null +++ b/docs/api/mne-compat.md @@ -0,0 +1,63 @@ +# MNE-Python compatibility (AMICAICA) + +`AMICAICA` fits AMICA directly from an [MNE-Python](https://mne.tools) +`Raw`/`Epochs` and hands the result back through the standard MNE ICA surface. +It is **additive**: the scikit-learn-style [`AMICA`](amica.md) interface and the +byte-identical [EEGLAB output](../guides/eeglab.md) are unchanged; this is a +second entry point for MNE users, not a replacement. + +MNE is an optional dependency, so `import pamica` never requires it. Install the +extra and import the wrapper explicitly: + +```bash +pip install pamica[mne] +``` + +```python +import mne +from pamica.mne_compat import AMICAICA + +raw = mne.io.read_raw_eeglab("subject.set", preload=True) + +ica = AMICAICA(n_mix=3, random_state=42).fit(raw, picks="eeg", max_iter=100) + +sources = ica.get_sources(raw) # an mne.io.RawArray of component activations +maps = ica.get_components() # scalp maps, shape (n_channels, n_components) +ica.plot_components() # native mne.viz topographies + +# Reconstruct with some components removed: +clean = ica.apply(raw.copy(), exclude=[0, 3]) +``` + +`fit` accepts a `Raw` or `Epochs` (epochs are concatenated along time, as MNE's +own ICA does), any MNE `picks` selector, and forwards remaining keywords +(`max_iter`, `lrate`, `do_newton`, ...) to [`AMICA.fit`](amica.md). + +## Interoperating with `mne.preprocessing.ICA` + +`to_mne_ica()` returns a fully-populated +[`mne.preprocessing.ICA`](https://mne.tools/stable/generated/mne.preprocessing.ICA.html), +so the entire MNE ICA ecosystem (plotting, `find_bads_eog`/`_ecg`, exclusion +workflows) works on an AMICA decomposition: + +```python +mne_ica = ica.to_mne_ica() +eog_idx, scores = mne_ica.find_bads_eog(raw) +mne_ica.plot_scores(scores) +``` + +The wrapper's `get_sources`, `apply`, `get_components` and `plot_components` +delegate to this object, so they reproduce `AMICA.transform` exactly: MNE +computes sources as `unmixing_matrix_ @ pca_components_ @ (X - pca_mean_)`, and +the export maps pamica's mean, symmetric-ZCA sphere and unmixing into those +matrices (writing the sphere as `V diag(1/√e) Vᵀ` with `V` orthonormal so MNE's +scalp maps come out in channel space). The equivalence +`to_mne_ica().get_sources(raw) == AMICA.transform(X)` is pinned by the test +suite on real sample EEG. + +!!! note "Single-model in this release" + `AMICAICA` covers the single-model case (`n_models == 1`). Multi-model + exposure for MNE users (per-model sources and model dominance) is tracked in + issue #141. + +::: pamica.mne_compat.AMICAICA diff --git a/docs/changelog.md b/docs/changelog.md index fc562a2..4f3188e 100644 --- a/docs/changelog.md +++ b/docs/changelog.md @@ -3,6 +3,26 @@ Release notes are also published on the [GitHub releases page](https://github.com/sccn/pAMICA/releases). +## 0.3.0 + +MNE-Python compatibility layer (epic #139), additive: the scikit-learn-style +`AMICA` API and the byte-identical EEGLAB I/O are unchanged. + +- `pamica.mne_compat.AMICAICA`, an MNE-facing wrapper that fits AMICA directly + from an `mne.io.Raw`/`Epochs` (`picks=...`, epochs concatenated along time like + MNE's own ICA) and interoperates with the standard MNE ICA consumer surface: + `get_sources`, `apply`, `get_components` and `plot_components`. `to_mne_ica()` + returns a fully-populated `mne.preprocessing.ICA`, so the whole MNE ICA + ecosystem (component plotting, `find_bads_eog`/`_ecg`, exclusion workflows) + works on an AMICA decomposition. The export maps pamica's mean, symmetric-ZCA + sphere and unmixing into MNE's `pca_mean_`/`pca_components_`/`unmixing_matrix_`, + writing the sphere as `V diag(1/sqrt(e)) V^T` with `V` orthonormal so MNE's + scalp maps are in channel space; `to_mne_ica().get_sources(raw)` reproduces + `AMICA.transform(X)` to float64 precision, pinned on real sample EEG. MNE is an + optional extra (`pip install pamica[mne]`); `import pamica` never requires it, + and a dedicated CI job runs the wrapper tests with the extra installed (phase 1, + single-model, #140). + ## 0.2.2 GitHub repository rename to pAMICA and a `__version__` fix. diff --git a/mkdocs.yml b/mkdocs.yml index bce8084..65307d3 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -102,6 +102,7 @@ nav: - MLX backend: api/mlx-backend.md - NumPy backend: api/numpy-backend.md - Native backend: api/native-backend.md + - MNE compatibility: api/mne-compat.md - Separation metrics: api/metrics.md - Visualization: api/viz.md - Development: From 77bb8141709f8b0671526e741a9e54c6590f05f6 Mon Sep 17 00:00:00 2001 From: Seyed Yahya Shirazi Date: Sat, 18 Jul 2026 07:22:34 -0700 Subject: [PATCH 4/7] Harden AMICAICA: atomic fit, guards, degenerate refusal --- pamica/mne_compat/core.py | 182 ++++++++++++++++++++++++++++---------- 1 file changed, 134 insertions(+), 48 deletions(-) diff --git a/pamica/mne_compat/core.py b/pamica/mne_compat/core.py index db808ce..de0ad0f 100644 --- a/pamica/mne_compat/core.py +++ b/pamica/mne_compat/core.py @@ -2,14 +2,14 @@ :class:`AMICAICA` fits AMICA directly from an :class:`mne.io.Raw` / :class:`mne.Epochs` and exposes the standard MNE ICA consumer surface -(``get_sources``, ``apply``, ``get_components``, ``plot_components``). Its core -value-add is :meth:`AMICAICA.to_mne_ica`, which maps pamica's fitted -mean/sphere/unmixing into a fully-populated :class:`mne.preprocessing.ICA`; the -other methods delegate to that object, so MNE performs the export/plot/back- -projection machinery unchanged. - -This phase covers the single-model case (``n_models == 1``). Multi-model -exposure is issue #141. +(``get_sources``, ``apply``, ``get_components``, ``plot_components``, +``plot_sources``). Its core value-add is :meth:`AMICAICA.to_mne_ica`, which maps +pamica's fitted mean/sphere/unmixing into a fully-populated +:class:`mne.preprocessing.ICA`; the other methods delegate to that object, so MNE +performs the export/plot/back-projection machinery unchanged. + +This phase covers the single-model, full-rank case (``n_models == 1``, no PCA +reduction). Multi-model exposure is issue #141. """ from typing import Optional, Union @@ -31,8 +31,8 @@ class AMICAICA: The wrapper fits pamica's natural-gradient AMICA backend on the data of an MNE :class:`~mne.io.Raw` or :class:`~mne.Epochs` and lets MNE consume the - result: :meth:`get_sources`, :meth:`apply`, :meth:`get_components` and - :meth:`plot_components` all delegate to a real + result: :meth:`get_sources`, :meth:`apply`, :meth:`get_components`, + :meth:`plot_components` and :meth:`plot_sources` all delegate to a real :class:`mne.preprocessing.ICA` built by :meth:`to_mne_ica`. Parameters @@ -59,12 +59,17 @@ class AMICAICA: n_components_ : int Number of ICA components (equals the number of fitted channels; AMICA keeps ``n_sources == n_channels``). + converged_ : bool + Whether the last fit ended usable (not degenerate). A degenerate fit is + kept for inspection but refused by the consumer methods (issue #50). + stop_reason_ : str or None + Why the backend fit stopped (e.g. ``"max_iter"``, ``"nan_ll"``). Notes ----- - The single-model AMICA transform is - ``S = W_fort @ sphere @ (X - mean)`` (the per-model data-space center ``c`` - is exactly zero for one model). MNE computes sources as + The single-model AMICA transform is ``S = W_fort @ sphere @ (X - mean)`` (the + per-model data-space center ``c`` is identically zero for one model, since + its update is gated to ``n_models > 1``). MNE computes sources as ``S = unmixing_matrix_ @ pca_components_ @ (X - pca_mean_)`` (with a unit pre-whitener). Writing the symmetric-ZCA ``sphere`` as ``V @ diag(1/sqrt(e)) @ V.T`` with ``V`` orthonormal, the exported ICA uses @@ -74,6 +79,13 @@ class AMICAICA: right, since MNE assumes orthonormal PCA rows. The mapping is pinned by a round-trip test (``to_mne_ica().get_sources(raw)`` equals ``amica_.transform(X)``). + + PCA reduction (``pcakeep``/``pcadb``) is unsupported: it leaves the sphere + rank-deficient, so the export (which assumes a full-rank whitening) would be + numerically invalid, and :meth:`fit` rejects it. Rank-deficient *data* (for + example average-referenced EEG) is the same hazard reached through data + conditioning rather than a keyword; such a fit typically diverges and is + refused as degenerate. """ def __init__( @@ -92,6 +104,9 @@ def __init__( self.info_ = None self.ch_names_: Optional[list] = None self.n_components_: Optional[int] = None + self.converged_: bool = False + self.stop_reason_: Optional[str] = None + self._n_samples: Optional[int] = None self._fit_kind: Optional[str] = None self._mne_ica: Optional[_MNEICA] = None @@ -118,41 +133,65 @@ def fit( ``"data"``, a name list). ``None`` selects all good data channels (bads excluded), matching MNE's ICA default. start, stop : int | None - Sample range (``Raw`` only) passed to ``get_data``; ``None`` uses - the full recording. + Sample range for ``Raw`` input, passed to ``get_data``; ``None`` uses + the full recording. Not supported for ``Epochs`` (raises). **fit_kwargs Forwarded to :meth:`AMICA.fit` (e.g. ``max_iter``, ``lrate``, ``do_newton``) and the backend constructor (e.g. ``block_size``). + ``pcakeep``/``pcadb`` (PCA reduction) are rejected, since they leave + the sphere rank-deficient and the MNE export invalid. Returns ------- self : AMICAICA + + Raises + ------ + TypeError + If ``inst`` is not an MNE ``Raw``/``Epochs``. + ValueError + If ``start``/``stop`` are given for ``Epochs``, ``stop`` exceeds the + recording length, the selected data is non-finite, or the fit used + PCA reduction. """ if not isinstance(inst, (mne.io.BaseRaw, mne.BaseEpochs)): raise TypeError( "AMICAICA.fit expects an mne.io.Raw or mne.Epochs, got " f"{type(inst).__name__}." ) + is_raw = isinstance(inst, mne.io.BaseRaw) + if not is_raw and (start is not None or stop is not None): + raise ValueError("start/stop are only supported for Raw input, not Epochs.") picked = inst.copy().pick("data" if picks is None else picks, exclude="bads") - self.info_ = picked.info - self.ch_names_ = list(picked.ch_names) - - if isinstance(inst, mne.io.BaseRaw): + ch_names = list(picked.ch_names) + + if is_raw: + if stop is not None and stop > inst.n_times: + raise ValueError( + f"stop={stop} exceeds the recording length ({inst.n_times} " + "samples)." + ) range_kw = {} if start is not None: range_kw["start"] = start if stop is not None: range_kw["stop"] = stop X = picked.get_data(**range_kw) - self._fit_kind = "raw" + fit_kind = "raw" else: # Concatenate epochs along time, as MNE's ICA does for fitting. X = np.hstack(picked.get_data()) - self._fit_kind = "epochs" + fit_kind = "epochs" X = np.ascontiguousarray(X, dtype=np.float64) - self.n_components_ = X.shape[0] + if not np.isfinite(X).all(): + bad = [ch_names[i] for i in np.flatnonzero(~np.isfinite(X).all(axis=1))] + raise ValueError( + "AMICAICA.fit: input contains non-finite (NaN/Inf) samples in " + f"channel(s) {bad}; clean bad segments/interpolation before " + "fitting." + ) if isinstance(self.random_state, (int, np.integer)): fit_kwargs.setdefault("seed", int(self.random_state)) @@ -162,7 +201,26 @@ def fit( ) amica.fit(X, **fit_kwargs) + if amica.model_ is not None and amica.model_._pca_reduced(): + raise ValueError( + "AMICAICA does not support PCA reduction (pcakeep/pcadb): it " + "leaves the sphere rank-deficient, so the MNE ICA export (which " + "assumes a full-rank whitening) would be numerically invalid. " + "Refit without pcakeep/pcadb." + ) + + # Publish to self only after every fallible step above succeeds, so a + # failed (re)fit leaves the previously fitted state intact rather than a + # mix of the old model and the new attempt's metadata (mirrors the + # local-first pattern in AMICA.fit). + self.info_ = picked.info + self.ch_names_ = ch_names + self.n_components_ = X.shape[0] + self._n_samples = X.shape[1] + self._fit_kind = fit_kind self.amica_ = amica + self.converged_ = amica.converged_ + self.stop_reason_ = amica.stop_reason_ self._mne_ica = None # invalidate any cached export return self @@ -173,33 +231,38 @@ def to_mne_ica(self) -> _MNEICA: """Build (and cache) a fully-populated :class:`mne.preprocessing.ICA`. The returned object is a genuine MNE ICA: ``get_sources``, ``apply``, - ``get_components`` and the ``mne.viz`` plotters operate on it natively. - See the class :class:`Notes ` for the mean/sphere/unmixing to - ``pca_mean_``/``pca_components_``/``unmixing_matrix_`` mapping. + ``get_components``, ``save`` and the ``mne.viz`` plotters operate on it + natively. See the class :class:`Notes ` for the + mean/sphere/unmixing to ``pca_mean_``/``pca_components_``/ + ``unmixing_matrix_`` mapping. + + The object is cached and returned by reference: mutating it (for example + setting ``.exclude``) persists across subsequent :meth:`apply`/ + :meth:`get_sources` calls until the next :meth:`fit`. Returns ------- ica : mne.preprocessing.ICA """ - self._check_fitted() + self._check_fitted("build the MNE ICA") if self._mne_ica is not None: return self._mne_ica - assert ( - self.amica_ is not None - and self.amica_.model_ is not None - and self.ch_names_ is not None - ) - backend = self.amica_.model_ - assert ( - backend.mean is not None - and backend.sphere is not None - and backend.c is not None - ) + amica = self.amica_ + if amica is None or amica.model_ is None or self.ch_names_ is None: + raise RuntimeError( + "AMICAICA: internal state is inconsistent; refit before to_mne_ica()." + ) + backend = amica.model_ + if backend.mean is None or backend.sphere is None: + raise RuntimeError( + "AMICAICA: the fitted backend is missing mean/sphere; refit " + "before to_mne_ica()." + ) + mean = backend.mean.cpu().numpy().ravel() sphere = backend.sphere.cpu().numpy() - w_fort = self.amica_.get_unmixing_matrix(model_idx=0) - c = backend.c.cpu().numpy()[:, 0] + w_fort = amica.get_unmixing_matrix(model_idx=0) n_ch = sphere.shape[0] # Orthonormal eigenbasis of the symmetric-ZCA sphere @@ -214,13 +277,9 @@ def to_mne_ica(self) -> _MNEICA: pca_components = v.T unmixing = w_fort @ sphere @ v - # c is the per-model center in sphered space; fold it into pca_mean via - # the data-space offset inv(sphere) @ c. For a single model c is exactly - # zero, so this leaves pca_mean == mean bit-for-bit (no eigh residual). - if np.any(c): - pca_mean = mean + np.linalg.solve(sphere, c) - else: - pca_mean = mean + # Single-model AMICA has the per-model center c identically zero (its + # update is gated to n_models>1), so no data-space offset enters pca_mean. + pca_mean = mean ica = _MNEICA( n_components=n_ch, @@ -228,6 +287,9 @@ def to_mne_ica(self) -> _MNEICA: max_iter="auto", random_state=self.random_state, ) + # `method` is inert here: the fitted attributes are hand-populated and + # ICA.fit (the only place `method` branches) is never called, but + # ICA.__init__ still requires a valid method name. ica.info = self.info_ ica.ch_names = list(self.ch_names_) ica.n_components_ = n_ch @@ -237,6 +299,11 @@ def to_mne_ica(self) -> _MNEICA: ica.unmixing_matrix_ = unmixing ica.pre_whitener_ = np.ones((n_ch, 1)) ica.n_iter_ = max(int(getattr(backend, "iteration", 0)), 1) + # MNE's own fit sets these; read_ica_eeglab (the precedent for building + # an ICA from an external decomposition) sets reject_=None. Without them + # ICA.save()/plot_properties raise AttributeError. + ica.reject_ = None + ica.n_samples_ = int(self._n_samples) if self._n_samples is not None else 0 ica._update_mixing_matrix() ica._update_ica_names() ica.current_fit = self._fit_kind @@ -274,13 +341,32 @@ def plot_sources(self, inst, *args, **kwargs): # ------------------------------------------------------------------ # Helpers # ------------------------------------------------------------------ - def _check_fitted(self) -> None: + def _check_fitted(self, action: str = "this call") -> None: + """Raise if no usable model is available. + + ``ValueError`` when never fitted (matching :class:`AMICA`'s convention); + ``RuntimeError`` when the fit ended degenerate (non-finite parameters, + issue #50), so the failure surfaces here rather than as opaque NaNs + downstream. + """ if self.amica_ is None: - raise RuntimeError("AMICAICA must be fitted before this call; run fit().") + raise ValueError(f"AMICAICA must be fitted before {action}; run fit().") + if not self.converged_: + raise RuntimeError( + f"Refusing to {action}: the AMICA fit ended degenerate " + f"(stop_reason={self.stop_reason_!r}), so it holds non-finite " + "parameters and would produce NaN output. Lower lrate, disable " + "Newton, or check data conditioning, then refit." + ) def __repr__(self) -> str: if self.amica_ is None: return f"" + if not self.converged_: + return ( + f"" + ) return ( f"" From 9d8fff34b3b64a484ccba85572d059adc767b9a0 Mon Sep 17 00:00:00 2001 From: Seyed Yahya Shirazi Date: Sat, 18 Jul 2026 07:22:35 -0700 Subject: [PATCH 5/7] Test eigen-order, guards, save-reload, refit safety --- pamica/tests/mne_tests/test_mne_wrapper.py | 140 ++++++++++++++++++++- 1 file changed, 139 insertions(+), 1 deletion(-) diff --git a/pamica/tests/mne_tests/test_mne_wrapper.py b/pamica/tests/mne_tests/test_mne_wrapper.py index c223b05..8b7f0a9 100644 --- a/pamica/tests/mne_tests/test_mne_wrapper.py +++ b/pamica/tests/mne_tests/test_mne_wrapper.py @@ -149,13 +149,151 @@ def test_plot_components_returns_figure(fitted): plt.close("all") +def test_pca_components_orthonormal_and_variance_ordered(raw, fitted): + """The eigenbasis ordering (invisible to the V-cancelling round trip) is real. + + get_sources/get_components reduce to W@sphere and inv(sphere)@inv(W) for ANY + orthonormal V, so they cannot see an eigenvector mis-ordering. Pin it here: + pca_components_ must be orthonormal and its rows ordered by descending data + variance, with pca_explained_variance_ equal to that per-row variance. + """ + ica = fitted.to_mne_ica() + p = ica.pca_components_ + np.testing.assert_allclose(p @ p.T, np.eye(p.shape[0]), atol=1e-10) + ev = ica.pca_explained_variance_ + assert np.all(np.diff(ev) <= 0), "explained variance must be descending" + # Independent oracle: project the mean-centered data onto pca_components_; + # each row's variance is that data-covariance eigenvalue, in the same order. + x = _picked_data(raw) + proj = p @ (x - ica.pca_mean_[:, None]) + np.testing.assert_allclose(proj.var(axis=1), ev, rtol=1e-5) + + +def test_fit_with_start_stop_subranges_raw(raw): + stop = 5000 + ica = AMICAICA(random_state=SEED, device="cpu", verbose=False).fit( + raw, start=100, stop=stop, max_iter=MAX_ITER + ) + assert ica.to_mne_ica().n_samples_ == stop - 100 + # A different range yields a different decomposition (the range is honored). + full = AMICAICA(random_state=SEED, device="cpu", verbose=False).fit( + raw, max_iter=MAX_ITER + ) + assert not np.allclose(ica.get_components(), full.get_components()) + + +def test_apply_on_epochs(raw): + epochs = mne.make_fixed_length_epochs(raw, duration=2.0, preload=True) + ica = AMICAICA(random_state=SEED, device="cpu", verbose=False).fit( + epochs, max_iter=MAX_ITER + ) + recon = ica.apply(epochs.copy()) + np.testing.assert_allclose( + recon.get_data(), epochs.get_data(), rtol=1e-6, atol=1e-12 + ) + cleaned = ica.apply(epochs.copy(), exclude=[0]) + assert not np.allclose(cleaned.get_data(), epochs.get_data()) + + +def test_exported_ica_saves_and_reloads(fitted, raw, tmp_path): + """reject_/n_samples_ must be set, else ICA.save() raises AttributeError.""" + ica = fitted.to_mne_ica() + fname = tmp_path / "amica-ica.fif" + ica.save(fname) + reloaded = mne.preprocessing.read_ica(fname) + np.testing.assert_allclose( + reloaded.get_sources(raw).get_data(), + ica.get_sources(raw).get_data(), + rtol=1e-6, + atol=1e-9, + ) + + +def test_plot_sources_returns_figure(fitted, raw): + import matplotlib + + matplotlib.use("Agg") + import matplotlib.figure + import matplotlib.pyplot as plt + + fig = fitted.plot_sources(raw, show=False) + assert isinstance(fig, matplotlib.figure.Figure) + plt.close("all") + + +def test_seed_determinism(raw): + a = AMICAICA(random_state=1, device="cpu", verbose=False).fit( + raw, max_iter=MAX_ITER + ) + b = AMICAICA(random_state=1, device="cpu", verbose=False).fit( + raw, max_iter=MAX_ITER + ) + np.testing.assert_allclose(a.get_components(), b.get_components()) + c = AMICAICA(random_state=2, device="cpu", verbose=False).fit( + raw, max_iter=MAX_ITER + ) + assert not np.allclose(a.get_components(), c.get_components()) + + # --- guards ----------------------------------------------------------------- def test_unfitted_calls_raise(): ica = AMICAICA() - with pytest.raises(RuntimeError, match="must be fitted"): + with pytest.raises(ValueError, match="must be fitted"): ica.to_mne_ica() def test_fit_rejects_non_mne_input(): with pytest.raises(TypeError, match="mne.io.Raw or mne.Epochs"): AMICAICA().fit(np.random.randn(8, 500)) + + +def test_fit_rejects_start_stop_on_epochs(raw): + epochs = mne.make_fixed_length_epochs(raw, duration=2.0, preload=True) + with pytest.raises(ValueError, match="only supported for Raw"): + AMICAICA().fit(epochs, stop=100) + + +def test_fit_rejects_out_of_range_stop(raw): + with pytest.raises(ValueError, match="exceeds the recording length"): + AMICAICA().fit(raw, stop=raw.n_times + 1) + + +def test_fit_rejects_non_finite_data(raw): + bad = raw.copy() + bad._data[3, 100:200] = np.nan + with pytest.raises(ValueError, match="non-finite"): + AMICAICA(device="cpu", verbose=False).fit(bad, max_iter=MAX_ITER) + + +def test_fit_rejects_pca_reduction(raw): + with pytest.raises(ValueError, match="PCA reduction"): + AMICAICA(device="cpu", verbose=False).fit(raw, max_iter=5, pcakeep=10) + + +def test_degenerate_fit_is_refused_and_labeled(raw): + """A degenerate (non-converged) fit is kept for inspection but refused. + + Real divergence is platform-dependent, so mark the wrapper's own + convergence signal (the documented contract) to exercise the guard + deterministically: `__repr__` labels it and the consumer methods refuse it. + """ + ica = AMICAICA(random_state=SEED, device="cpu", verbose=False).fit( + raw, max_iter=MAX_ITER + ) + ica.converged_ = False + ica.stop_reason_ = "nan_ll" + assert "degenerate" in repr(ica) + with pytest.raises(RuntimeError, match="degenerate"): + ica.to_mne_ica() + + +def test_failed_refit_leaves_prior_state_intact(raw): + """An exception mid-refit must not corrupt a previously fitted wrapper.""" + ica = AMICAICA(random_state=SEED, device="cpu", verbose=False).fit( + raw, max_iter=MAX_ITER + ) + before = ica.get_components() + with pytest.raises(ValueError): # stop past the end aborts before publishing + ica.fit(raw, picks=raw.ch_names[:10], stop=raw.n_times + 1, max_iter=MAX_ITER) + assert ica.n_components_ == before.shape[0] # still the 32-channel fit + np.testing.assert_array_equal(ica.get_components(), before) From 8365412d234dfb2a38261e6ecadd6dcdeec4bfc8 Mon Sep 17 00:00:00 2001 From: Seyed Yahya Shirazi Date: Sat, 18 Jul 2026 07:22:36 -0700 Subject: [PATCH 6/7] Document plot_sources, PCA-reduction limit, guards --- docs/api/mne-compat.md | 10 +++++++--- docs/changelog.md | 11 ++++++++--- 2 files changed, 15 insertions(+), 6 deletions(-) diff --git a/docs/api/mne-compat.md b/docs/api/mne-compat.md index e44c37b..7750a0b 100644 --- a/docs/api/mne-compat.md +++ b/docs/api/mne-compat.md @@ -31,7 +31,10 @@ clean = ica.apply(raw.copy(), exclude=[0, 3]) `fit` accepts a `Raw` or `Epochs` (epochs are concatenated along time, as MNE's own ICA does), any MNE `picks` selector, and forwards remaining keywords -(`max_iter`, `lrate`, `do_newton`, ...) to [`AMICA.fit`](amica.md). +(`max_iter`, `lrate`, `do_newton`, ...) to [`AMICA.fit`](amica.md). It rejects +non-finite input and PCA reduction (`pcakeep`/`pcadb`, which leaves the sphere +rank-deficient so the full-rank export would be invalid), and a degenerate +(diverged) fit is refused by the consumer methods rather than emitting NaNs. ## Interoperating with `mne.preprocessing.ICA` @@ -46,8 +49,9 @@ eog_idx, scores = mne_ica.find_bads_eog(raw) mne_ica.plot_scores(scores) ``` -The wrapper's `get_sources`, `apply`, `get_components` and `plot_components` -delegate to this object, so they reproduce `AMICA.transform` exactly: MNE +The wrapper's `get_sources`, `apply`, `get_components`, `plot_components` and +`plot_sources` delegate to this object, so they reproduce `AMICA.transform` +exactly: MNE computes sources as `unmixing_matrix_ @ pca_components_ @ (X - pca_mean_)`, and the export maps pamica's mean, symmetric-ZCA sphere and unmixing into those matrices (writing the sphere as `V diag(1/√e) Vᵀ` with `V` orthonormal so MNE's diff --git a/docs/changelog.md b/docs/changelog.md index 4f3188e..55d35bd 100644 --- a/docs/changelog.md +++ b/docs/changelog.md @@ -11,14 +11,19 @@ MNE-Python compatibility layer (epic #139), additive: the scikit-learn-style - `pamica.mne_compat.AMICAICA`, an MNE-facing wrapper that fits AMICA directly from an `mne.io.Raw`/`Epochs` (`picks=...`, epochs concatenated along time like MNE's own ICA) and interoperates with the standard MNE ICA consumer surface: - `get_sources`, `apply`, `get_components` and `plot_components`. `to_mne_ica()` - returns a fully-populated `mne.preprocessing.ICA`, so the whole MNE ICA + `get_sources`, `apply`, `get_components`, `plot_components` and `plot_sources`. + `to_mne_ica()` + returns a fully-populated `mne.preprocessing.ICA` (including `reject_`/ + `n_samples_`, so `ICA.save` and `plot_properties` work), so the whole MNE ICA ecosystem (component plotting, `find_bads_eog`/`_ecg`, exclusion workflows) works on an AMICA decomposition. The export maps pamica's mean, symmetric-ZCA sphere and unmixing into MNE's `pca_mean_`/`pca_components_`/`unmixing_matrix_`, writing the sphere as `V diag(1/sqrt(e)) V^T` with `V` orthonormal so MNE's scalp maps are in channel space; `to_mne_ica().get_sources(raw)` reproduces - `AMICA.transform(X)` to float64 precision, pinned on real sample EEG. MNE is an + `AMICA.transform(X)` to float64 precision, pinned on real sample EEG. `fit` + rejects PCA reduction (`pcakeep`/`pcadb`, which leaves the sphere rank-deficient + and the export invalid) and non-finite input, and a degenerate fit is refused + by the consumer methods rather than emitting NaNs. MNE is an optional extra (`pip install pamica[mne]`); `import pamica` never requires it, and a dedicated CI job runs the wrapper tests with the extra installed (phase 1, single-model, #140). From 0f1a97d1654cd3f83273b5f4c9aec9ef98626d84 Mon Sep 17 00:00:00 2001 From: Seyed Yahya Shirazi Date: Sat, 18 Jul 2026 07:24:32 -0700 Subject: [PATCH 7/7] Fix typo flagged by typos CI --- pamica/tests/mne_tests/test_mne_wrapper.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pamica/tests/mne_tests/test_mne_wrapper.py b/pamica/tests/mne_tests/test_mne_wrapper.py index 8b7f0a9..7ea150a 100644 --- a/pamica/tests/mne_tests/test_mne_wrapper.py +++ b/pamica/tests/mne_tests/test_mne_wrapper.py @@ -153,7 +153,7 @@ def test_pca_components_orthonormal_and_variance_ordered(raw, fitted): """The eigenbasis ordering (invisible to the V-cancelling round trip) is real. get_sources/get_components reduce to W@sphere and inv(sphere)@inv(W) for ANY - orthonormal V, so they cannot see an eigenvector mis-ordering. Pin it here: + orthonormal V, so they cannot see a wrong eigenvector order. Pin it here: pca_components_ must be orthonormal and its rows ordered by descending data variance, with pca_explained_variance_ equal to that per-row variance. """