AMICA

The main scikit-learn-style interface. It wraps a natural-gradient expectation-maximization (EM) backend: AMICATorchNG by default (float64, Fortran parity), or the Apple-GPU AMICAMLXNG with backend="mlx" (float32, issue #313). Every method works the same on both; see Selecting a backend for the rules and the precision of each.

from pamica import AMICA

model = AMICA(backend="mlx").fit(X, pcakeep=X.shape[0] - 1, seed=42)
model.save("model.pt")          # records the backend
model = AMICA.load("model.pt")  # restored on the MLX backend

device and a dtype fit keyword apply to the PyTorch backend only and raise ValueError with backend="mlx"; backend="mlx" without MLX installed raises ImportError.

fit's max_iter, lrate, do_mean, do_sphere and do_newton default to the selected backend's own values (issue #354), so AMICA().fit(X) runs the same fit as the backend class with its defaults: lrate=0.1, the compiled amica15 default, where EEGLAB's runamica15.m uses 0.05. Default settings compares every default with the compiled binary's and EEGLAB's.

get_sphere(), get_mean() and get_model_center(model_idx) return the fitted preprocessing as float64 arrays on either backend, so the transform can be composed by hand: transform(X) is W @ (sphere @ (X - mean[:, None]) - c[:, None]) with W = get_unmixing_matrix(). get_sensor_mixing_matrix() gives the scalp maps in input-channel space, which is the only valid back-map after rank reduction.

save writes format_version 2, which records the backend; load restores the model on that backend and still reads version 1 files, which predate backend selection and always hold a PyTorch model. A model saved before issue #334, which changed how the mixing matrix is stored, is converted on load without loss; one in which share_comps had merged components raises ValueError and must be refit.

pamica.AMICA

Adaptive Mixture ICA over pamica's natural-gradient EM backends.

This is the main interface for pamica, providing a scikit-learn style API over :class:AMICATorchNG, the natural-gradient EM implementation that matches the Fortran reference (Newton, exact-EM mixture updates, symmetric-ZCA sphere, Jacobian LL), or, with backend="mlx", over its Apple-GPU port :class:pamica.mlx_impl.AMICAMLXNG (issue #313). Every method below behaves the same on both backends; only precision differs.

Parameters:

Name	Type	Description	Default
n_models	int	Number of ICA models to learn	1
n_mix	int	Number of mixture components per source	3
device	str or device	Device to use ('cuda', 'mps', 'cpu', or None for auto), passed to :class:AMICATorchNG. With None (auto), the backend moves an auto-selected MPS device to the CPU with a logged warning, because it computes in float64 for Fortran parity and MPS cannot represent it; pass dtype=torch.float32 (with device="mps") to run on MPS instead. PyTorch backend only: with backend="mlx" it must stay None.	None
verbose	bool	Whether to show progress during fitting	True
backend	(torch, mlx)	Which backend :meth:fit builds. "torch" is :class:AMICATorchNG, float64 by default and the Fortran-parity path. "mlx" is :class:pamica.mlx_impl.AMICAMLXNG, the fastest option on Apple Silicon, which computes in float32 only (about 7 significant digits, not float64 parity), so it takes neither device nor dtype. MLX is an optional dependency: "mlx" without it raises ImportError here, and import pamica never imports it.	"torch"

Raises:

Type	Description
ValueError	If backend is not "torch" or "mlx", or device is set with backend="mlx".
ImportError	If backend="mlx" and MLX is not installed.

Attributes:

Name	Type	Description
model_	AMICATorchNG or AMICAMLXNG	The underlying backend model
is_fitted_	bool	Whether a usable model is available. fit sets this True only when the fit ended on a non-degenerate stop (a convergence stop or max_iter); a degenerate fit (see converged_) leaves it False, and every output method (transform, the get_* accessors, write_amica_output, save) refuses such a model (issue #50).
converged_	bool	Whether the last fit ended on a usable stop rather than a degenerate one (stop_reason_ not in the backend's _DEGENERATE_STOP_REASONS: nan_ll/singular_ll/nan_direction/nan_params, or restart_error under best-of-N restarts). A max_iter stop counts as usable. A degenerate fit holds non-finite parameters, or stopped before applying a non-finite step, and would produce NaN sources (issue #50).
stop_reason_	str or None	Why the last fit stopped (the backend stop_reason): "max_iter", "lrate_floor", "grad_norm_floor", "min_dll", "grad_norm", "nan_ll", "singular_ll", "nan_direction" or "nan_params". "lrate_floor" to "grad_norm" are Fortran-faithful convergence stops (issue #207: lrate_floor/ grad_norm_floor fire together as two halves of the same likelihood-decrease branch; min_dll/grad_norm are separate, unconditional per-iteration checks). The last four are degenerate (see converged_): a non-finite log-likelihood (nan_ll/ singular_ll), a non-finite update direction caught before it is applied (nan_direction), or non-finite parameters right after an update (nan_params); PyTorch and MLX use the same set (issue #339 review), and under best-of-N restarts a restart that raised is recorded as "restart_error", also degenerate. None of these checks short- circuits on an earlier one in the same iteration, so under the shipped use_grad_norm=True default "grad_norm" always takes precedence over "grad_norm_floor" when both would apply -- "grad_norm_floor" only surfaces as this value when use_grad_norm=False (see AMICATorchNG's use_grad_norm docstring for the full explanation).
ll_history_	list	Log-likelihood history during training (the true per-iteration trajectory; may dip below its peak on a late overshoot): entry i is the likelihood of the parameters iteration i started from, recorded before that iteration's checks and update.
final_ll_	float	Log-likelihood of the fitted parameters (issue #51). Use this, not ll_history_[-1], as the model's log-likelihood: with the best-iterate safeguard the returned parameters can be an earlier, higher-LL iterate. Exact after a convergence stop, which exits before that iteration's update, as the reference does (issue #339); after max_iter it is the likelihood one update before the returned parameters, also as in the reference.
mir_history_	list	Mutual Information Reduction (MIR) waypoint trajectory (issue #137), populated when fit is called with mir_step > 0: (iteration, mir_nats, variance) tuples from the mid-fit W/sphere. Like ll_history_, a keep_best restore does not rewrite it -- use :meth:mir on the fitted model for the value of the returned parameters, not mir_history_[-1]. Not index-aligned with ll_history_: entry i is computed after iteration i's update, while ll_history_[i] is the likelihood of the parameters before it, so the two are one update apart (issue #161). An iteration that ends the fit on a stop (a convergence check or a degenerate value) takes no update and so records no waypoint: the last waypoint then belongs to the iteration before it (issue #339).
restart_seeds_	list	The seed each restart ran from (issue #198). One entry for a default n_restarts=1 fit, n_restarts entries otherwise.
restart_lls_	list	Each restart's returned log-likelihood, index-aligned with restart_seeds_; NaN for a restart that ended degenerate (those are recorded but excluded from the selection). final_ll_ is the winner.
restart_stop_reasons_	list	Each restart's stop_reason, index-aligned with restart_seeds_.

Examples:

>>> from pamica import AMICA
>>> import numpy as np
>>>
>>> # Generate sample data
>>> X = np.random.randn(32, 10000)  # 32 channels, 10000 samples
>>>
>>> # Fit AMICA model
>>> amica = AMICA(n_models=1, n_mix=3)
>>> amica.fit(X, max_iter=100)
>>>
>>> # Transform data to sources
>>> S = amica.transform(X)
>>>
>>> # Get mixing matrix
>>> A = amica.get_mixing_matrix()
>>>
>>> # The same on the Apple GPU (requires the mlx extra)
>>> S_mlx = AMICA(backend="mlx").fit(X, max_iter=100).transform(X)

Source code in pamica/amica.py

class AMICA:
    """
    Adaptive Mixture ICA over pamica's natural-gradient EM backends.

    This is the main interface for pamica, providing a scikit-learn style
    API over :class:`AMICATorchNG`, the natural-gradient EM implementation
    that matches the Fortran reference (Newton, exact-EM mixture updates,
    symmetric-ZCA sphere, Jacobian LL), or, with ``backend="mlx"``, over its
    Apple-GPU port :class:`pamica.mlx_impl.AMICAMLXNG` (issue #313). Every
    method below behaves the same on both backends; only precision differs.

    Parameters
    ----------
    n_models : int, default=1
        Number of ICA models to learn
    n_mix : int, default=3
        Number of mixture components per source
    device : str or torch.device, optional
        Device to use ('cuda', 'mps', 'cpu', or None for auto), passed to
        :class:`AMICATorchNG`. With ``None`` (auto), the backend moves an
        auto-selected MPS device to the CPU with a logged warning, because it
        computes in float64 for Fortran parity and MPS cannot represent it;
        pass ``dtype=torch.float32`` (with ``device="mps"``) to run on MPS
        instead. PyTorch backend only: with ``backend="mlx"`` it must stay
        ``None``.
    verbose : bool, default=True
        Whether to show progress during fitting
    backend : {"torch", "mlx"}, default="torch"
        Which backend :meth:`fit` builds. ``"torch"`` is
        :class:`AMICATorchNG`, float64 by default and the Fortran-parity
        path. ``"mlx"`` is :class:`pamica.mlx_impl.AMICAMLXNG`, the fastest
        option on Apple Silicon, which computes in float32 only (about 7
        significant digits, not float64 parity), so it takes neither
        ``device`` nor ``dtype``. MLX is an optional dependency: ``"mlx"``
        without it raises ``ImportError`` here, and ``import pamica`` never
        imports it.

    Raises
    ------
    ValueError
        If ``backend`` is not ``"torch"`` or ``"mlx"``, or ``device`` is set
        with ``backend="mlx"``.
    ImportError
        If ``backend="mlx"`` and MLX is not installed.

    Attributes
    ----------
    model_ : AMICATorchNG or AMICAMLXNG
        The underlying backend model
    is_fitted_ : bool
        Whether a *usable* model is available. ``fit`` sets this True only when
        the fit ended on a non-degenerate stop (a convergence stop or
        ``max_iter``); a degenerate fit (see ``converged_``) leaves it False,
        and every output method (``transform``, the ``get_*`` accessors,
        ``write_amica_output``, ``save``) refuses such a model (issue #50).
    converged_ : bool
        Whether the last ``fit`` ended on a usable stop rather than a degenerate
        one (``stop_reason_`` not in the backend's ``_DEGENERATE_STOP_REASONS``:
        ``nan_ll``/``singular_ll``/``nan_direction``/``nan_params``, or
        ``restart_error`` under best-of-N restarts). A ``max_iter`` stop counts
        as usable. A degenerate fit holds non-finite parameters, or stopped
        before applying a non-finite step, and would produce NaN sources
        (issue #50).
    stop_reason_ : str or None
        Why the last ``fit`` stopped (the backend ``stop_reason``):
        ``"max_iter"``, ``"lrate_floor"``, ``"grad_norm_floor"``, ``"min_dll"``,
        ``"grad_norm"``, ``"nan_ll"``, ``"singular_ll"``, ``"nan_direction"``
        or ``"nan_params"``. ``"lrate_floor"`` to ``"grad_norm"`` are
        Fortran-faithful convergence stops (issue #207: ``lrate_floor``/
        ``grad_norm_floor`` fire together as two halves of the same
        likelihood-decrease branch; ``min_dll``/``grad_norm`` are separate,
        unconditional per-iteration checks). The last four are degenerate (see
        ``converged_``): a non-finite log-likelihood (``nan_ll``/
        ``singular_ll``), a non-finite update direction caught before it is
        applied (``nan_direction``), or non-finite parameters right after an
        update (``nan_params``); PyTorch and MLX use the same set (issue #339
        review), and under best-of-N restarts a restart that raised is
        recorded as ``"restart_error"``, also degenerate. None of these checks short-
        circuits on an earlier one in the same iteration, so under the
        shipped ``use_grad_norm=True`` default ``"grad_norm"`` always takes
        precedence over ``"grad_norm_floor"`` when both would apply --
        ``"grad_norm_floor"`` only surfaces as this value when
        ``use_grad_norm=False`` (see ``AMICATorchNG``'s ``use_grad_norm``
        docstring for the full explanation).
    ll_history_ : list
        Log-likelihood history during training (the true per-iteration
        trajectory; may dip below its peak on a late overshoot): entry ``i`` is
        the likelihood of the parameters iteration ``i`` started from, recorded
        before that iteration's checks and update.
    final_ll_ : float
        Log-likelihood of the *fitted* parameters (issue #51). Use this, not
        ``ll_history_[-1]``, as the model's log-likelihood: with the best-iterate
        safeguard the returned parameters can be an earlier, higher-LL iterate.
        Exact after a convergence stop, which exits before that iteration's
        update, as the reference does (issue #339); after ``max_iter`` it is
        the likelihood one update before the returned parameters, also as in
        the reference.
    mir_history_ : list
        Mutual Information Reduction (MIR) waypoint trajectory (issue #137),
        populated when ``fit`` is called
        with ``mir_step > 0``: ``(iteration, mir_nats, variance)`` tuples from
        the mid-fit ``W``/``sphere``. Like ``ll_history_``, a ``keep_best``
        restore does not rewrite it -- use :meth:`mir` on the fitted model for
        the value of the *returned* parameters, not ``mir_history_[-1]``.
        Not index-aligned with ``ll_history_``: entry ``i`` is computed after
        iteration ``i``'s update, while ``ll_history_[i]`` is the likelihood of
        the parameters before it, so the two are one update apart (issue #161).
        An iteration that ends the fit on a stop (a convergence check or a
        degenerate value) takes no update and so records no waypoint: the last
        waypoint then belongs to the iteration before it (issue #339).
    restart_seeds_ : list
        The seed each restart ran from (issue #198). One entry for a default
        ``n_restarts=1`` fit, ``n_restarts`` entries otherwise.
    restart_lls_ : list
        Each restart's returned log-likelihood, index-aligned with
        ``restart_seeds_``; NaN for a restart that ended degenerate (those are
        recorded but excluded from the selection). ``final_ll_`` is the winner.
    restart_stop_reasons_ : list
        Each restart's ``stop_reason``, index-aligned with ``restart_seeds_``.

    Examples
    --------
    >>> from pamica import AMICA
    >>> import numpy as np
    >>>
    >>> # Generate sample data
    >>> X = np.random.randn(32, 10000)  # 32 channels, 10000 samples
    >>>
    >>> # Fit AMICA model
    >>> amica = AMICA(n_models=1, n_mix=3)
    >>> amica.fit(X, max_iter=100)
    >>>
    >>> # Transform data to sources
    >>> S = amica.transform(X)
    >>>
    >>> # Get mixing matrix
    >>> A = amica.get_mixing_matrix()
    >>>
    >>> # The same on the Apple GPU (requires the mlx extra)
    >>> S_mlx = AMICA(backend="mlx").fit(X, max_iter=100).transform(X)
    """

    def __init__(
        self,
        n_models: int = 1,
        n_mix: int = 3,
        device: Optional[Union[str, torch.device]] = None,
        verbose: bool = True,
        backend: str = "torch",
    ):
        _check_backend(backend, device)
        self.n_models = n_models
        self.n_mix = n_mix
        self.device = device
        self.verbose = verbose
        self.backend = backend

        self.model_: Optional[Union[AMICATorchNG, "AMICAMLXNG"]] = None
        self.is_fitted_ = False
        self.ll_history_ = []
        self.final_ll_ = None
        self.stop_reason_ = None
        self.converged_ = False
        self.mir_history_ = []
        # Best-of-N restart records (issue #198), mirrored off the backend by
        # fit()/load(): index-aligned lists of the seed each restart ran from,
        # the log-likelihood it returned (NaN for a degenerate restart) and why
        # it stopped. One entry even for the default single-restart fit.
        self.restart_seeds_ = []
        self.restart_lls_ = []
        self.restart_stop_reasons_ = []
        # Set by from_params_file (issue #132 review item 2): the full
        # translated parameter-file dict, applied by fit() as per-call
        # defaults (an explicitly passed fit()/backend kwarg always wins).
        # None for an instance built directly via AMICA(...).
        self._file_params: Optional[dict] = None

    def fit(
        self,
        X: np.ndarray,
        max_iter=_UNSET,
        lrate=_UNSET,
        do_mean=_UNSET,
        do_sphere=_UNSET,
        do_newton=_UNSET,
        mir_step: int = 0,
        **kwargs,
    ) -> "AMICA":
        """
        Fit AMICA model to data.

        The defaults of ``max_iter``, ``lrate``, ``do_mean``, ``do_sphere``
        and ``do_newton`` are the selected backend's own, read from its
        signatures (issue #354), so ``AMICA().fit(X)`` fits exactly as the
        backend class does with its defaults. Both backends default to the
        values listed below.

        Parameters
        ----------
        X : np.ndarray
            Input data of shape (n_channels, n_samples)
        max_iter : int, default=100
            Maximum number of iterations
        lrate : float, default=0.1
            Initial and maximum natural-gradient learning rate, the compiled
            amica15 default. EEGLAB's ``runamica15.m`` uses 0.05; see the
            defaults table in ``docs/guides/amica-differences.md``.
        do_mean : bool, default=True
            Whether to remove mean from data
        do_sphere : bool, default=True
            Whether to sphere (whiten) the data
        do_newton : bool, default=False
            Whether to enable the Fortran-parity Newton preconditioner (tune
            via ``newt_start``/``newtrate`` in ``**kwargs``).
        mir_step : int, default=0
            If > 0, compute MIR every ``mir_step`` iterations during training
            and record it in ``mir_history_`` (issue #137). ``0`` (default)
            disables the waypoints; see :meth:`AMICATorchNG.fit` for details
            and the interaction with ``keep_best``.
        **kwargs
            Additional parameters passed to the backend constructor,
            :class:`AMICATorchNG` or :class:`pamica.mlx_impl.AMICAMLXNG`
            (the two take the same keywords, except that ``dtype`` is
            PyTorch-only and raises ``ValueError`` with ``backend="mlx"``,
            which is float32-only; e.g. ``block_size``, ``rho0``, ``seed``,
            ``pcakeep``, ``use_min_dll``, ``min_dll``, ``maxincs``, ``use_grad_norm``,
            ``min_nd`` -- the issue #207 convergence stops, Fortran-faithful
            defaults ``True``/``1e-9``/``5``/``True``/``1e-7`` -- or
            ``do_opt_block``/``blk_min``/``blk_max``/``blk_step``, the
            issue #232 block-size search, off by default) -- the backend's tunables
            are constructor arguments, not fit() kwargs. A keyword the
            selected backend's constructor does not take raises
            ``TypeError``.

            ``n_restarts`` (default 1) runs the fit from that many seeds and
            keeps the highest-likelihood one (issue #198), recording every
            restart in ``restart_seeds_``/``restart_lls_``/
            ``restart_stop_reasons_``; it needs a base ``seed`` (or explicit
            ``restart_seeds``) and costs proportionally more time, since
            restarts run serially. ``n_restarts=1`` is bit-identical to a fit
            that never heard of restarts.

            Rank-deficient input (Maxwell-filtered MEG, average-referenced or
            interpolated EEG) is handled by ``mineig``/``mineig_rel``
            (issue #223): the model is sized to the detected numerical rank and
            :meth:`get_sensor_mixing_matrix` maps components back to input
            channels. ``mineig`` is an absolute eigenvalue floor and so
            unit-dependent; pass ``mineig_rel`` for data far from unit scale.

            When the instance was built via :meth:`from_params_file`
            (issue #132), any of the parameters above -- named or in ``**kwargs`` --
            left unset here falls back to that file's translated value instead
            of the backend's default; an explicitly passed argument always
            wins over the file. Settings the file carries that match neither a
            named ``fit()`` parameter nor a constructor keyword of the
            selected backend (data-location metadata like ``files``/
            ``outdir``/``data_dim``, or a setting with no equivalent in that
            backend) are not applied and are named in a single
            ``logger.warning``.

        Returns
        -------
        self : AMICA
            Fitted model

        Raises
        ------
        ValueError
            If ``X`` is not 2D, or ``dtype``/``device`` is passed with
            ``backend="mlx"``.
        TypeError
            If a keyword is neither a ``fit()`` parameter nor a constructor
            keyword of the selected backend.
        """
        # Validate input
        if X.ndim != 2:
            raise ValueError(f"X must be 2D array, got shape {X.shape}")

        n_channels, n_samples = X.shape

        # Re-checked here, not only in __init__, so a backend/device attribute
        # changed after construction cannot slip past the validation.
        _check_backend(self.backend, self.device)
        backend_cls = _backend_class(self.backend)
        ctor_params = _ctor_params(backend_cls)
        # Both categories computed before either is raised (PR #347 review
        # item 4, the same fix as AMICA_NumPy's _reject_unknown_kwargs): a
        # first version of this check raised on torch_only alone, so
        # backend='mlx' with fit(dtype=..., blocksize=...) named only
        # 'dtype' and silently dropped 'blocksize' from the message.
        # unknown excludes torch_only so a torch-only name (which is not in
        # ctor_params for backend='mlx' either) is not double-counted.
        torch_only = (
            sorted(set(kwargs) & set(_TORCH_ONLY_PARAMS))
            if self.backend == "mlx"
            else []
        )
        unknown = sorted(set(kwargs) - ctor_params - set(torch_only))
        if torch_only and unknown:
            raise TypeError(
                f"AMICA.fit got invalid keyword argument(s): {torch_only} "
                "apply only to backend='torch' (the MLX backend computes in "
                "float32 only, on its default device; use backend='torch' "
                f"for float64 Fortran-parity runs), and {unknown} is "
                "neither a fit() parameter nor a constructor keyword of "
                f"{backend_cls.__name__}."
            )
        if torch_only:
            raise ValueError(
                f"{torch_only} apply only to backend='torch'. The MLX "
                "backend computes in float32 only (Apple GPUs have no "
                "float64) on MLX's default device; use backend='torch' "
                "for float64 Fortran-parity runs."
            )
        if unknown:
            raise TypeError(
                f"AMICA.fit got unexpected keyword argument(s) {unknown}: "
                "neither a fit() parameter nor a constructor keyword of "
                f"{backend_cls.__name__}."
            )

        # Apply from_params_file's translated dict as per-call defaults
        # (issue #132 review item 2): an explicitly passed argument here
        # always wins, whether named (max_iter/lrate/do_mean/do_sphere/
        # do_newton, via the _UNSET sentinel) or in **kwargs (the selected
        # backend's constructor keywords, via plain dict membership). Settings
        # the file carries that apply to neither surface are named in one
        # warning rather than silently discarded.
        file_params = self._file_params or {}
        backend_defaults = _fit_defaults(backend_cls)

        def _file_default(explicit, name):
            if explicit is not _UNSET:
                return explicit
            return file_params.get(name, backend_defaults[name])

        max_iter = _file_default(max_iter, "max_iter")
        lrate = _file_default(lrate, "lrate")
        do_mean = _file_default(do_mean, "do_mean")
        do_sphere = _file_default(do_sphere, "do_sphere")
        do_newton = _file_default(do_newton, "do_newton")

        if file_params:
            for key, value in file_params.items():
                if key in ctor_params and key not in kwargs:
                    kwargs[key] = value
            handled = _FIT_NAMED_PARAMS | ctor_params | {"num_models", "num_mix"}
            unhandled = sorted(set(file_params) - handled)
            if unhandled:
                logger.warning(
                    "AMICA.fit: %d parameter-file setting(s) match neither a "
                    "fit()/%s parameter and were NOT applied "
                    "(informational only -- data-location metadata like "
                    "files/outdir/data_dim/num_comps is expected here; "
                    "anything else means this backend has no equivalent): %s",
                    len(unhandled),
                    backend_cls.__name__,
                    unhandled,
                )

        if self.verbose:
            print(f"Fitting AMICA with {n_channels} channels, {n_samples} samples")
            print(f"Models: {self.n_models}, Mixture components: {self.n_mix}")

        # The torch device as given (None picks one inside AMICATorchNG,
        # which also moves an automatic MPS pick for a float64 fit to the
        # CPU, issue #354); MLX always runs on MLX's default device.
        placement: dict[str, Any] = {}
        if self.backend == "torch":
            placement["device"] = self.device

        # Build and train the backend on a LOCAL reference first, and only
        # publish it to self (and derive the fitted-state attributes) once
        # fit() returns. If the backend constructor or fit() raises mid-training
        # (a numerical crash, OOM, singular sphere, interrupt, ...), self is left
        # untouched: a first fit keeps model_ is None (so the output methods
        # raise a clean "not fitted"), and a refit keeps the previous, known-good
        # model rather than a half-trained one falsely marked usable (issue #50
        # silent-failure review).
        backend = backend_cls(
            n_channels=n_channels,
            n_models=self.n_models,
            n_mix=self.n_mix,
            lrate=lrate,
            do_mean=do_mean,
            do_sphere=do_sphere,
            do_newton=do_newton,
            **placement,
            **kwargs,
        )
        backend.fit(X, max_iter=max_iter, verbose=self.verbose, mir_step=mir_step)

        self.model_ = backend
        self._mirror_backend()
        if not self.converged_:
            logger.warning(
                "AMICA.fit ended degenerate (stop_reason=%r) at iteration %d: the "
                "model holds non-finite parameters and cannot transform. Inspect "
                "stop_reason_/ll_history_; lower lrate, disable Newton, or check "
                "data conditioning, then refit.",
                self.stop_reason_,
                backend.iteration,
            )

        return self

    def _mirror_backend(self) -> None:
        """Copy the fitted backend's record onto the wrapper.

        Shared by :meth:`fit` and :meth:`load`, so both backends and both
        routes expose the same surface. A degenerate stop (one of the backend
        class's own ``_DEGENERATE_STOP_REASONS``) holds non-finite parameters
        and would return NaN sources, so it is not a usable model:
        ``is_fitted_`` stays False and the output methods refuse it
        (issue #50), while ``stop_reason_``/``converged_`` stay set for inspection.
        """
        backend = self.model_
        assert backend is not None
        self.ll_history_ = backend.ll_history
        self.final_ll_ = backend.final_ll_
        self.stop_reason_ = backend.stop_reason
        # mir_history_ is not persisted in state_dict() (a diagnostic
        # trajectory, not a fitted parameter), so a loaded model's is always
        # empty; exposed anyway for attribute-surface consistency.
        self.mir_history_ = backend.mir_history_
        # Restart records are persisted by state_dict (issue #198), so a loaded
        # best-of-N model can still say how its parameters were chosen.
        self.restart_seeds_ = backend.restart_seeds_
        self.restart_lls_ = backend.restart_lls_
        self.restart_stop_reasons_ = backend.restart_stop_reasons_
        self.converged_ = self.stop_reason_ not in backend._DEGENERATE_STOP_REASONS
        self.is_fitted_ = self.converged_

    def _check_usable(self, action: str) -> None:
        """Raise if the model cannot produce valid output: either never fitted,
        or the fit ended degenerate (a stop in the backend class's own
        ``_DEGENERATE_STOP_REASONS``, such as ``nan_ll``/``singular_ll``),
        leaving non-finite parameters that would yield NaN sources. This
        mirrors the backends' ``state_dict`` refusal to serialize a degenerate
        model (issue #50): a diverged fit fails loudly here instead of silently
        returning garbage."""
        if self.model_ is None:
            raise ValueError(f"Model must be fitted before {action}.")
        if self.stop_reason_ in self.model_._DEGENERATE_STOP_REASONS:
            raise RuntimeError(
                f"Refusing to {action}: fit ended degenerate "
                f"(stop_reason={self.stop_reason_!r}), so the model holds "
                f"non-finite parameters and would produce NaN output. Lower "
                f"lrate, disable Newton, or check data conditioning, then refit."
            )

    def __repr__(self) -> str:
        config = (
            f"backend={self.backend!r}, n_models={self.n_models}, n_mix={self.n_mix}"
        )
        if self.model_ is None:
            return f"<AMICA (unfitted, {config})>"
        if not self.converged_:
            return (
                f"<AMICA (degenerate fit, stop_reason={self.stop_reason_!r}, {config})>"
            )
        return f"<AMICA (fitted: {self.model_.n_channels} sources, {config})>"

    def transform(self, X: np.ndarray, model_idx: int = 0) -> np.ndarray:
        """
        Transform data to source space.

        Parameters
        ----------
        X : np.ndarray
            Input data of shape (n_channels, n_samples)
        model_idx : int, default=0
            Which model to use for transformation

        Returns
        -------
        S : np.ndarray
            Sources of shape (n_sources, n_samples): float64 from a default
            PyTorch fit, float32 from an MLX fit (its only precision).
        """
        self._check_usable("transform")
        assert self.model_ is not None

        return self.model_.transform(X, model_idx=model_idx)

    def fit_transform(self, X: np.ndarray, **fit_params) -> np.ndarray:
        """
        Fit model and transform data.

        Parameters
        ----------
        X : np.ndarray
            Input data of shape (n_channels, n_samples)
        **fit_params
            Parameters passed to fit()

        Returns
        -------
        S : np.ndarray
            Sources of shape (n_sources, n_samples)
        """
        self.fit(X, **fit_params)
        return self.transform(X)

    def get_mixing_matrix(self, model_idx: int = 0) -> np.ndarray:
        """
        Get the mixing matrix A in the sphered space.

        Column ``i`` is source ``i``'s mixing vector after sphering, the
        reference's ``A(:, comp_list(:, h))``. For scalp maps in input-channel
        space use :meth:`get_sensor_mixing_matrix`.

        Parameters
        ----------
        model_idx : int, default=0
            Which model's mixing matrix to return

        Returns
        -------
        A : np.ndarray
            Mixing matrix of shape (n_sources, n_sources)
        """
        self._check_usable("get the mixing matrix")
        assert self.model_ is not None

        return self.model_.get_mixing_matrix(model_idx=model_idx)

    def get_unmixing_matrix(self, model_idx: int = 0) -> np.ndarray:
        """
        Get the unmixing matrix W, which acts on sphered data.

        :meth:`transform` applies it after centering and sphering:
        ``S = W @ (get_sphere() @ (X - mean[:, None]) - c[:, None])``, with
        ``mean = get_mean()`` and ``c = get_model_center(model_idx)``, so
        ``W @ get_sphere()`` is only the linear part of that map.

        Parameters
        ----------
        model_idx : int, default=0
            Which model's unmixing matrix to return

        Returns
        -------
        W : np.ndarray
            Unmixing matrix of shape (n_sources, n_sources)
        """
        self._check_usable("get the unmixing matrix")
        assert self.model_ is not None

        return self.model_.get_unmixing_matrix(model_idx=model_idx)

    def get_sensor_mixing_matrix(self, model_idx: int = 0) -> np.ndarray:
        """
        Get the mixing matrix in input-channel space, ``pinv(sphere) @ A``.

        These are the scalp maps. Unlike :meth:`get_mixing_matrix` (sphered
        space), they stay valid after rank reduction, where the sphere is
        non-square (issue #223).

        Parameters
        ----------
        model_idx : int, default=0
            Which model's maps to return

        Returns
        -------
        A_sensor : np.ndarray
            Mixing matrix of shape (n_channels_in, n_sources)

        Raises
        ------
        ValueError
            If the model is unfitted, or ``model_idx`` is out of range.
        TypeError
            If ``model_idx`` is not an integer.
        RuntimeError
            If the fit ended degenerate (issue #50).
        """
        self._check_usable("get the sensor mixing matrix")
        assert self.model_ is not None

        return self.model_.get_sensor_mixing_matrix(model_idx=model_idx)

    def get_sphere(self) -> np.ndarray:
        """
        Get the fitted sphering matrix (issue #313).

        With :meth:`get_mean`, :meth:`get_model_center` and
        :meth:`get_unmixing_matrix` it composes :meth:`transform`:
        ``S = W @ (sphere @ (X - mean) - c)``.

        Returns
        -------
        sphere : np.ndarray of float64
            Shape (n_sources, n_channels_in), square unless the fit was
            rank-reduced. An MLX fit returns the float64 sphere its float32
            GPU copy was cast from.

        Raises
        ------
        ValueError
            If the model is unfitted.
        RuntimeError
            If the fit ended degenerate (issue #50).
        """
        self._check_usable("get the sphere")
        assert self.model_ is not None

        return self.model_.get_sphere()

    def get_mean(self) -> np.ndarray:
        """
        Get the per-channel mean removed before sphering (issue #313).

        Returns
        -------
        mean : np.ndarray of float64
            Shape (n_channels_in,); zeros for a ``do_mean=False`` fit, and
            float32 values for an MLX fit.

        Raises
        ------
        ValueError
            If the model is unfitted.
        RuntimeError
            If the fit ended degenerate (issue #50).
        """
        self._check_usable("get the mean")
        assert self.model_ is not None

        return self.model_.get_mean()

    def get_model_center(self, model_idx: int = 0) -> np.ndarray:
        """
        Get model ``model_idx``'s center ``c`` in the sphered space (issue #313).

        Parameters
        ----------
        model_idx : int, default=0
            Which model's center to return

        Returns
        -------
        c : np.ndarray of float64
            Shape (n_sources,); zeros for a single-model fit, and float32
            values for an MLX fit.

        Raises
        ------
        ValueError
            If the model is unfitted, or ``model_idx`` is out of range.
        TypeError
            If ``model_idx`` is not an integer.
        RuntimeError
            If the fit ended degenerate (issue #50).
        """
        self._check_usable("get the model center")
        assert self.model_ is not None

        return self.model_.get_model_center(model_idx=model_idx)

    def mir(
        self, X: np.ndarray, model_idx: int = 0, nbins: Optional[int] = None
    ) -> tuple:
        """
        Mutual Information Reduction (issue #137) of the fitted unmixing on ``X``.

        Composes the linear part of the raw-data-to-sources transform
        (unmixing @ sphere) and delegates to :func:`pamica.metrics.mir`. MIR
        is shift-invariant, so the mean and center that :meth:`transform`
        subtracts do not change it.

        Parameters
        ----------
        X : np.ndarray
            Raw (unpreprocessed) data of shape (n_channels, n_samples)
        model_idx : int, default=0
            Which model's unmixing to use
        nbins : int, optional
            Histogram bin count; see :func:`pamica.metrics.mir`

        Returns
        -------
        mir_nats : float
            Mutual information removed, in nats.
        variance : float
            Variance of the estimate.

        Raises
        ------
        ValueError
            If the model is unfitted; or if the fitted sphere is rank-reduced
            (explicit ``pcakeep``/``pcadb`` or automatic ``mineig``/
            ``mineig_rel`` detection), which leaves it rank-deficient so
            MIR's log-Jacobian term is undefined; or if ``X`` is non-finite or
            has a constant channel. See the backend's ``mir`` (:meth:`AMICATorchNG.mir`) and
            :func:`pamica.metrics.mir`.
        RuntimeError
            If the fit ended degenerate (issue #50), since the parameters are
            non-finite and any metric from them would be meaningless.
        """
        self._check_usable("compute MIR")
        assert self.model_ is not None

        return self.model_.mir(X, model_idx=model_idx, nbins=nbins)

    def pmi(
        self, X: np.ndarray, model_idx: int = 0, nbins: Optional[int] = None
    ) -> np.ndarray:
        """
        Pairwise Mutual Information (issue #137) between the fitted sources on ``X``.

        Delegates to :func:`pamica.metrics.pairwise_mi` on
        ``transform(X, model_idx)``.

        Parameters
        ----------
        X : np.ndarray
            Raw (unpreprocessed) data of shape (n_channels, n_samples)
        model_idx : int, default=0
            Which model's sources to use
        nbins : int, optional
            Histogram bin count; see :func:`pamica.metrics.pairwise_mi`

        Returns
        -------
        mi_matrix : np.ndarray of shape (n_sources, n_sources)
            Symmetric pairwise mutual information, in nats. The diagonal is each
            source's own entropy, not a mutual information; see
            :func:`pamica.viz.plot_pmi_heatmap`, which masks it by default.

        Raises
        ------
        ValueError
            If the model is unfitted; or if ``X`` is non-finite or has a
            constant channel; or if ``nbins`` is too large for the sample count
            (a sparse joint histogram fabricates mutual information from noise).
            See :func:`pamica.metrics.pairwise_mi`.
        RuntimeError
            If the fit ended degenerate (issue #50).
        """
        self._check_usable("compute PMI")
        assert self.model_ is not None

        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 the backend's ``model_loglik`` (:meth:`AMICATorchNG.model_loglik`,
        the same on either backend). 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 the backend's ``model_probability`` (see
        :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 get_pdftype(self, model_idx: int = 0) -> np.ndarray:
        """Per-source density-family code for model ``model_idx`` (issue #142).

        Delegates to the backend's ``get_pdftype`` (see
        :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 the backend's ``get_rho`` (see :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 the backend's ``shared_components`` (see
        :meth:`AMICATorchNG.shared_components`): one group of
        ``(model_idx, source_idx)`` pairs per shared component; 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]:
        """
        Component order by EEGLAB back-projected variance (IC1 = highest).

        Reports the display order EEGLAB's ``loadmodout15.m`` applies on load,
        without mutating the fitted parameters. Apply it to the columns of
        :meth:`get_mixing_matrix` (or rows of :meth:`get_unmixing_matrix`) to get
        EEGLAB-ordered components in Python.

        Parameters
        ----------
        model_idx : int, default=0
            Which model's components to order.
        return_svar : bool, default=False
            If True, also return the per-source variance sorted to ``order``.

        Returns
        -------
        order : np.ndarray of int
            Source indices, highest back-projected variance first.
        """
        self._check_usable("compute the variance order")
        assert self.model_ is not None

        return self.model_.variance_order(model_idx=model_idx, return_svar=return_svar)

    def write_amica_output(self, outdir: str) -> None:
        """
        Write the fitted model as an EEGLAB-readable AMICA output directory.

        Emits the raw binary files that EEGLAB's ``loadmodout15.m`` reads (``W``,
        ``S``, ``gm``, ``mean``, ``c``, ``alpha``, ``mu``, ``sbeta``, ``rho``,
        ``comp_list``, ``LL``), plus the sphered-space mixing matrix ``A``, so
        a pamica fit drops directly into an EEGLAB workflow
        (``mod = loadmodout15(outdir)``). ``loadmodout15`` applies the
        variance-ordering and normalization on load, so no manual re-ordering or
        sign-flipping is needed. Every file is in the Fortran reference's
        layout for any number of models, so single-model output is
        byte-compatible with the reference (issue #92).

        Also writes ``LLt`` (the per-sample/per-model log-likelihood,
        issue #155) for a model that was just fit in this process, taken from the
        E-step stash (issue #157); a model restored via :meth:`load` carries no
        stash, so ``LLt`` is omitted for it (a warning is logged). As in the
        reference, ``LLt`` is the E-step that produced ``final_ll_``: after a
        fit that ran to ``max_iter`` it is therefore one M-step older than the
        ``W``/``A`` written beside it, and after a convergence stop, which
        exits before that iteration's update, it belongs to them -- see
        ``docs/guides/amica-differences.md``. Use :meth:`model_loglik` for the
        log-likelihood of the written parameters.

        Parameters
        ----------
        outdir : str
            Destination directory (created if absent).
        """
        self._check_usable("write EEGLAB output")
        assert self.model_ is not None

        self.model_.write_amica_output(outdir)

    def save(self, filepath: str) -> None:
        """
        Save the fitted model to ``filepath`` via ``torch.save``.

        Persists the backend's ``state_dict`` (config, fitted arrays and fit
        record) plus the wrapper's own configuration, including which backend
        built the model, so :meth:`load` can fully reconstruct a
        transform-ready model of either backend. Everything written is a
        tensor or a plain Python primitive, so the file reloads with
        ``torch.load(weights_only=True)``: an MLX model's numpy arrays are
        stored as CPU tensors of the same dtype, and numpy scalars in the
        config or fit record (a ``seed=np.int64(...)``, say) as the equivalent
        Python numbers.

        The file is ``format_version`` 2 (issue #313), which records the
        backend; :meth:`load` still reads version 1 files, written before
        backend selection existed.

        Parameters
        ----------
        filepath : str
            Destination path (a ``.pt`` file by convention).

        Raises
        ------
        TypeError
            If the backend's state holds a value that ``weights_only`` loading
            could not read back (anything but a tensor, an array or a plain
            Python primitive), rather than writing a file :meth:`load` cannot
            open.
        """
        self._check_usable("save")
        assert self.model_ is not None

        backend = "torch" if isinstance(self.model_, AMICATorchNG) else "mlx"
        payload = {
            "format_version": _SAVE_FORMAT_VERSION,
            "wrapper": {
                "n_models": self.n_models,
                "n_mix": self.n_mix,
                "verbose": self.verbose,
                "backend": backend,
            },
            "backend": _state_to_payload(self.model_.state_dict()),
        }
        torch.save(payload, filepath)

    @classmethod
    def load(
        cls, filepath: str, device: Optional[Union[str, torch.device]] = None
    ) -> "AMICA":
        """
        Load a fitted model saved by :meth:`save`.

        The file records which backend built the model (``format_version``
        2), and the model comes back on that backend. A ``format_version`` 1
        file, written before backend selection existed (issue #313), holds a
        PyTorch model by construction and still loads. A backend payload
        saved before issue #334 (components as columns of ``A``) is converted
        on load, or refused with a request to refit when ``share_comps`` had
        merged components (see ``AMICATorchNG.from_state_dict``).

        Parameters
        ----------
        filepath : str
            Path to a file written by :meth:`save`.
        device : str or torch.device, optional
            Device to place a restored PyTorch model on. With ``None`` (auto),
            the backend constructor chooses, as in :meth:`fit`, so a float64
            parity model never lands on MPS. An MLX model always loads onto
            MLX's default device, so it must stay ``None`` for one.

        Returns
        -------
        amica : AMICA
            A fitted model ready for :meth:`transform` / :meth:`get_mixing_matrix`.

        Raises
        ------
        ValueError
            If the file's ``format_version`` is not 1 or 2, a section is
            missing, ``device`` is set for an MLX model, or the model was
            saved before issue #334 with merged components.
        ImportError
            If the file holds an MLX model and MLX is not installed.
        """
        payload = torch.load(filepath, weights_only=True)
        version = payload.get("format_version")
        if version not in _LOADABLE_FORMAT_VERSIONS:
            raise ValueError(
                f"unsupported AMICA save format_version: {version!r} (expected "
                f"one of {list(_LOADABLE_FORMAT_VERSIONS)})"
            )
        for key in ("wrapper", "backend"):
            if key not in payload:
                raise ValueError(
                    f"malformed AMICA save file {filepath!r}: missing {key!r} "
                    f"(format_version={version}); the file may be truncated or "
                    f"corrupted."
                )

        wrapper = payload["wrapper"]
        # Version 1 predates backend selection: the wrapper could only build
        # AMICATorchNG then, so every version 1 file is a torch payload.
        backend = "torch" if version == 1 else wrapper.get("backend")
        if backend not in _BACKENDS:
            raise ValueError(
                f"malformed AMICA save file {filepath!r}: wrapper['backend'] is "
                f"{backend!r}, expected one of {list(_BACKENDS)} "
                f"(format_version={version})."
            )
        if backend == "mlx":
            if device is not None:
                raise ValueError(
                    f"{filepath!r} holds an MLX-backend model, which always "
                    f"loads onto MLX's default device; device={device!r} "
                    "applies only to PyTorch-backend models, so pass "
                    "device=None."
                )
            if importlib.util.find_spec("mlx") is None:
                raise ImportError(
                    f"{filepath!r} holds an MLX-backend model (backend='mlx'), "
                    f"but MLX is not installed; {_MLX_INSTALL_HINT}."
                )

        model = cls(
            n_models=wrapper["n_models"],
            n_mix=wrapper["n_mix"],
            device=device,
            verbose=wrapper["verbose"],
            backend=backend,
        )
        state = payload["backend"]
        if backend == "torch":
            # device=None lets the constructor choose from the saved dtype, so
            # a float64 model's automatic MPS pick moves to the CPU as in fit().
            model.model_ = AMICATorchNG.from_state_dict(state, device=device)
        else:
            model.model_ = _mlx_backend_class().from_state_dict(
                _mlx_state_from_payload(state, filepath)
            )
        # state_dict() refuses to serialize a degenerate model, so a loaded
        # model is always usable; its stop_reason is carried through for
        # inspection anyway.
        model._mirror_backend()
        return model

    @classmethod
    def from_params_file(cls, params_file: str, **kwargs) -> "AMICA":
        """
        Create AMICA instance from a parameter file.

        Accepts two formats, auto-detected from ``params_file``'s *content*
        (issue #132): pamica's own JSON schema (``sample_data/
        sample_params.json``) and the literal Fortran ``input.param`` text
        format (``sample_data/input.param``), so the same file that drives
        the reference binary can drive pamica too. Detection always sniffs
        the content (JSON if it starts with ``{``/``[``, Fortran text
        otherwise) rather than trusting the extension -- a compact JSON file
        saved with a ``.param`` extension must still parse as JSON, not be
        silently misread as garbled Fortran text (issue #132 review item 3).
        Both formats are read through :func:`pamica.fortran_params.
        read_params_file` (issue #304), which also applies pamica's JSON
        schema's own alias spellings (``min_grad_norm``/``max_decs``/
        ``share_int``/...) to the canonical/constructor names, so a
        ``sample_params.json`` fit applies its ``max_decs``/
        ``min_grad_norm``/``share_int`` settings as ``maxdecs``/``min_nd``/
        ``share_iter``. See that function and
        :func:`pamica.fortran_params.read_fortran_param_file` for the
        Fortran-side key-mapping table and the deliberately-unmapped keys
        they warn about rather than silently drop.

        The full translated dict (beyond the ``n_models``/``n_mix`` used to
        size the instance here) is stashed on the returned instance and
        applied by :meth:`fit` as per-call defaults -- see ``fit``'s
        docstring for the precedence rule. Which settings apply is decided
        by the selected backend's own constructor signature, so
        ``from_params_file(path, backend="mlx")`` drives
        :class:`pamica.mlx_impl.AMICAMLXNG` from the same file (issue #313),
        and :meth:`fit` names any setting that backend cannot take in its
        "not applied" warning.

        Parameters
        ----------
        params_file : str
            Path to a JSON or Fortran-format parameter file.
        **kwargs
            Constructor arguments (``n_models``, ``n_mix``, ``device``,
            ``verbose``, ``backend``); ``n_models``/``n_mix`` given here
            override the file's.

        Returns
        -------
        amica : AMICA
            Configured AMICA instance
        """
        from .fortran_params import read_params_file

        params = read_params_file(Path(params_file))

        # Extract relevant parameters
        n_models = params.get("num_models", 1)
        n_mix = params.get("num_mix", 3)

        # Override with kwargs
        n_models = kwargs.pop("n_models", n_models)
        n_mix = kwargs.pop("n_mix", n_mix)

        model = cls(n_models=n_models, n_mix=n_mix, **kwargs)
        model._file_params = params
        return model

fit(X, max_iter=_UNSET, lrate=_UNSET, do_mean=_UNSET, do_sphere=_UNSET, do_newton=_UNSET, mir_step=0, **kwargs)

Fit AMICA model to data.

The defaults of max_iter, lrate, do_mean, do_sphere and do_newton are the selected backend's own, read from its signatures (issue #354), so AMICA().fit(X) fits exactly as the backend class does with its defaults. Both backends default to the values listed below.

Parameters:

Name	Type	Description	Default
X	ndarray	Input data of shape (n_channels, n_samples)	required
max_iter	int	Maximum number of iterations	100
lrate	float	Initial and maximum natural-gradient learning rate, the compiled amica15 default. EEGLAB's runamica15.m uses 0.05; see the defaults table in docs/guides/amica-differences.md.	0.1
do_mean	bool	Whether to remove mean from data	True
do_sphere	bool	Whether to sphere (whiten) the data	True
do_newton	bool	Whether to enable the Fortran-parity Newton preconditioner (tune via newt_start/newtrate in **kwargs).	False
mir_step	int	If > 0, compute MIR every mir_step iterations during training and record it in mir_history_ (issue #137). 0 (default) disables the waypoints; see :meth:AMICATorchNG.fit for details and the interaction with keep_best.	0
**kwargs		Additional parameters passed to the backend constructor, :class:AMICATorchNG or :class:pamica.mlx_impl.AMICAMLXNG (the two take the same keywords, except that dtype is PyTorch-only and raises ValueError with backend="mlx", which is float32-only; e.g. block_size, rho0, seed, pcakeep, use_min_dll, min_dll, maxincs, use_grad_norm, min_nd -- the issue #207 convergence stops, Fortran-faithful defaults True/1e-9/5/True/1e-7 -- or do_opt_block/blk_min/blk_max/blk_step, the issue #232 block-size search, off by default) -- the backend's tunables are constructor arguments, not fit() kwargs. A keyword the selected backend's constructor does not take raises TypeError. n_restarts (default 1) runs the fit from that many seeds and keeps the highest-likelihood one (issue #198), recording every restart in restart_seeds_/restart_lls_/ restart_stop_reasons_; it needs a base seed (or explicit restart_seeds) and costs proportionally more time, since restarts run serially. n_restarts=1 is bit-identical to a fit that never heard of restarts. Rank-deficient input (Maxwell-filtered MEG, average-referenced or interpolated EEG) is handled by mineig/mineig_rel (issue #223): the model is sized to the detected numerical rank and :meth:get_sensor_mixing_matrix maps components back to input channels. mineig is an absolute eigenvalue floor and so unit-dependent; pass mineig_rel for data far from unit scale. When the instance was built via :meth:from_params_file (issue #132), any of the parameters above -- named or in **kwargs -- left unset here falls back to that file's translated value instead of the backend's default; an explicitly passed argument always wins over the file. Settings the file carries that match neither a named fit() parameter nor a constructor keyword of the selected backend (data-location metadata like files/ outdir/data_dim, or a setting with no equivalent in that backend) are not applied and are named in a single logger.warning.	{}

Returns:

Name	Type	Description
self	AMICA	Fitted model

Raises:

Type	Description
ValueError	If X is not 2D, or dtype/device is passed with backend="mlx".
TypeError	If a keyword is neither a fit() parameter nor a constructor keyword of the selected backend.

Source code in pamica/amica.py

def fit(
    self,
    X: np.ndarray,
    max_iter=_UNSET,
    lrate=_UNSET,
    do_mean=_UNSET,
    do_sphere=_UNSET,
    do_newton=_UNSET,
    mir_step: int = 0,
    **kwargs,
) -> "AMICA":
    """
    Fit AMICA model to data.

    The defaults of ``max_iter``, ``lrate``, ``do_mean``, ``do_sphere``
    and ``do_newton`` are the selected backend's own, read from its
    signatures (issue #354), so ``AMICA().fit(X)`` fits exactly as the
    backend class does with its defaults. Both backends default to the
    values listed below.

    Parameters
    ----------
    X : np.ndarray
        Input data of shape (n_channels, n_samples)
    max_iter : int, default=100
        Maximum number of iterations
    lrate : float, default=0.1
        Initial and maximum natural-gradient learning rate, the compiled
        amica15 default. EEGLAB's ``runamica15.m`` uses 0.05; see the
        defaults table in ``docs/guides/amica-differences.md``.
    do_mean : bool, default=True
        Whether to remove mean from data
    do_sphere : bool, default=True
        Whether to sphere (whiten) the data
    do_newton : bool, default=False
        Whether to enable the Fortran-parity Newton preconditioner (tune
        via ``newt_start``/``newtrate`` in ``**kwargs``).
    mir_step : int, default=0
        If > 0, compute MIR every ``mir_step`` iterations during training
        and record it in ``mir_history_`` (issue #137). ``0`` (default)
        disables the waypoints; see :meth:`AMICATorchNG.fit` for details
        and the interaction with ``keep_best``.
    **kwargs
        Additional parameters passed to the backend constructor,
        :class:`AMICATorchNG` or :class:`pamica.mlx_impl.AMICAMLXNG`
        (the two take the same keywords, except that ``dtype`` is
        PyTorch-only and raises ``ValueError`` with ``backend="mlx"``,
        which is float32-only; e.g. ``block_size``, ``rho0``, ``seed``,
        ``pcakeep``, ``use_min_dll``, ``min_dll``, ``maxincs``, ``use_grad_norm``,
        ``min_nd`` -- the issue #207 convergence stops, Fortran-faithful
        defaults ``True``/``1e-9``/``5``/``True``/``1e-7`` -- or
        ``do_opt_block``/``blk_min``/``blk_max``/``blk_step``, the
        issue #232 block-size search, off by default) -- the backend's tunables
        are constructor arguments, not fit() kwargs. A keyword the
        selected backend's constructor does not take raises
        ``TypeError``.

        ``n_restarts`` (default 1) runs the fit from that many seeds and
        keeps the highest-likelihood one (issue #198), recording every
        restart in ``restart_seeds_``/``restart_lls_``/
        ``restart_stop_reasons_``; it needs a base ``seed`` (or explicit
        ``restart_seeds``) and costs proportionally more time, since
        restarts run serially. ``n_restarts=1`` is bit-identical to a fit
        that never heard of restarts.

        Rank-deficient input (Maxwell-filtered MEG, average-referenced or
        interpolated EEG) is handled by ``mineig``/``mineig_rel``
        (issue #223): the model is sized to the detected numerical rank and
        :meth:`get_sensor_mixing_matrix` maps components back to input
        channels. ``mineig`` is an absolute eigenvalue floor and so
        unit-dependent; pass ``mineig_rel`` for data far from unit scale.

        When the instance was built via :meth:`from_params_file`
        (issue #132), any of the parameters above -- named or in ``**kwargs`` --
        left unset here falls back to that file's translated value instead
        of the backend's default; an explicitly passed argument always
        wins over the file. Settings the file carries that match neither a
        named ``fit()`` parameter nor a constructor keyword of the
        selected backend (data-location metadata like ``files``/
        ``outdir``/``data_dim``, or a setting with no equivalent in that
        backend) are not applied and are named in a single
        ``logger.warning``.

    Returns
    -------
    self : AMICA
        Fitted model

    Raises
    ------
    ValueError
        If ``X`` is not 2D, or ``dtype``/``device`` is passed with
        ``backend="mlx"``.
    TypeError
        If a keyword is neither a ``fit()`` parameter nor a constructor
        keyword of the selected backend.
    """
    # Validate input
    if X.ndim != 2:
        raise ValueError(f"X must be 2D array, got shape {X.shape}")

    n_channels, n_samples = X.shape

    # Re-checked here, not only in __init__, so a backend/device attribute
    # changed after construction cannot slip past the validation.
    _check_backend(self.backend, self.device)
    backend_cls = _backend_class(self.backend)
    ctor_params = _ctor_params(backend_cls)
    # Both categories computed before either is raised (PR #347 review
    # item 4, the same fix as AMICA_NumPy's _reject_unknown_kwargs): a
    # first version of this check raised on torch_only alone, so
    # backend='mlx' with fit(dtype=..., blocksize=...) named only
    # 'dtype' and silently dropped 'blocksize' from the message.
    # unknown excludes torch_only so a torch-only name (which is not in
    # ctor_params for backend='mlx' either) is not double-counted.
    torch_only = (
        sorted(set(kwargs) & set(_TORCH_ONLY_PARAMS))
        if self.backend == "mlx"
        else []
    )
    unknown = sorted(set(kwargs) - ctor_params - set(torch_only))
    if torch_only and unknown:
        raise TypeError(
            f"AMICA.fit got invalid keyword argument(s): {torch_only} "
            "apply only to backend='torch' (the MLX backend computes in "
            "float32 only, on its default device; use backend='torch' "
            f"for float64 Fortran-parity runs), and {unknown} is "
            "neither a fit() parameter nor a constructor keyword of "
            f"{backend_cls.__name__}."
        )
    if torch_only:
        raise ValueError(
            f"{torch_only} apply only to backend='torch'. The MLX "
            "backend computes in float32 only (Apple GPUs have no "
            "float64) on MLX's default device; use backend='torch' "
            "for float64 Fortran-parity runs."
        )
    if unknown:
        raise TypeError(
            f"AMICA.fit got unexpected keyword argument(s) {unknown}: "
            "neither a fit() parameter nor a constructor keyword of "
            f"{backend_cls.__name__}."
        )

    # Apply from_params_file's translated dict as per-call defaults
    # (issue #132 review item 2): an explicitly passed argument here
    # always wins, whether named (max_iter/lrate/do_mean/do_sphere/
    # do_newton, via the _UNSET sentinel) or in **kwargs (the selected
    # backend's constructor keywords, via plain dict membership). Settings
    # the file carries that apply to neither surface are named in one
    # warning rather than silently discarded.
    file_params = self._file_params or {}
    backend_defaults = _fit_defaults(backend_cls)

    def _file_default(explicit, name):
        if explicit is not _UNSET:
            return explicit
        return file_params.get(name, backend_defaults[name])

    max_iter = _file_default(max_iter, "max_iter")
    lrate = _file_default(lrate, "lrate")
    do_mean = _file_default(do_mean, "do_mean")
    do_sphere = _file_default(do_sphere, "do_sphere")
    do_newton = _file_default(do_newton, "do_newton")

    if file_params:
        for key, value in file_params.items():
            if key in ctor_params and key not in kwargs:
                kwargs[key] = value
        handled = _FIT_NAMED_PARAMS | ctor_params | {"num_models", "num_mix"}
        unhandled = sorted(set(file_params) - handled)
        if unhandled:
            logger.warning(
                "AMICA.fit: %d parameter-file setting(s) match neither a "
                "fit()/%s parameter and were NOT applied "
                "(informational only -- data-location metadata like "
                "files/outdir/data_dim/num_comps is expected here; "
                "anything else means this backend has no equivalent): %s",
                len(unhandled),
                backend_cls.__name__,
                unhandled,
            )

    if self.verbose:
        print(f"Fitting AMICA with {n_channels} channels, {n_samples} samples")
        print(f"Models: {self.n_models}, Mixture components: {self.n_mix}")

    # The torch device as given (None picks one inside AMICATorchNG,
    # which also moves an automatic MPS pick for a float64 fit to the
    # CPU, issue #354); MLX always runs on MLX's default device.
    placement: dict[str, Any] = {}
    if self.backend == "torch":
        placement["device"] = self.device

    # Build and train the backend on a LOCAL reference first, and only
    # publish it to self (and derive the fitted-state attributes) once
    # fit() returns. If the backend constructor or fit() raises mid-training
    # (a numerical crash, OOM, singular sphere, interrupt, ...), self is left
    # untouched: a first fit keeps model_ is None (so the output methods
    # raise a clean "not fitted"), and a refit keeps the previous, known-good
    # model rather than a half-trained one falsely marked usable (issue #50
    # silent-failure review).
    backend = backend_cls(
        n_channels=n_channels,
        n_models=self.n_models,
        n_mix=self.n_mix,
        lrate=lrate,
        do_mean=do_mean,
        do_sphere=do_sphere,
        do_newton=do_newton,
        **placement,
        **kwargs,
    )
    backend.fit(X, max_iter=max_iter, verbose=self.verbose, mir_step=mir_step)

    self.model_ = backend
    self._mirror_backend()
    if not self.converged_:
        logger.warning(
            "AMICA.fit ended degenerate (stop_reason=%r) at iteration %d: the "
            "model holds non-finite parameters and cannot transform. Inspect "
            "stop_reason_/ll_history_; lower lrate, disable Newton, or check "
            "data conditioning, then refit.",
            self.stop_reason_,
            backend.iteration,
        )

    return self

transform(X, model_idx=0)

Transform data to source space.

Parameters:

Name	Type	Description	Default
X	ndarray	Input data of shape (n_channels, n_samples)	required
model_idx	int	Which model to use for transformation	0

Returns:

Name	Type	Description
S	ndarray	Sources of shape (n_sources, n_samples): float64 from a default PyTorch fit, float32 from an MLX fit (its only precision).

Source code in pamica/amica.py

def transform(self, X: np.ndarray, model_idx: int = 0) -> np.ndarray:
    """
    Transform data to source space.

    Parameters
    ----------
    X : np.ndarray
        Input data of shape (n_channels, n_samples)
    model_idx : int, default=0
        Which model to use for transformation

    Returns
    -------
    S : np.ndarray
        Sources of shape (n_sources, n_samples): float64 from a default
        PyTorch fit, float32 from an MLX fit (its only precision).
    """
    self._check_usable("transform")
    assert self.model_ is not None

    return self.model_.transform(X, model_idx=model_idx)

fit_transform(X, **fit_params)

Fit model and transform data.

Parameters:

Name	Type	Description	Default
X	ndarray	Input data of shape (n_channels, n_samples)	required
**fit_params		Parameters passed to fit()	{}

Returns:

Name	Type	Description
S	ndarray	Sources of shape (n_sources, n_samples)

Source code in pamica/amica.py

def fit_transform(self, X: np.ndarray, **fit_params) -> np.ndarray:
    """
    Fit model and transform data.

    Parameters
    ----------
    X : np.ndarray
        Input data of shape (n_channels, n_samples)
    **fit_params
        Parameters passed to fit()

    Returns
    -------
    S : np.ndarray
        Sources of shape (n_sources, n_samples)
    """
    self.fit(X, **fit_params)
    return self.transform(X)

get_mixing_matrix(model_idx=0)

Get the mixing matrix A in the sphered space.

Column i is source i's mixing vector after sphering, the reference's A(:, comp_list(:, h)). For scalp maps in input-channel space use :meth:get_sensor_mixing_matrix.

Parameters:

Name	Type	Description	Default
model_idx	int	Which model's mixing matrix to return	0

Returns:

Name	Type	Description
A	ndarray	Mixing matrix of shape (n_sources, n_sources)

Source code in pamica/amica.py

def get_mixing_matrix(self, model_idx: int = 0) -> np.ndarray:
    """
    Get the mixing matrix A in the sphered space.

    Column ``i`` is source ``i``'s mixing vector after sphering, the
    reference's ``A(:, comp_list(:, h))``. For scalp maps in input-channel
    space use :meth:`get_sensor_mixing_matrix`.

    Parameters
    ----------
    model_idx : int, default=0
        Which model's mixing matrix to return

    Returns
    -------
    A : np.ndarray
        Mixing matrix of shape (n_sources, n_sources)
    """
    self._check_usable("get the mixing matrix")
    assert self.model_ is not None

    return self.model_.get_mixing_matrix(model_idx=model_idx)

get_unmixing_matrix(model_idx=0)

Get the unmixing matrix W, which acts on sphered data.

:meth:transform applies it after centering and sphering: S = W @ (get_sphere() @ (X - mean[:, None]) - c[:, None]), with mean = get_mean() and c = get_model_center(model_idx), so W @ get_sphere() is only the linear part of that map.

Parameters:

Name	Type	Description	Default
model_idx	int	Which model's unmixing matrix to return	0

Returns:

Name	Type	Description
W	ndarray	Unmixing matrix of shape (n_sources, n_sources)

Source code in pamica/amica.py

def get_unmixing_matrix(self, model_idx: int = 0) -> np.ndarray:
    """
    Get the unmixing matrix W, which acts on sphered data.

    :meth:`transform` applies it after centering and sphering:
    ``S = W @ (get_sphere() @ (X - mean[:, None]) - c[:, None])``, with
    ``mean = get_mean()`` and ``c = get_model_center(model_idx)``, so
    ``W @ get_sphere()`` is only the linear part of that map.

    Parameters
    ----------
    model_idx : int, default=0
        Which model's unmixing matrix to return

    Returns
    -------
    W : np.ndarray
        Unmixing matrix of shape (n_sources, n_sources)
    """
    self._check_usable("get the unmixing matrix")
    assert self.model_ is not None

    return self.model_.get_unmixing_matrix(model_idx=model_idx)

get_sensor_mixing_matrix(model_idx=0)

Get the mixing matrix in input-channel space, pinv(sphere) @ A.

These are the scalp maps. Unlike :meth:get_mixing_matrix (sphered space), they stay valid after rank reduction, where the sphere is non-square (issue #223).

Parameters:

Name	Type	Description	Default
model_idx	int	Which model's maps to return	0

Returns:

Name	Type	Description
A_sensor	ndarray	Mixing matrix of shape (n_channels_in, n_sources)

Raises:

Type	Description
ValueError	If the model is unfitted, or model_idx is out of range.
TypeError	If model_idx is not an integer.
RuntimeError	If the fit ended degenerate (issue #50).

Source code in pamica/amica.py

def get_sensor_mixing_matrix(self, model_idx: int = 0) -> np.ndarray:
    """
    Get the mixing matrix in input-channel space, ``pinv(sphere) @ A``.

    These are the scalp maps. Unlike :meth:`get_mixing_matrix` (sphered
    space), they stay valid after rank reduction, where the sphere is
    non-square (issue #223).

    Parameters
    ----------
    model_idx : int, default=0
        Which model's maps to return

    Returns
    -------
    A_sensor : np.ndarray
        Mixing matrix of shape (n_channels_in, n_sources)

    Raises
    ------
    ValueError
        If the model is unfitted, or ``model_idx`` is out of range.
    TypeError
        If ``model_idx`` is not an integer.
    RuntimeError
        If the fit ended degenerate (issue #50).
    """
    self._check_usable("get the sensor mixing matrix")
    assert self.model_ is not None

    return self.model_.get_sensor_mixing_matrix(model_idx=model_idx)

get_sphere()

Get the fitted sphering matrix (issue #313).

With :meth:get_mean, :meth:get_model_center and :meth:get_unmixing_matrix it composes :meth:transform: S = W @ (sphere @ (X - mean) - c).

Returns:

Name	Type	Description
sphere	np.ndarray of float64	Shape (n_sources, n_channels_in), square unless the fit was rank-reduced. An MLX fit returns the float64 sphere its float32 GPU copy was cast from.

Raises:

Type	Description
ValueError	If the model is unfitted.
RuntimeError	If the fit ended degenerate (issue #50).

Source code in pamica/amica.py

def get_sphere(self) -> np.ndarray:
    """
    Get the fitted sphering matrix (issue #313).

    With :meth:`get_mean`, :meth:`get_model_center` and
    :meth:`get_unmixing_matrix` it composes :meth:`transform`:
    ``S = W @ (sphere @ (X - mean) - c)``.

    Returns
    -------
    sphere : np.ndarray of float64
        Shape (n_sources, n_channels_in), square unless the fit was
        rank-reduced. An MLX fit returns the float64 sphere its float32
        GPU copy was cast from.

    Raises
    ------
    ValueError
        If the model is unfitted.
    RuntimeError
        If the fit ended degenerate (issue #50).
    """
    self._check_usable("get the sphere")
    assert self.model_ is not None

    return self.model_.get_sphere()

get_mean()

Get the per-channel mean removed before sphering (issue #313).

Returns:

Name	Type	Description
mean	np.ndarray of float64	Shape (n_channels_in,); zeros for a do_mean=False fit, and float32 values for an MLX fit.

Raises:

Type	Description
ValueError	If the model is unfitted.
RuntimeError	If the fit ended degenerate (issue #50).

Source code in pamica/amica.py

def get_mean(self) -> np.ndarray:
    """
    Get the per-channel mean removed before sphering (issue #313).

    Returns
    -------
    mean : np.ndarray of float64
        Shape (n_channels_in,); zeros for a ``do_mean=False`` fit, and
        float32 values for an MLX fit.

    Raises
    ------
    ValueError
        If the model is unfitted.
    RuntimeError
        If the fit ended degenerate (issue #50).
    """
    self._check_usable("get the mean")
    assert self.model_ is not None

    return self.model_.get_mean()

get_model_center(model_idx=0)

Get model model_idx's center c in the sphered space (issue #313).

Parameters:

Name	Type	Description	Default
model_idx	int	Which model's center to return	0

Returns:

Name	Type	Description
c	np.ndarray of float64	Shape (n_sources,); zeros for a single-model fit, and float32 values for an MLX fit.

Raises:

Type	Description
ValueError	If the model is unfitted, or model_idx is out of range.
TypeError	If model_idx is not an integer.
RuntimeError	If the fit ended degenerate (issue #50).

Source code in pamica/amica.py

def get_model_center(self, model_idx: int = 0) -> np.ndarray:
    """
    Get model ``model_idx``'s center ``c`` in the sphered space (issue #313).

    Parameters
    ----------
    model_idx : int, default=0
        Which model's center to return

    Returns
    -------
    c : np.ndarray of float64
        Shape (n_sources,); zeros for a single-model fit, and float32
        values for an MLX fit.

    Raises
    ------
    ValueError
        If the model is unfitted, or ``model_idx`` is out of range.
    TypeError
        If ``model_idx`` is not an integer.
    RuntimeError
        If the fit ended degenerate (issue #50).
    """
    self._check_usable("get the model center")
    assert self.model_ is not None

    return self.model_.get_model_center(model_idx=model_idx)

mir(X, model_idx=0, nbins=None)

Mutual Information Reduction (issue #137) of the fitted unmixing on X.

Composes the linear part of the raw-data-to-sources transform (unmixing @ sphere) and delegates to :func:pamica.metrics.mir. MIR is shift-invariant, so the mean and center that :meth:transform subtracts do not change it.

Parameters:

Name	Type	Description	Default
X	ndarray	Raw (unpreprocessed) data of shape (n_channels, n_samples)	required
model_idx	int	Which model's unmixing to use	0
nbins	int	Histogram bin count; see :func:pamica.metrics.mir	None

Returns:

Name	Type	Description
mir_nats	float	Mutual information removed, in nats.
variance	float	Variance of the estimate.

Raises:

Type	Description
ValueError	If the model is unfitted; or if the fitted sphere is rank-reduced (explicit pcakeep/pcadb or automatic mineig/ mineig_rel detection), which leaves it rank-deficient so MIR's log-Jacobian term is undefined; or if X is non-finite or has a constant channel. See the backend's mir (:meth:AMICATorchNG.mir) and :func:pamica.metrics.mir.
RuntimeError	If the fit ended degenerate (issue #50), since the parameters are non-finite and any metric from them would be meaningless.

Source code in pamica/amica.py

def mir(
    self, X: np.ndarray, model_idx: int = 0, nbins: Optional[int] = None
) -> tuple:
    """
    Mutual Information Reduction (issue #137) of the fitted unmixing on ``X``.

    Composes the linear part of the raw-data-to-sources transform
    (unmixing @ sphere) and delegates to :func:`pamica.metrics.mir`. MIR
    is shift-invariant, so the mean and center that :meth:`transform`
    subtracts do not change it.

    Parameters
    ----------
    X : np.ndarray
        Raw (unpreprocessed) data of shape (n_channels, n_samples)
    model_idx : int, default=0
        Which model's unmixing to use
    nbins : int, optional
        Histogram bin count; see :func:`pamica.metrics.mir`

    Returns
    -------
    mir_nats : float
        Mutual information removed, in nats.
    variance : float
        Variance of the estimate.

    Raises
    ------
    ValueError
        If the model is unfitted; or if the fitted sphere is rank-reduced
        (explicit ``pcakeep``/``pcadb`` or automatic ``mineig``/
        ``mineig_rel`` detection), which leaves it rank-deficient so
        MIR's log-Jacobian term is undefined; or if ``X`` is non-finite or
        has a constant channel. See the backend's ``mir`` (:meth:`AMICATorchNG.mir`) and
        :func:`pamica.metrics.mir`.
    RuntimeError
        If the fit ended degenerate (issue #50), since the parameters are
        non-finite and any metric from them would be meaningless.
    """
    self._check_usable("compute MIR")
    assert self.model_ is not None

    return self.model_.mir(X, model_idx=model_idx, nbins=nbins)

pmi(X, model_idx=0, nbins=None)

Pairwise Mutual Information (issue #137) between the fitted sources on X.

Delegates to :func:pamica.metrics.pairwise_mi on transform(X, model_idx).

Parameters:

Name	Type	Description	Default
X	ndarray	Raw (unpreprocessed) data of shape (n_channels, n_samples)	required
model_idx	int	Which model's sources to use	0
nbins	int	Histogram bin count; see :func:pamica.metrics.pairwise_mi	None

Returns:

Name	Type	Description
mi_matrix	np.ndarray of shape (n_sources, n_sources)	Symmetric pairwise mutual information, in nats. The diagonal is each source's own entropy, not a mutual information; see :func:pamica.viz.plot_pmi_heatmap, which masks it by default.

Raises:

Type	Description
ValueError	If the model is unfitted; or if X is non-finite or has a constant channel; or if nbins is too large for the sample count (a sparse joint histogram fabricates mutual information from noise). See :func:pamica.metrics.pairwise_mi.
RuntimeError	If the fit ended degenerate (issue #50).

Source code in pamica/amica.py

def pmi(
    self, X: np.ndarray, model_idx: int = 0, nbins: Optional[int] = None
) -> np.ndarray:
    """
    Pairwise Mutual Information (issue #137) between the fitted sources on ``X``.

    Delegates to :func:`pamica.metrics.pairwise_mi` on
    ``transform(X, model_idx)``.

    Parameters
    ----------
    X : np.ndarray
        Raw (unpreprocessed) data of shape (n_channels, n_samples)
    model_idx : int, default=0
        Which model's sources to use
    nbins : int, optional
        Histogram bin count; see :func:`pamica.metrics.pairwise_mi`

    Returns
    -------
    mi_matrix : np.ndarray of shape (n_sources, n_sources)
        Symmetric pairwise mutual information, in nats. The diagonal is each
        source's own entropy, not a mutual information; see
        :func:`pamica.viz.plot_pmi_heatmap`, which masks it by default.

    Raises
    ------
    ValueError
        If the model is unfitted; or if ``X`` is non-finite or has a
        constant channel; or if ``nbins`` is too large for the sample count
        (a sparse joint histogram fabricates mutual information from noise).
        See :func:`pamica.metrics.pairwise_mi`.
    RuntimeError
        If the fit ended degenerate (issue #50).
    """
    self._check_usable("compute PMI")
    assert self.model_ is not None

    return self.model_.pmi(X, model_idx=model_idx, nbins=nbins)

model_loglik(X)

Per-model, per-sample log-likelihood Lht on X (issue #141).

Delegates to the backend's model_loglik (:meth:AMICATorchNG.model_loglik, the same on either backend). 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:

Name	Type	Description	Default
X	np.ndarray of shape (n_channels, n_samples)	Raw (unpreprocessed) data.	required

Returns:

Name	Type	Description
Lht	np.ndarray of shape (n_models, n_samples)

Raises:

Type	Description
ValueError	If the model is unfitted, or if X is non-finite.
RuntimeError	If the fit ended degenerate (issue #50).

Source code in pamica/amica.py

def model_loglik(self, X: np.ndarray) -> np.ndarray:
    """Per-model, per-sample log-likelihood ``Lht`` on ``X`` (issue #141).

    Delegates to the backend's ``model_loglik`` (:meth:`AMICATorchNG.model_loglik`,
    the same on either backend). 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)

model_probability(X)

Per-sample posterior probability of each model (issue #141).

Delegates to the backend's model_probability (see :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:

Name	Type	Description	Default
X	np.ndarray of shape (n_channels, n_samples)	Raw (unpreprocessed) data.	required

Returns:

Name	Type	Description
prob	np.ndarray of shape (n_models, n_samples)

Raises:

Type	Description
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).

Source code in pamica/amica.py

def model_probability(self, X: np.ndarray) -> np.ndarray:
    """Per-sample posterior probability of each model (issue #141).

    Delegates to the backend's ``model_probability`` (see
    :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)

get_pdftype(model_idx=0)

Per-source density-family code for model model_idx (issue #142).

Delegates to the backend's get_pdftype (see :meth:AMICATorchNG.get_pdftype). One integer per source component (0-4; see :data:pamica.torch_impl.PDFTYPE_NAMES).

Returns:

Type	Description
np.ndarray of int, shape (n_sources,)

Source code in pamica/amica.py

def get_pdftype(self, model_idx: int = 0) -> np.ndarray:
    """Per-source density-family code for model ``model_idx`` (issue #142).

    Delegates to the backend's ``get_pdftype`` (see
    :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)

get_rho(model_idx=0)

Generalized-Gaussian shape rho for model model_idx (issue #142).

Delegates to the backend's get_rho (see :meth:AMICATorchNG.get_rho).

Returns:

Type	Description
np.ndarray of float, shape (n_mix, n_sources)

Source code in pamica/amica.py

def get_rho(self, model_idx: int = 0) -> np.ndarray:
    """Generalized-Gaussian shape ``rho`` for model ``model_idx`` (issue #142).

    Delegates to the backend's ``get_rho`` (see :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)

shared_components()

Components shared across models by share_comps (issue #142).

Delegates to the backend's shared_components (see :meth:AMICATorchNG.shared_components): one group of (model_idx, source_idx) pairs per shared component; empty when nothing is shared.

Source code in pamica/amica.py

def shared_components(self) -> list:
    """Components shared across models by ``share_comps`` (issue #142).

    Delegates to the backend's ``shared_components`` (see
    :meth:`AMICATorchNG.shared_components`): one group of
    ``(model_idx, source_idx)`` pairs per shared component; empty when
    nothing is shared.
    """
    self._check_usable("get the shared components")
    assert self.model_ is not None

    return self.model_.shared_components()

variance_order(model_idx=0, return_svar=False)

Component order by EEGLAB back-projected variance (IC1 = highest).

Reports the display order EEGLAB's loadmodout15.m applies on load, without mutating the fitted parameters. Apply it to the columns of :meth:get_mixing_matrix (or rows of :meth:get_unmixing_matrix) to get EEGLAB-ordered components in Python.

Parameters:

Name	Type	Description	Default
model_idx	int	Which model's components to order.	0
return_svar	bool	If True, also return the per-source variance sorted to order.	False

Returns:

Name	Type	Description
order	np.ndarray of int	Source indices, highest back-projected variance first.

Source code in pamica/amica.py

def variance_order(
    self, model_idx: int = 0, return_svar: bool = False
) -> Union[np.ndarray, tuple]:
    """
    Component order by EEGLAB back-projected variance (IC1 = highest).

    Reports the display order EEGLAB's ``loadmodout15.m`` applies on load,
    without mutating the fitted parameters. Apply it to the columns of
    :meth:`get_mixing_matrix` (or rows of :meth:`get_unmixing_matrix`) to get
    EEGLAB-ordered components in Python.

    Parameters
    ----------
    model_idx : int, default=0
        Which model's components to order.
    return_svar : bool, default=False
        If True, also return the per-source variance sorted to ``order``.

    Returns
    -------
    order : np.ndarray of int
        Source indices, highest back-projected variance first.
    """
    self._check_usable("compute the variance order")
    assert self.model_ is not None

    return self.model_.variance_order(model_idx=model_idx, return_svar=return_svar)

write_amica_output(outdir)

Write the fitted model as an EEGLAB-readable AMICA output directory.

Emits the raw binary files that EEGLAB's loadmodout15.m reads (W, S, gm, mean, c, alpha, mu, sbeta, rho, comp_list, LL), plus the sphered-space mixing matrix A, so a pamica fit drops directly into an EEGLAB workflow (mod = loadmodout15(outdir)). loadmodout15 applies the variance-ordering and normalization on load, so no manual re-ordering or sign-flipping is needed. Every file is in the Fortran reference's layout for any number of models, so single-model output is byte-compatible with the reference (issue #92).

Also writes LLt (the per-sample/per-model log-likelihood, issue #155) for a model that was just fit in this process, taken from the E-step stash (issue #157); a model restored via :meth:load carries no stash, so LLt is omitted for it (a warning is logged). As in the reference, LLt is the E-step that produced final_ll_: after a fit that ran to max_iter it is therefore one M-step older than the W/A written beside it, and after a convergence stop, which exits before that iteration's update, it belongs to them -- see docs/guides/amica-differences.md. Use :meth:model_loglik for the log-likelihood of the written parameters.

Parameters:

Name	Type	Description	Default
outdir	str	Destination directory (created if absent).	required

Source code in pamica/amica.py

def write_amica_output(self, outdir: str) -> None:
    """
    Write the fitted model as an EEGLAB-readable AMICA output directory.

    Emits the raw binary files that EEGLAB's ``loadmodout15.m`` reads (``W``,
    ``S``, ``gm``, ``mean``, ``c``, ``alpha``, ``mu``, ``sbeta``, ``rho``,
    ``comp_list``, ``LL``), plus the sphered-space mixing matrix ``A``, so
    a pamica fit drops directly into an EEGLAB workflow
    (``mod = loadmodout15(outdir)``). ``loadmodout15`` applies the
    variance-ordering and normalization on load, so no manual re-ordering or
    sign-flipping is needed. Every file is in the Fortran reference's
    layout for any number of models, so single-model output is
    byte-compatible with the reference (issue #92).

    Also writes ``LLt`` (the per-sample/per-model log-likelihood,
    issue #155) for a model that was just fit in this process, taken from the
    E-step stash (issue #157); a model restored via :meth:`load` carries no
    stash, so ``LLt`` is omitted for it (a warning is logged). As in the
    reference, ``LLt`` is the E-step that produced ``final_ll_``: after a
    fit that ran to ``max_iter`` it is therefore one M-step older than the
    ``W``/``A`` written beside it, and after a convergence stop, which
    exits before that iteration's update, it belongs to them -- see
    ``docs/guides/amica-differences.md``. Use :meth:`model_loglik` for the
    log-likelihood of the written parameters.

    Parameters
    ----------
    outdir : str
        Destination directory (created if absent).
    """
    self._check_usable("write EEGLAB output")
    assert self.model_ is not None

    self.model_.write_amica_output(outdir)

save(filepath)

Save the fitted model to filepath via torch.save.

Persists the backend's state_dict (config, fitted arrays and fit record) plus the wrapper's own configuration, including which backend built the model, so :meth:load can fully reconstruct a transform-ready model of either backend. Everything written is a tensor or a plain Python primitive, so the file reloads with torch.load(weights_only=True): an MLX model's numpy arrays are stored as CPU tensors of the same dtype, and numpy scalars in the config or fit record (a seed=np.int64(...), say) as the equivalent Python numbers.

The file is format_version 2 (issue #313), which records the backend; :meth:load still reads version 1 files, written before backend selection existed.

Parameters:

Name	Type	Description	Default
filepath	str	Destination path (a .pt file by convention).	required

Raises:

Type	Description
TypeError	If the backend's state holds a value that weights_only loading could not read back (anything but a tensor, an array or a plain Python primitive), rather than writing a file :meth:load cannot open.

Source code in pamica/amica.py

def save(self, filepath: str) -> None:
    """
    Save the fitted model to ``filepath`` via ``torch.save``.

    Persists the backend's ``state_dict`` (config, fitted arrays and fit
    record) plus the wrapper's own configuration, including which backend
    built the model, so :meth:`load` can fully reconstruct a
    transform-ready model of either backend. Everything written is a
    tensor or a plain Python primitive, so the file reloads with
    ``torch.load(weights_only=True)``: an MLX model's numpy arrays are
    stored as CPU tensors of the same dtype, and numpy scalars in the
    config or fit record (a ``seed=np.int64(...)``, say) as the equivalent
    Python numbers.

    The file is ``format_version`` 2 (issue #313), which records the
    backend; :meth:`load` still reads version 1 files, written before
    backend selection existed.

    Parameters
    ----------
    filepath : str
        Destination path (a ``.pt`` file by convention).

    Raises
    ------
    TypeError
        If the backend's state holds a value that ``weights_only`` loading
        could not read back (anything but a tensor, an array or a plain
        Python primitive), rather than writing a file :meth:`load` cannot
        open.
    """
    self._check_usable("save")
    assert self.model_ is not None

    backend = "torch" if isinstance(self.model_, AMICATorchNG) else "mlx"
    payload = {
        "format_version": _SAVE_FORMAT_VERSION,
        "wrapper": {
            "n_models": self.n_models,
            "n_mix": self.n_mix,
            "verbose": self.verbose,
            "backend": backend,
        },
        "backend": _state_to_payload(self.model_.state_dict()),
    }
    torch.save(payload, filepath)

load(filepath, device=None) classmethod

Load a fitted model saved by :meth:save.

The file records which backend built the model (format_version 2), and the model comes back on that backend. A format_version 1 file, written before backend selection existed (issue #313), holds a PyTorch model by construction and still loads. A backend payload saved before issue #334 (components as columns of A) is converted on load, or refused with a request to refit when share_comps had merged components (see AMICATorchNG.from_state_dict).

Parameters:

Name	Type	Description	Default
filepath	str	Path to a file written by :meth:save.	required
device	str or device	Device to place a restored PyTorch model on. With None (auto), the backend constructor chooses, as in :meth:fit, so a float64 parity model never lands on MPS. An MLX model always loads onto MLX's default device, so it must stay None for one.	None

Returns:

Name	Type	Description
amica	AMICA	A fitted model ready for :meth:transform / :meth:get_mixing_matrix.

Raises:

Type	Description
ValueError	If the file's format_version is not 1 or 2, a section is missing, device is set for an MLX model, or the model was saved before issue #334 with merged components.
ImportError	If the file holds an MLX model and MLX is not installed.

Source code in pamica/amica.py

@classmethod
def load(
    cls, filepath: str, device: Optional[Union[str, torch.device]] = None
) -> "AMICA":
    """
    Load a fitted model saved by :meth:`save`.

    The file records which backend built the model (``format_version``
    2), and the model comes back on that backend. A ``format_version`` 1
    file, written before backend selection existed (issue #313), holds a
    PyTorch model by construction and still loads. A backend payload
    saved before issue #334 (components as columns of ``A``) is converted
    on load, or refused with a request to refit when ``share_comps`` had
    merged components (see ``AMICATorchNG.from_state_dict``).

    Parameters
    ----------
    filepath : str
        Path to a file written by :meth:`save`.
    device : str or torch.device, optional
        Device to place a restored PyTorch model on. With ``None`` (auto),
        the backend constructor chooses, as in :meth:`fit`, so a float64
        parity model never lands on MPS. An MLX model always loads onto
        MLX's default device, so it must stay ``None`` for one.

    Returns
    -------
    amica : AMICA
        A fitted model ready for :meth:`transform` / :meth:`get_mixing_matrix`.

    Raises
    ------
    ValueError
        If the file's ``format_version`` is not 1 or 2, a section is
        missing, ``device`` is set for an MLX model, or the model was
        saved before issue #334 with merged components.
    ImportError
        If the file holds an MLX model and MLX is not installed.
    """
    payload = torch.load(filepath, weights_only=True)
    version = payload.get("format_version")
    if version not in _LOADABLE_FORMAT_VERSIONS:
        raise ValueError(
            f"unsupported AMICA save format_version: {version!r} (expected "
            f"one of {list(_LOADABLE_FORMAT_VERSIONS)})"
        )
    for key in ("wrapper", "backend"):
        if key not in payload:
            raise ValueError(
                f"malformed AMICA save file {filepath!r}: missing {key!r} "
                f"(format_version={version}); the file may be truncated or "
                f"corrupted."
            )

    wrapper = payload["wrapper"]
    # Version 1 predates backend selection: the wrapper could only build
    # AMICATorchNG then, so every version 1 file is a torch payload.
    backend = "torch" if version == 1 else wrapper.get("backend")
    if backend not in _BACKENDS:
        raise ValueError(
            f"malformed AMICA save file {filepath!r}: wrapper['backend'] is "
            f"{backend!r}, expected one of {list(_BACKENDS)} "
            f"(format_version={version})."
        )
    if backend == "mlx":
        if device is not None:
            raise ValueError(
                f"{filepath!r} holds an MLX-backend model, which always "
                f"loads onto MLX's default device; device={device!r} "
                "applies only to PyTorch-backend models, so pass "
                "device=None."
            )
        if importlib.util.find_spec("mlx") is None:
            raise ImportError(
                f"{filepath!r} holds an MLX-backend model (backend='mlx'), "
                f"but MLX is not installed; {_MLX_INSTALL_HINT}."
            )

    model = cls(
        n_models=wrapper["n_models"],
        n_mix=wrapper["n_mix"],
        device=device,
        verbose=wrapper["verbose"],
        backend=backend,
    )
    state = payload["backend"]
    if backend == "torch":
        # device=None lets the constructor choose from the saved dtype, so
        # a float64 model's automatic MPS pick moves to the CPU as in fit().
        model.model_ = AMICATorchNG.from_state_dict(state, device=device)
    else:
        model.model_ = _mlx_backend_class().from_state_dict(
            _mlx_state_from_payload(state, filepath)
        )
    # state_dict() refuses to serialize a degenerate model, so a loaded
    # model is always usable; its stop_reason is carried through for
    # inspection anyway.
    model._mirror_backend()
    return model

from_params_file(params_file, **kwargs) classmethod

Create AMICA instance from a parameter file.

Accepts two formats, auto-detected from params_file's content (issue #132): pamica's own JSON schema (sample_data/ sample_params.json) and the literal Fortran input.param text format (sample_data/input.param), so the same file that drives the reference binary can drive pamica too. Detection always sniffs the content (JSON if it starts with {/[, Fortran text otherwise) rather than trusting the extension -- a compact JSON file saved with a .param extension must still parse as JSON, not be silently misread as garbled Fortran text (issue #132 review item 3). Both formats are read through :func:pamica.fortran_params. read_params_file (issue #304), which also applies pamica's JSON schema's own alias spellings (min_grad_norm/max_decs/ share_int/...) to the canonical/constructor names, so a sample_params.json fit applies its max_decs/ min_grad_norm/share_int settings as maxdecs/min_nd/ share_iter. See that function and :func:pamica.fortran_params.read_fortran_param_file for the Fortran-side key-mapping table and the deliberately-unmapped keys they warn about rather than silently drop.

The full translated dict (beyond the n_models/n_mix used to size the instance here) is stashed on the returned instance and applied by :meth:fit as per-call defaults -- see fit's docstring for the precedence rule. Which settings apply is decided by the selected backend's own constructor signature, so from_params_file(path, backend="mlx") drives :class:pamica.mlx_impl.AMICAMLXNG from the same file (issue #313), and :meth:fit names any setting that backend cannot take in its "not applied" warning.

Parameters:

Name	Type	Description	Default
params_file	str	Path to a JSON or Fortran-format parameter file.	required
**kwargs		Constructor arguments (n_models, n_mix, device, verbose, backend); n_models/n_mix given here override the file's.	{}

Returns:

Name	Type	Description
amica	AMICA	Configured AMICA instance

Source code in pamica/amica.py

@classmethod
def from_params_file(cls, params_file: str, **kwargs) -> "AMICA":
    """
    Create AMICA instance from a parameter file.

    Accepts two formats, auto-detected from ``params_file``'s *content*
    (issue #132): pamica's own JSON schema (``sample_data/
    sample_params.json``) and the literal Fortran ``input.param`` text
    format (``sample_data/input.param``), so the same file that drives
    the reference binary can drive pamica too. Detection always sniffs
    the content (JSON if it starts with ``{``/``[``, Fortran text
    otherwise) rather than trusting the extension -- a compact JSON file
    saved with a ``.param`` extension must still parse as JSON, not be
    silently misread as garbled Fortran text (issue #132 review item 3).
    Both formats are read through :func:`pamica.fortran_params.
    read_params_file` (issue #304), which also applies pamica's JSON
    schema's own alias spellings (``min_grad_norm``/``max_decs``/
    ``share_int``/...) to the canonical/constructor names, so a
    ``sample_params.json`` fit applies its ``max_decs``/
    ``min_grad_norm``/``share_int`` settings as ``maxdecs``/``min_nd``/
    ``share_iter``. See that function and
    :func:`pamica.fortran_params.read_fortran_param_file` for the
    Fortran-side key-mapping table and the deliberately-unmapped keys
    they warn about rather than silently drop.

    The full translated dict (beyond the ``n_models``/``n_mix`` used to
    size the instance here) is stashed on the returned instance and
    applied by :meth:`fit` as per-call defaults -- see ``fit``'s
    docstring for the precedence rule. Which settings apply is decided
    by the selected backend's own constructor signature, so
    ``from_params_file(path, backend="mlx")`` drives
    :class:`pamica.mlx_impl.AMICAMLXNG` from the same file (issue #313),
    and :meth:`fit` names any setting that backend cannot take in its
    "not applied" warning.

    Parameters
    ----------
    params_file : str
        Path to a JSON or Fortran-format parameter file.
    **kwargs
        Constructor arguments (``n_models``, ``n_mix``, ``device``,
        ``verbose``, ``backend``); ``n_models``/``n_mix`` given here
        override the file's.

    Returns
    -------
    amica : AMICA
        Configured AMICA instance
    """
    from .fortran_params import read_params_file

    params = read_params_file(Path(params_file))

    # Extract relevant parameters
    n_models = params.get("num_models", 1)
    n_mix = params.get("num_mix", 3)

    # Override with kwargs
    n_models = kwargs.pop("n_models", n_models)
    n_mix = kwargs.pop("n_mix", n_mix)

    model = cls(n_models=n_models, n_mix=n_mix, **kwargs)
    model._file_params = params
    return model