ENH: Add JAMICA as an ICA method

.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()

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`.

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

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},

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"

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 <https://snesmaeili.github.io/jamica/>`__ 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

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."""

mne/utils/config.py

@@ -877,6 +877,7 @@ def sys_info(
         "dipy",
         "openmeeg",
         "python-picard",
+        "jamica",
         "cupy",
         "pandas",
         "h5io",

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",

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 <https://snesmaeili.github.io/jamica/>`__ 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 <https://scikit-learn.org/stable/>`__: some general parameters