diff --git a/docs/api/mne-compat.md b/docs/api/mne-compat.md index c4afba5..6f240b7 100644 --- a/docs/api/mne-compat.md +++ b/docs/api/mne-compat.md @@ -86,4 +86,27 @@ 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. +## Inspecting pamica-specific metadata + +An `mne.preprocessing.ICA` has no field for AMICA's adaptive source densities or +component sharing, so rather than drop them, the wrapper exposes them directly: + +```python +from pamica.mne_compat import AMICAICA, PDFTYPE_NAMES + +ica = AMICAICA(n_models=2, random_state=42).fit(raw, max_iter=100) + +families = ica.get_pdftype(model_idx=0) # (n_components,) codes 0-4 +names = [PDFTYPE_NAMES[c] for c in families] # e.g. "generalized_gaussian" +rho = ica.get_rho(model_idx=0) # (n_mix, n_components) GG shape +shared = ica.shared_components() # [(model, comp), ...] groups +``` + +`get_pdftype` returns each component's density family (0 generalized Gaussian, +1 super-Gaussian cosh, 2 Gaussian, 3 logistic, 4 sub-Gaussian cosh; they differ +per component only under the adaptive switcher `pdftype=1`). `get_rho` is the +generalized-Gaussian shape (meaningful for `pdftype=0`). `shared_components` +lists components merged across models by `share_comps` (empty otherwise). The +same three accessors exist on the scikit-learn-style [`AMICA`](amica.md). + ::: pamica.mne_compat.AMICAICA diff --git a/docs/changelog.md b/docs/changelog.md index 1669e00..f7caebe 100644 --- a/docs/changelog.md +++ b/docs/changelog.md @@ -42,6 +42,13 @@ MNE-Python compatibility layer (epic #139), additive: the scikit-learn-style 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). +- pamica-specific fitted metadata is inspectable through the MNE wrapper rather + than silently dropped by the `mne.preprocessing.ICA` export: `get_pdftype(model_idx=...)` + returns each component's source-density family code (0-4, named by + `pamica.mne_compat.PDFTYPE_NAMES`), `get_rho(model_idx=...)` the + generalized-Gaussian shape parameters, and `shared_components()` the components + merged across models by `share_comps`. The same accessors are added to + `AMICA`/`AMICATorchNG` (phase 3, #142). ## 0.2.2 diff --git a/pamica/amica.py b/pamica/amica.py index 6a9f66f..b264107 100644 --- a/pamica/amica.py +++ b/pamica/amica.py @@ -477,6 +477,47 @@ def model_probability(self, X: np.ndarray) -> np.ndarray: return self.model_.model_probability(X) + def get_pdftype(self, model_idx: int = 0) -> np.ndarray: + """Per-source density-family code for model ``model_idx`` (issue #142). + + Delegates to :meth:`AMICATorchNG.get_pdftype`. One integer per source + component (0-4; see :data:`pamica.torch_impl.PDFTYPE_NAMES`). + + Returns + ------- + np.ndarray of int, shape (n_sources,) + """ + self._check_usable("get the density family") + assert self.model_ is not None + + return self.model_.get_pdftype(model_idx=model_idx) + + def get_rho(self, model_idx: int = 0) -> np.ndarray: + """Generalized-Gaussian shape ``rho`` for model ``model_idx`` (issue #142). + + Delegates to :meth:`AMICATorchNG.get_rho`. + + Returns + ------- + np.ndarray of float, shape (n_mix, n_sources) + """ + self._check_usable("get rho") + assert self.model_ is not None + + return self.model_.get_rho(model_idx=model_idx) + + def shared_components(self) -> list: + """Components shared across models by ``share_comps`` (issue #142). + + Delegates to :meth:`AMICATorchNG.shared_components`: one group of + ``(model_idx, source_idx)`` pairs per shared column; empty when nothing + is shared. + """ + self._check_usable("get the shared components") + assert self.model_ is not None + + return self.model_.shared_components() + def variance_order( self, model_idx: int = 0, return_svar: bool = False ) -> Union[np.ndarray, tuple]: diff --git a/pamica/mne_compat/__init__.py b/pamica/mne_compat/__init__.py index 26a5ad9..2bc1ef3 100644 --- a/pamica/mne_compat/__init__.py +++ b/pamica/mne_compat/__init__.py @@ -10,6 +10,6 @@ pattern as the optional ``mlx`` backend). """ -from .core import AMICAICA +from .core import AMICAICA, PDFTYPE_NAMES -__all__ = ["AMICAICA"] +__all__ = ["AMICAICA", "PDFTYPE_NAMES"] diff --git a/pamica/mne_compat/core.py b/pamica/mne_compat/core.py index 05a7e67..6132b70 100644 --- a/pamica/mne_compat/core.py +++ b/pamica/mne_compat/core.py @@ -26,6 +26,9 @@ from mne.preprocessing import ICA as _MNEICA # ty: ignore[unresolved-import] from ..amica import AMICA +from ..torch_impl import PDFTYPE_NAMES + +__all__ = ["AMICAICA", "PDFTYPE_NAMES"] class AMICAICA: @@ -415,6 +418,48 @@ def plot_model_probability(self, inst, *, srate: Optional[float] = None, **kwarg srate = float(self.info_["sfreq"]) return _plot_model_probability(lht=lht, srate=srate, **kwargs) + # ------------------------------------------------------------------ + # pamica-specific metadata (issue #142) + # + # MNE's ``ICA`` carries no source-density family, GG shape, or + # component-sharing state, so these are exposed here rather than silently + # dropped by the ``mne.preprocessing.ICA`` export. + # ------------------------------------------------------------------ + def get_pdftype(self, *, model_idx: int = 0) -> np.ndarray: + """Per-component source-density family code for model ``model_idx``. + + One integer per ICA component (0-4); map to names with + :data:`pamica.mne_compat.PDFTYPE_NAMES`. All components share one family + unless the adaptive switcher (``pdftype=1``) moved them (issue #26). + """ + self._check_fitted("get the density family") + model_idx = self._check_model_idx(model_idx) + assert self.amica_ is not None + return self.amica_.get_pdftype(model_idx=model_idx) + + def get_rho(self, *, model_idx: int = 0) -> np.ndarray: + """Generalized-Gaussian shape ``rho`` for model ``model_idx``. + + Shape ``(n_mix, n_components)``; ``rho == 2`` is Gaussian-shaped, + ``rho == 1`` Laplacian, ``rho < 1`` heavier-tailed. Meaningful only for + the generalized-Gaussian family (``pdftype=0``). + """ + self._check_fitted("get rho") + model_idx = self._check_model_idx(model_idx) + assert self.amica_ is not None + return self.amica_.get_rho(model_idx=model_idx) + + def shared_components(self) -> list: + """Components shared across models by ``share_comps`` (issue #60). + + One group of ``(model_idx, component_idx)`` pairs per shared column; + empty when nothing is shared (always so for a single model or a default + multi-model fit with ``share_comps`` off). + """ + self._check_fitted("get the shared components") + assert self.amica_ is not None + return self.amica_.shared_components() + # ------------------------------------------------------------------ # Helpers # ------------------------------------------------------------------ diff --git a/pamica/tests/mne_tests/test_mne_wrapper.py b/pamica/tests/mne_tests/test_mne_wrapper.py index da89323..902f935 100644 --- a/pamica/tests/mne_tests/test_mne_wrapper.py +++ b/pamica/tests/mne_tests/test_mne_wrapper.py @@ -18,7 +18,7 @@ mne = pytest.importorskip("mne") -from pamica.mne_compat import AMICAICA # noqa: E402 (after importorskip) +from pamica.mne_compat import AMICAICA, PDFTYPE_NAMES # noqa: E402 (after importorskip) mne.set_log_level("ERROR") @@ -434,3 +434,34 @@ def test_degenerate_fit_refused_on_dominance_methods(raw): ica.get_model_probability(raw) with pytest.raises(RuntimeError, match="degenerate"): ica.plot_model_probability(raw) + + +# --- pamica-specific metadata (issue #142) ---------------------------------- +def test_get_pdftype_and_rho_per_model(fitted_2m): + for h in range(2): + pdf = fitted_2m.get_pdftype(model_idx=h) + assert pdf.shape == (fitted_2m.n_components_,) + assert np.array_equal(np.unique(pdf), [0]) # default generalized Gaussian + assert PDFTYPE_NAMES[int(pdf[0])] == "generalized_gaussian" + rho = fitted_2m.get_rho(model_idx=h) + assert rho.shape == (fitted_2m.n_mix, fitted_2m.n_components_) + assert np.all(np.isfinite(rho)) + + +def test_shared_components_empty_by_default(fitted_2m): + assert fitted_2m.shared_components() == [] + + +def test_metadata_respects_model_idx_bounds(fitted_2m): + with pytest.raises(ValueError, match="out of range"): + fitted_2m.get_pdftype(model_idx=2) + with pytest.raises(ValueError, match="out of range"): + fitted_2m.get_rho(model_idx=5) + + +def test_metadata_requires_fit(): + ica = AMICAICA() + with pytest.raises(ValueError, match="must be fitted"): + ica.get_pdftype() + with pytest.raises(ValueError, match="must be fitted"): + ica.shared_components() diff --git a/pamica/tests/torch_tests/test_amica_ng_wrapper.py b/pamica/tests/torch_tests/test_amica_ng_wrapper.py index 4d4cb9f..0e30e30 100644 --- a/pamica/tests/torch_tests/test_amica_ng_wrapper.py +++ b/pamica/tests/torch_tests/test_amica_ng_wrapper.py @@ -643,6 +643,12 @@ def test_unfitted_output_raises_not_fitted(): model.mir(np.zeros((NW, 16))) with pytest.raises(ValueError, match="fitted"): model.pmi(np.zeros((NW, 16))) + with pytest.raises(ValueError, match="fitted"): + model.get_pdftype() + with pytest.raises(ValueError, match="fitted"): + model.get_rho() + with pytest.raises(ValueError, match="fitted"): + model.shared_components() def test_degenerate_fit_refuses_output(real_data, tmp_path, caplog): @@ -675,6 +681,9 @@ def test_degenerate_fit_refuses_output(real_data, tmp_path, caplog): lambda: model.variance_order(), lambda: model.mir(real_data[:, :512]), lambda: model.pmi(real_data[:, :512]), + lambda: model.get_pdftype(), + lambda: model.get_rho(), + lambda: model.shared_components(), ): with pytest.raises(RuntimeError, match="degenerate.*nan_ll"): action() diff --git a/pamica/tests/torch_tests/test_ng_metadata.py b/pamica/tests/torch_tests/test_ng_metadata.py new file mode 100644 index 0000000..968c14a --- /dev/null +++ b/pamica/tests/torch_tests/test_ng_metadata.py @@ -0,0 +1,145 @@ +"""Fitted-parameter metadata accessors (issue #142). + +``get_pdftype``/``get_rho``/``shared_components`` expose the source-density +family, GG shape, and component-sharing state that the EEGLAB/MNE exports do not +carry. Real sample EEG only (no synthetic data). +""" + +from pathlib import Path + +import numpy as np +import pytest + +from pamica import AMICATorchNG +from pamica.torch_impl import PDFTYPE_NAMES +from pamica.torch_impl.utils import load_eeglab_data + +SAMPLE_DIR = Path(__file__).resolve().parents[2] / "sample_data" +DATA_FILE = SAMPLE_DIR / "eeglab_data.fdt" +NW = 32 +FIELD = 30504 +NMIX = 3 +SEED = 42 + +pytestmark = pytest.mark.skipif(not DATA_FILE.exists(), reason="sample data missing") + + +@pytest.fixture(scope="module") +def real_data(): + return load_eeglab_data(str(DATA_FILE), data_dim=NW, field_dim=FIELD).astype( + np.float64 + ) + + +def _fit(X, n_models=1, max_iter=8, **kw): + m = AMICATorchNG( + n_channels=NW, n_models=n_models, n_mix=NMIX, seed=SEED, device="cpu", **kw + ) + m.fit(X, max_iter=max_iter, verbose=False) + return m + + +def test_pdftype_names_cover_all_codes(): + assert set(PDFTYPE_NAMES) == {0, 1, 2, 3, 4} + + +def test_get_pdftype_default_is_generalized_gaussian(real_data): + m = _fit(real_data) + pdf = m.get_pdftype() + assert pdf.shape == (NW,) + assert np.array_equal(np.unique(pdf), [0]) # pdftype=0 default + assert PDFTYPE_NAMES[int(pdf[0])] == "generalized_gaussian" + + +def test_get_pdftype_reflects_fixed_family(real_data): + m = _fit(real_data, pdftype=2) # Gaussian + assert np.array_equal(np.unique(m.get_pdftype()), [2]) + + +def test_get_rho_shape_and_bounds(real_data): + m = _fit(real_data) + rho = m.get_rho() + assert rho.shape == (NMIX, NW) + assert np.all(rho >= m.minrho) and np.all(rho <= m.maxrho) + + +def test_shared_components_empty_without_sharing(real_data): + assert _fit(real_data, n_models=1).shared_components() == [] + assert _fit(real_data, n_models=2).shared_components() == [] # share off default + + +def test_shared_components_reports_partial_merges(real_data): + """A PARTIAL merge (some columns single-model) exercises both the grouping + and the exclude-filter: shared groups appear, non-shared columns do not.""" + m = _fit( + real_data[:, :4096], + n_models=2, + max_iter=25, + block_size=1024, + share_comps=True, + share_start=8, + share_iter=10, + do_newton=True, + ) + groups = m.shared_components() + # Partial: not zero (something merged), and fewer than n_sources (some + # columns stay single-model and are correctly excluded from the groups). + assert 0 < len(groups) < NW + for group in groups: + models = {h for h, _ in group} + assert len(models) >= 2 # a shared column spans >= 2 models + assert all(0 <= h < 2 and 0 <= i < NW for h, i in group) + + +def test_get_rho_differs_per_model(real_data): + """get_rho indexes comp_list per model, so a 2-model fit's rho differs.""" + m = _fit(real_data, n_models=2, max_iter=10) + assert not np.allclose(m.get_rho(model_idx=0), m.get_rho(model_idx=1)) + + +def test_get_pdftype_adaptive_switcher_is_mixed(real_data): + """The adaptive switcher (pdftype=1) moves sources into families {1, 4}; + get_pdftype must surface those per-source codes, not a uniform array.""" + m = AMICATorchNG( + n_channels=NW, + n_models=1, + n_mix=1, # adaptive mode is single-component + seed=SEED, + device="cpu", + pdftype=1, + kurt_start=3, + num_kurt=5, + kurt_int=1, + ) + m.fit(real_data, max_iter=20, verbose=False) + codes = set(np.unique(m.get_pdftype()).tolist()) + assert codes.issubset({1, 4}) and codes # only the two cosh families + + +def test_metadata_accessors_reject_bad_model_idx(real_data): + m = _fit(real_data, n_models=2, max_iter=5) + for call in (m.get_pdftype, m.get_rho): + with pytest.raises(ValueError, match="out of range"): + call(model_idx=2) + with pytest.raises(ValueError, match="out of range"): + call(model_idx=-1) # negative would silently wrap without the guard + + +def test_metadata_accessors_require_fit(): + m = AMICATorchNG(n_channels=NW, n_models=1, n_mix=NMIX, seed=SEED, device="cpu") + for call in (m.get_pdftype, m.get_rho, m.shared_components): + with pytest.raises(RuntimeError, match="fitted"): + call() + + +def test_amica_wrapper_delegates_metadata(real_data): + """The scikit-learn-style AMICA wrapper exposes the same accessors.""" + from pamica import AMICA + + a = AMICA(n_models=2, n_mix=NMIX, device="cpu", verbose=False) + a.fit(real_data, max_iter=8, seed=SEED) + assert a.get_pdftype(model_idx=1).shape == (NW,) + assert a.get_rho(model_idx=1).shape == (NMIX, NW) + assert a.shared_components() == [] # share off + with pytest.raises(ValueError, match="out of range"): + a.get_pdftype(model_idx=-1) # guard propagates through the wrapper diff --git a/pamica/torch_impl/__init__.py b/pamica/torch_impl/__init__.py index 065c84a..09f6512 100644 --- a/pamica/torch_impl/__init__.py +++ b/pamica/torch_impl/__init__.py @@ -5,11 +5,12 @@ with support for CUDA, ROCm, and Apple Silicon (MPS) backends. """ -from .core import AMICATorchNG +from .core import AMICATorchNG, PDFTYPE_NAMES from .utils import setup_device, check_numerical_stability __all__ = [ "AMICATorchNG", + "PDFTYPE_NAMES", "setup_device", "check_numerical_stability", ] diff --git a/pamica/torch_impl/core.py b/pamica/torch_impl/core.py index ff80d3f..8f4f89d 100644 --- a/pamica/torch_impl/core.py +++ b/pamica/torch_impl/core.py @@ -58,6 +58,17 @@ logger = logging.getLogger(__name__) +# Human-readable names for the ``pdftype``/``pdtype`` source-density family codes +# (issue #26; amica15.f90). Exposed alongside the numeric codes so a fitted +# model's per-source density family is inspectable (issue #142). +PDFTYPE_NAMES = { + 0: "generalized_gaussian", + 1: "super_gaussian_cosh", + 2: "gaussian", + 3: "logistic", + 4: "sub_gaussian_cosh", +} + _LOG2 = math.log(2.0) _LOG4 = math.log(4.0) # logistic-family normalizer (amica15.f90:1328) _HALF_LOG_PI = 0.5 * math.log(math.pi) @@ -1953,6 +1964,23 @@ def _reject_outliers(self, ll_vec: torch.Tensor): int(self.good_idx.numel()), ) + def _check_model_idx(self, model_idx: int) -> None: + """Validate a model index against the fitted ``n_models``. + + Raises a clear ``ValueError`` (rejecting negatives, which torch's + negative indexing would otherwise turn into a silent wrong-model result) + instead of an opaque tensor ``IndexError``. + """ + if not isinstance(model_idx, (int, np.integer)): + raise TypeError( + f"model_idx must be an int, got {type(model_idx).__name__}." + ) + if not (0 <= model_idx < self.n_models): + raise ValueError( + f"model_idx={model_idx} out of range for a {self.n_models}-model " + f"fit (valid: 0..{self.n_models - 1})." + ) + def transform(self, X: np.ndarray, model_idx: int = 0) -> np.ndarray: """Apply the learned unmixing matrix to (new) data. @@ -1966,6 +1994,7 @@ def transform(self, X: np.ndarray, model_idx: int = 0) -> np.ndarray: raise RuntimeError( "AMICATorchNG.transform() requires a fitted model; call fit() first." ) + self._check_model_idx(model_idx) X_t = torch.from_numpy(np.ascontiguousarray(X)).to(self.device, self.dtype) X_t = self.sphere @ (X_t - self.mean) # c is the per-model data-space center: unmix as W(x - c) (issue #27). @@ -1979,6 +2008,7 @@ def get_mixing_matrix(self, model_idx: int = 0) -> np.ndarray: "AMICATorchNG.get_mixing_matrix() requires a fitted model; call " "fit() first." ) + self._check_model_idx(model_idx) return self.A[:, self.comp_list[:, model_idx]].T.cpu().numpy() def get_unmixing_matrix(self, model_idx: int = 0) -> np.ndarray: @@ -1988,6 +2018,7 @@ def get_unmixing_matrix(self, model_idx: int = 0) -> np.ndarray: "AMICATorchNG.get_unmixing_matrix() requires a fitted model; call " "fit() first." ) + self._check_model_idx(model_idx) return self.W[:, :, model_idx].T.cpu().numpy() def _pca_reduced(self) -> bool: @@ -2031,6 +2062,7 @@ def mir( raise RuntimeError( "AMICATorchNG.mir() requires a fitted model; call fit() first." ) + self._check_model_idx(model_idx) if self._pca_reduced(): raise ValueError( "mir() is incompatible with PCA reduction (pcakeep/pcadb): " @@ -2165,6 +2197,92 @@ def model_probability(self, X: np.ndarray) -> np.ndarray: ex = np.exp(Lht - col_max) return ex / ex.sum(axis=0, keepdims=True) + # ------------------------------------------------------------------ + # Fitted-parameter metadata (issue #142) + # ------------------------------------------------------------------ + def get_pdftype(self, model_idx: int = 0) -> np.ndarray: + """Per-source density-family code for model ``model_idx``. + + One integer per source component (0-4; see + :data:`pamica.torch_impl.PDFTYPE_NAMES`): 0 generalized Gaussian, 1 + super-Gaussian cosh, 2 Gaussian, 3 logistic, 4 sub-Gaussian cosh. All + sources share ``pdftype`` unless the adaptive switcher (``pdftype=1``) + moved them individually (issue #26). + + Returns + ------- + np.ndarray of int, shape (n_sources,) + """ + if self.pdtype is None: + raise RuntimeError( + "AMICATorchNG.get_pdftype() requires a fitted model; call fit() first." + ) + self._check_model_idx(model_idx) + return self.pdtype[:, model_idx].detach().cpu().numpy() + + def get_rho(self, model_idx: int = 0) -> np.ndarray: + """Generalized-Gaussian shape parameter ``rho`` for model ``model_idx``. + + One value per (mixture component, source): ``rho == 2`` is Gaussian- + shaped, ``rho == 1`` Laplacian, ``rho < 1`` heavier-tailed. Only the + generalized-Gaussian family (``pdftype=0``) updates ``rho``; for every + non-zero code (1-4, the fixed and adaptive cosh families) it stays frozen + at ``rho0`` and does not describe the fitted density. + + Returns + ------- + np.ndarray of float, shape (n_mix, n_sources) + """ + if self.rho is None or self.comp_list is None: + raise RuntimeError( + "AMICATorchNG.get_rho() requires a fitted model; call fit() first." + ) + self._check_model_idx(model_idx) + # Defense-in-depth, matching state_dict()'s isfinite sweep: a degenerate + # multi-model fit can leave one model's rho non-finite without the + # aggregate LL tripping nan_ll, and _check_usable only inspects + # stop_reason. Refuse rather than return a silent NaN. + if not torch.isfinite(self.rho).all(): + raise RuntimeError( + "AMICATorchNG.get_rho(): rho holds non-finite values (a " + "degenerate fit); inspect stop_reason and refit." + ) + idx = self.comp_list[:, model_idx] + return self.rho[:, idx].detach().cpu().numpy() + + def shared_components(self) -> list: + """Components shared across models by ``share_comps`` (issue #60). + + ``share_comps`` folds near-collinear components of different models onto + one shared mixing column + density, recorded as a repeated index in + ``comp_list``. Returns one group per shared column: a list of + ``(model_idx, source_idx)`` pairs that all reference it. Empty when no + component is shared across two or more models (always for one model, and + for a default multi-model fit with ``share_comps`` off). + + Note that a merge synchronizes only the mixture parameters routed through + ``comp_list`` (``mu``/``alpha``/``beta``/``rho``); the per-source density + *family* code ``pdtype`` is a separate tensor and is not synchronized, so + under the adaptive switcher (``pdftype=1``) a shared pair can still report + different :meth:`get_pdftype` codes. + + Returns + ------- + list of list of tuple(int, int) + """ + if self.comp_list is None: + raise RuntimeError( + "AMICATorchNG.shared_components() requires a fitted model; call " + "fit() first." + ) + cl = self.comp_list.detach().cpu().numpy() # (n_sources, n_models) + groups = [] + for col in np.unique(cl): + src, mdl = np.where(cl == col) + if np.unique(mdl).size >= 2: + groups.append([(int(h), int(i)) for i, h in zip(src, mdl)]) + return groups + # ------------------------------------------------------------------ # EEGLAB drop-in output (issue #92) # ------------------------------------------------------------------ @@ -2210,6 +2328,7 @@ def variance_order( "AMICATorchNG.variance_order() requires a fitted model; call " "fit() first." ) + self._check_model_idx(model_idx) cl = self.comp_list[:, model_idx].cpu().numpy() alpha = self.alpha[:, cl].cpu().numpy() mu = self.mu[:, cl].cpu().numpy()