Phase 2: multi-model exposure for MNE users

docs/api/mne-compat.md

@@ -59,9 +59,31 @@ 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.
+## Multi-model fits
+
+AMICA can learn a mixture of ICA models (`n_models > 1`). MNE's `ICA` represents
+only one unmixing matrix, so each model is exported as its own single-model
+`mne.preprocessing.ICA`, and the per-sample *model dominance* (which model best
+explains each timepoint) is exposed directly, since MNE has no concept for it:
+
+```python
+ica = AMICAICA(n_models=2, random_state=42).fit(raw, max_iter=100)
+
+# Per-model: model_idx selects the model on every consumer method.
+sources_m1 = ica.get_sources(raw, model_idx=1)
+ica.plot_components(model_idx=1)
+model1 = ica.to_mne_ica(model_idx=1)   # a standard ICA for model 1
+
+# Model dominance over time (P(model | sample), columns sum to 1):
+prob = ica.get_model_probability(raw)  # (n_models, n_samples)
+ica.plot_model_probability(raw)        # per-model probability + best-model LL
+```
+
+`get_model_probability`/`plot_model_probability` build on the public
+`AMICA.model_loglik`/`model_probability` accessors, which score arbitrary data
+through the fitted sphere and mean. Each per-model export folds that model's
+data-space center into `pca_mean_`, so `to_mne_ica(model_idx=h).get_sources(raw)`
+reproduces `AMICA.transform(X, model_idx=h)` (with `X` the picked channel array)
+for every model, not just the first.
 
 ::: pamica.mne_compat.AMICAICA

docs/changelog.md

@@ -27,6 +27,21 @@ MNE-Python compatibility layer (epic #139), additive: the scikit-learn-style
   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).
+- Multi-model exposure through the MNE wrapper: `AMICAICA(n_models=...)` fits a
+  mixture of ICA models, and since MNE's `ICA` represents only one unmixing,
+  each model is exported as its own single-model `mne.preprocessing.ICA` via
+  `to_mne_ica(model_idx=...)` (and the `model_idx` argument on `get_sources`/
+  `apply`/`get_components`/`plot_components`/`plot_sources`). The per-sample model
+  dominance MNE cannot represent is exposed directly: `get_model_probability(inst)`
+  returns `P(model | sample)` (`(n_models, n_samples)`, columns sum to 1) and
+  `plot_model_probability(inst)` draws the per-model probability plus best-model
+  log-likelihood over time. These build on a new public live accessor,
+  `AMICA.model_loglik`/`model_probability` (and the `AMICATorchNG` equivalents),
+  which score arbitrary data through the stored sphere/mean; the training-data
+  path (without `do_reject`) is pinned bit-for-bit against the E-step's own `Lht`. The per-model export
+  folds each model's data-space center `c` into `pca_mean_`, so the round trip
+  holds for the multi-model case too. `pamica.viz.plot_model_probability` now also
+  accepts a live `lht` array, not only a written `AmicaOutput` (phase 2, #141).
 
 ## 0.2.2
 

pamica/amica.py

@@ -419,6 +419,64 @@ def pmi(
 
         return self.model_.pmi(X, model_idx=model_idx, nbins=nbins)
 
+    def model_loglik(self, X: np.ndarray) -> np.ndarray:
+        """Per-model, per-sample log-likelihood ``Lht`` on ``X`` (issue #141).
+
+        Delegates to :meth:`AMICATorchNG.model_loglik`. For a multi-model fit
+        this is the joint log-likelihood of each model at each sample, from
+        which the per-sample model posterior (dominance) is
+        ``softmax(Lht, axis=0)``; see :meth:`model_probability`.
+
+        Parameters
+        ----------
+        X : np.ndarray of shape (n_channels, n_samples)
+            Raw (unpreprocessed) data.
+
+        Returns
+        -------
+        Lht : np.ndarray of shape (n_models, n_samples)
+
+        Raises
+        ------
+        ValueError
+            If the model is unfitted, or if ``X`` is non-finite.
+        RuntimeError
+            If the fit ended degenerate (issue #50).
+        """
+        self._check_usable("compute the model log-likelihood")
+        assert self.model_ is not None
+
+        return self.model_.model_loglik(X)
+
+    def model_probability(self, X: np.ndarray) -> np.ndarray:
+        """Per-sample posterior probability of each model (issue #141).
+
+        Delegates to :meth:`AMICATorchNG.model_probability`: the column-wise
+        ``softmax`` over models of :meth:`model_loglik`, i.e. ``P(model h |
+        x_t)``. Each column sums to 1; all ones for a single model.
+
+        Parameters
+        ----------
+        X : np.ndarray of shape (n_channels, n_samples)
+            Raw (unpreprocessed) data.
+
+        Returns
+        -------
+        prob : np.ndarray of shape (n_models, n_samples)
+
+        Raises
+        ------
+        ValueError
+            If the model is unfitted, if ``X`` is non-finite, or if every model
+            underflows to ``-inf`` log-likelihood at some sample.
+        RuntimeError
+            If the fit ended degenerate (issue #50).
+        """
+        self._check_usable("compute the model probability")
+        assert self.model_ is not None
+
+        return self.model_.model_probability(X)
+
     def variance_order(
         self, model_idx: int = 0, return_svar: bool = False
     ) -> Union[np.ndarray, tuple]: