diff --git a/.github/workflows/autofix.yml b/.github/workflows/autofix.yml index eb46084c14f..78c93ba5e0d 100644 --- a/.github/workflows/autofix.yml +++ b/.github/workflows/autofix.yml @@ -94,7 +94,7 @@ jobs: # e.g., `full` includes `mne[full-no-qt]`, but the action won't recognise this, so we need to specify `full-no-qt` explicitly # ffmpeg is pinned to work around a conda solver issue (gh-14023) additional-dependencies: pip,mamba,conda,nomkl,noqt5,ffmpeg==8.1.2 - pip-dependencies: pymef + pip-dependencies: jamica,pymef requirements-overrides: PySide6==6.11.1,vtk==9.6.2 - uses: autofix-ci/action@c5b2d67aa2274e7b5a18224e8171550871fc7e4a if: success() || failure() diff --git a/doc/changes/dev/14247.newfeature.rst b/doc/changes/dev/14247.newfeature.rst new file mode 100644 index 00000000000..3e62c655687 --- /dev/null +++ b/doc/changes/dev/14247.newfeature.rst @@ -0,0 +1 @@ +Add single-model JAMICA support to :class:`mne.preprocessing.ICA` via ``method="jamica"``, by :newcontrib:`Sina Esmaeili`. diff --git a/doc/changes/names.inc b/doc/changes/names.inc index 6856ca7e2da..6f48b057789 100644 --- a/doc/changes/names.inc +++ b/doc/changes/names.inc @@ -425,6 +425,7 @@ .. _Simon Kornblith: https://simonster.com .. _Simon M. Hofmann: https://github.com/SHEscher .. _Simon-Shlomo Poil: https://github.com/simon-shlomo +.. _Sina Esmaeili: https://github.com/snesmaeili .. _Sondre Foslien: https://github.com/sondrfos .. _Sophie Herbst: https://github.com/SophieHerbst .. _Sourav Singh: https://github.com/souravsingh diff --git a/doc/references.bib b/doc/references.bib index 4a9a0b89f02..30c3cc8d930 100644 --- a/doc/references.bib +++ b/doc/references.bib @@ -57,6 +57,14 @@ @article{AblinEtAl2018 year = {2018} } +@techreport{PalmerEtAl2011, + author = {Palmer, Jason A. and Kreutz-Delgado, Ken and Makeig, Scott}, + institution = {Swartz Center for Computational Neuroscience}, + title = {{AMICA}: An Adaptive Mixture of Independent Component Analyzers with Shared Components}, + url = {https://sccn.ucsd.edu/~jason/amica_a.pdf}, + year = {2011} +} + @article{AllenEtAl2000, author = {Allen, Philip J. and Josephs, Oliver and Turner, Robert}, doi = {10.1006/nimg.2000.0599}, diff --git a/environment.yml b/environment.yml index cfd278d4599..ec0996b7593 100644 --- a/environment.yml +++ b/environment.yml @@ -69,5 +69,6 @@ dependencies: - vtk ==9.6.2 - xlrd - pip: + - jamica>=0.3.0 - pymef - pyobjc-framework-Cocoa>=5.2.0; platform_system == "Darwin" diff --git a/mne/preprocessing/ica.py b/mne/preprocessing/ica.py index 1c706979460..c8aebe7dbc5 100644 --- a/mne/preprocessing/ica.py +++ b/mne/preprocessing/ica.py @@ -72,6 +72,7 @@ _pl, _reject_data_segments, _require_version, + _soft_import, _validate_type, check_fname, check_random_state, @@ -188,7 +189,7 @@ def _check_for_unsupported_ica_channels(picks, info, allow_ref_meg=False): ) -_KNOWN_ICA_METHODS = ("fastica", "infomax", "picard") +_KNOWN_ICA_METHODS = ("fastica", "infomax", "jamica", "picard") def _rng_to_seed(rng): @@ -249,23 +250,27 @@ class ICA(ContainsMixin): type prior to the whitening by PCA. %(rng)s %(random_state_rng)s - method : 'fastica' | 'infomax' | 'picard' + method : 'fastica' | 'infomax' | 'jamica' | 'picard' The ICA method to use in the fit method. Use the ``fit_params`` argument to set additional parameters. Specifically, if you want Extended Infomax, set ``method='infomax'`` and ``fit_params=dict(extended=True)`` (this also works for ``method='picard'``). Defaults to ``'fastica'``. - For reference, see :footcite:`Hyvarinen1999,BellSejnowski1995,LeeEtAl1999,AblinEtAl2018`. + ``method='jamica'`` fits a single ICA model. For multi-model adaptive + mixture ICA and other advanced functionality, use the + `jamica package `__ directly. For + reference, see :footcite:`Hyvarinen1999,BellSejnowski1995,LeeEtAl1999,AblinEtAl2018,PalmerEtAl2011`. fit_params : dict | None Additional parameters passed to the ICA estimator as specified by ``method``. Allowed entries are determined by the various algorithm implementations: see :class:`~sklearn.decomposition.FastICA`, - :func:`~picard.picard`, :func:`~mne.preprocessing.infomax`. + :func:`~picard.picard`, :func:`~mne.preprocessing.infomax`, and the + ``jamica.amica`` function from the ``jamica`` package. max_iter : int | 'auto' Maximum number of iterations during fit. If ``'auto'``, it will set maximum iterations to ``1000`` for ``'fastica'`` - and to ``500`` for ``'infomax'`` or ``'picard'``. The actual number of - iterations it took :meth:`ICA.fit` to complete will be stored in the - ``n_iter_`` attribute. + and to ``500`` for ``'infomax'``, ``'jamica'``, or ``'picard'``. The + actual number of iterations it took :meth:`ICA.fit` to complete will be + stored in the ``n_iter_`` attribute. allow_ref_meg : bool Allow ICA on MEG reference channels. Defaults to False. @@ -501,7 +506,7 @@ def __init__( _check_option("max_iter", max_iter, ("auto",), "when str") if method == "fastica": max_iter = 1000 - elif method in ["infomax", "picard"]: + elif method in ["infomax", "jamica", "picard"]: max_iter = 500 fit_params.setdefault("max_iter", max_iter) self.max_iter = max_iter @@ -517,7 +522,9 @@ def _get_infos_for_repr(self): @dataclass class _InfosForRepr: fit_on: Literal["raw data", "epochs"] | None - fit_method: Literal["fastica", "infomax", "extended-infomax", "picard"] + fit_method: Literal[ + "fastica", "infomax", "extended-infomax", "jamica", "picard" + ] fit_params: dict[str, str | float] fit_n_iter: int | None fit_n_samples: int | None @@ -678,7 +685,6 @@ def fit( for method, mod in req_map.items(): if self.method == method: _require_version(mod, f"use method={repr(method)}") - _validate_type(inst, (BaseRaw, BaseEpochs), "inst", "Raw or Epochs") if np.isclose(inst.info["highpass"], 0.0): @@ -1008,6 +1014,21 @@ def _fit(self, data, fit_type): self.unmixing_matrix_ = W self.n_iter_ = n_iter + 1 # picard() starts counting at 0 del _, n_iter + elif self.method == "jamica": + jamica = _soft_import( + "jamica", "fitting ICA with method='jamica'", min_version="0.3.0" + ) + + _, W, _, n_iter = jamica.amica( + data[:, sel].T, + whiten=False, + return_n_iter=True, + random_state=_rng_to_seed(rng), + **self.fit_params, + ) + self.unmixing_matrix_ = W + self.n_iter_ = n_iter + del _, n_iter assert self.unmixing_matrix_.shape == (self.n_components_,) * 2 norms = self.pca_explained_variance_ stable = norms / norms[0] > 1e-6 # to be stable during pinv diff --git a/mne/preprocessing/tests/test_ica.py b/mne/preprocessing/tests/test_ica.py index dacca6dccac..d095154936c 100644 --- a/mne/preprocessing/tests/test_ica.py +++ b/mne/preprocessing/tests/test_ica.py @@ -4,8 +4,10 @@ import os import shutil +import sys from contextlib import nullcontext from pathlib import Path +from unittest.mock import MagicMock import matplotlib.pyplot as plt import numpy as np @@ -182,6 +184,94 @@ def test_ica_simple(method): assert amari_distance < 0.1 +def test_ica_jamica_multimodel_error(): + """Test that JAMICA's single-model error propagates through ICA.""" + pytest.importorskip("jamica", minversion="0.3.0") + rng = np.random.default_rng(0) + info = create_info(2, 100.0, "mag") + with info._unlock(): + info["highpass"] = 1.0 + raw = RawArray(rng.standard_normal((2, 100)), info, verbose=False) + ica = ICA( + n_components=2, + method="jamica", + max_iter=1, + fit_params=dict(num_models=2), + ) + with pytest.raises(ValueError, match=r"single AMICA model.*jamica\.AmicaICA"): + ica.fit(raw, verbose=False) + + +def test_ica_jamica(tmp_path, monkeypatch): + """Test the single-model JAMICA solver boundary and I/O.""" + jamica = pytest.importorskip("jamica", minversion="0.3.0") + original_amica = jamica.amica + amica_mock = MagicMock(wraps=original_amica) + monkeypatch.setattr(jamica, "amica", amica_mock) + + rng = np.random.default_rng(7) + sources = rng.laplace(size=(3, 600)) + mixing = np.array([[1.0, 0.4, -0.2], [0.3, 1.2, 0.5], [-0.4, 0.2, 0.9]]) + data = mixing @ sources + info = create_info(3, 200.0, "mag") + with info._unlock(): + info["highpass"] = 1.0 + raw = RawArray(data, info, verbose=False) + fit_params = dict(do_newton=False, chunk_size=None) + ica = ICA( + n_components=3, + method="jamica", + rng=42, + max_iter=2, + fit_params=fit_params, + ) + with pytest.warns(jamica.JamicaConvergenceWarning, match="did not converge"): + ica.fit(raw, verbose=False) + + assert amica_mock.call_count == 1 + (solver_x,), solver_kwargs = amica_mock.call_args + assert solver_kwargs == dict( + whiten=False, + return_n_iter=True, + random_state=42, + max_iter=2, + **fit_params, + ) + + # Explicitly reconstruct the matrix MNE passes to external ICA solvers. + expected_x = ica._pre_whiten(data.copy()) + expected_x = ica.pca_components_ @ (expected_x - ica.pca_mean_[:, np.newaxis]) + pca_norms = np.sqrt(ica.pca_explained_variance_) + expected_x /= pca_norms[:, np.newaxis] + assert_allclose(solver_x, expected_x[: ica.n_components_], rtol=1e-12, atol=1e-12) + + with pytest.warns(jamica.JamicaConvergenceWarning, match="did not converge"): + _, direct_w, _, direct_n_iter = original_amica(solver_x, **solver_kwargs) + # The same input and seed must produce the same operator, so ICA ambiguity + # does not apply to this direct-call equivalence check. + mne_solver_w = ica.unmixing_matrix_ * pca_norms[: ica.n_components_] + assert_allclose(mne_solver_w, direct_w, rtol=1e-12, atol=1e-12) + assert ica.n_iter_ == direct_n_iter == 2 + assert_allclose(ica.mixing_matrix_, linalg.pinv(ica.unmixing_matrix_)) + assert ica.get_sources(raw).get_data().shape == data.shape + + ica.exclude = [0] + fname = tmp_path / "jamica-ica.fif" + ica.save(fname) + with monkeypatch.context() as m: + m.setitem(sys.modules, "jamica", None) + ica_read = read_ica(fname) + assert ica_read.method == "jamica" + assert_allclose(ica_read.unmixing_matrix_, ica.unmixing_matrix_) + assert_allclose(ica_read.mixing_matrix_, ica.mixing_matrix_) + assert_allclose( + ica_read.get_sources(raw).get_data(), ica.get_sources(raw).get_data() + ) + assert_allclose( + ica_read.apply(raw.copy()).get_data(), ica.apply(raw.copy()).get_data() + ) + + def test_warnings(): """Test that ICA warns on certain input data conditions.""" raw = read_raw_fif(raw_fname).crop(0, 5).load_data() @@ -273,7 +363,8 @@ def test_ica_noop(n_components, n_pca_components, tmp_path): @pytest.mark.parametrize( - "method, max_iter_default", [("fastica", 1000), ("infomax", 500), ("picard", 500)] + "method, max_iter_default", + [("fastica", 1000), ("infomax", 500), ("jamica", 500), ("picard", 500)], ) def test_ica_max_iter_(method, max_iter_default): """Test that ICA.max_iter is set to the right defaults.""" diff --git a/mne/utils/config.py b/mne/utils/config.py index 5664b00c5ce..305b1112460 100644 --- a/mne/utils/config.py +++ b/mne/utils/config.py @@ -877,6 +877,7 @@ def sys_info( "dipy", "openmeeg", "python-picard", + "jamica", "cupy", "pandas", "h5io", diff --git a/pyproject.toml b/pyproject.toml index 2ea1ff516dc..711b4339d26 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -178,6 +178,7 @@ full-no-qt = [ "ipympl", "ipython >= 2.0, != 8.7.0", # for notebook backend; also in "doc" and "test" "ipywidgets", + "jamica >= 0.3.0", "joblib >= 0.8", "jupyter", "mffpy >= 0.11.0", diff --git a/tutorials/preprocessing/40_artifact_correction_ica.py b/tutorials/preprocessing/40_artifact_correction_ica.py index 4afc3d415ec..3d2fc6f7923 100644 --- a/tutorials/preprocessing/40_artifact_correction_ica.py +++ b/tutorials/preprocessing/40_artifact_correction_ica.py @@ -95,13 +95,18 @@ # dimensionality, set ``n_components=n`` during initialization and pass # ``n_pca_components=n`` to `~mne.preprocessing.ICA.apply`. # -# MNE-Python implements three different ICA algorithms: ``fastica`` (the -# default), ``picard``, and ``infomax``. FastICA and Infomax are both in fairly -# widespread use; Picard is a newer (2017) algorithm that is expected to -# converge faster than FastICA and Infomax, and is more robust than other -# algorithms in cases where the sources are not completely independent, which -# typically happens with real EEG/MEG data. See -# :footcite:`AblinEtAl2018` for more information. +# MNE-Python implements four different ICA algorithms: ``fastica`` (the +# default), ``picard``, ``infomax``, and ``jamica``. FastICA and Infomax are +# both in fairly widespread use; Picard is a newer (2017) algorithm that is +# expected to converge faster than FastICA and Infomax, and is more robust than +# other algorithms in cases where the sources are not completely independent, +# which typically happens with real EEG/MEG data. See +# :footcite:`AblinEtAl2018` for more information. JAMICA is an optional Python +# implementation of Adaptive Mixture ICA (AMICA). Within +# `~mne.preprocessing.ICA`, ``method='jamica'`` fits one ICA model; use the +# `jamica package `__ directly for +# multi-model decompositions and its other advanced functionality. See +# :footcite:`PalmerEtAl2011` for the AMICA method. # # The ICA interface in MNE-Python is similar to the interface in # `scikit-learn `__: some general parameters