Native backend (AMICANative)

AMICANative runs the real AMICA Fortran reference binary as a fourth backend, alongside the PyTorch, MLX and NumPy backends. Because it executes the literal reference implementation, it is the strongest possible parity oracle: rather than checking pamica against Fortran, you can run Fortran as a pamica backend, with no separate toolchain to build or drive.

The reference is built dependency-free (a single-rank MPI shim removes the Open MPI runtime, on top of the no-MKL recipe) and released as a self-contained binary per platform. The engine resolves the binary for your host, downloads it from the GitHub release on first use (verifying its SHA-256), caches it, and runs it.

import numpy as np
from pamica import AMICANative

# X is (n_channels, n_samples) of real EEG/EMG.
model = AMICANative(n_models=1, n_mix=3).fit(X)

out = model.output_             # an AmicaOutput (W, A, mixing/unmixing, LL, ...)
sources = model.transform(X)    # source activations, EEGLAB variance order

fit writes the data and a full input.param (so it does not depend on an installed sample_data), runs the binary, and loads the result with loadmodout, exposed as model.output_ (a pamica.numpy_impl.load.AmicaOutput). n_models and n_mix are friendly aliases; any Fortran input.param field (max_iter, lrate, pdftype, do_newton, ...) can be passed as a keyword. A collapsed fit (non-finite weights) is raised as a clear degenerate-fit error rather than an opaque SVD failure.

The input.param it writes carries pamica's shared defaults (issue #354), the values the PyTorch, MLX and NumPy backends default to, under the binary's keywords: lrate 0.1, Newton off, max_iter 100, do_opt_block off, pdftype 0, and so on (Default settings has the full table). The binary's block_size counts one thread's share of a block, and it processes no block at all, writing an all-NaN fit, when max_threads * block_size exceeds the samples (issue #292). Unless you pass block_size, the engine writes pamica's 8192-sample block, capped at the data's length, divided by max_threads (10 unless given), so a default run keeps pamica's block and never meets that case; a block_size you pass is written as is. pamica settings the binary has no keyword for, such as keep_best and mineig_rel, are not written, so the binary returns its last iterate and applies the absolute mineig floor. Keys only the binary reads (max_threads, byte_size, writestep, the load_* and update_* switches) keep the bundled input.param's values. Before issue #354 the engine wrote the bundled pamica/sample_data/input.param's settings instead (lrate 0.05, Newton on from iteration 50, max_iter 2000, block_size 512). To run as that file or an EEGLAB-written one configures the binary, pass the file's settings as keywords (Reproducing an EEGLAB run shows how).

Binary resolution and caching

On the first run for a given host, the engine downloads the matching release asset and caches it, so later runs are offline:

• Cache location: ~/.cache/pamica/bin/<version>/ (or $XDG_CACHE_HOME/pamica/bin/; override with PAMICA_NATIVE_CACHE).
• Integrity: the download is staged in a temporary directory, verified against its .sha256 release asset, marked executable, and only then moved atomically into the cache. A file therefore exists at the cache path only once it has passed its checksum; a failed or tampered download never runs.
• Platforms: prebuilt binaries are attached to each release for macOS arm64, Linux x64, Linux arm64 and Windows x64. Windows arm64 has no native Fortran toolchain yet (#173); it maps to the x64 binary, which runs under Windows 11 ARM's x64 emulation.

Install the binary explicitly (for example to pre-populate the cache in an offline or CI environment):

python -m pamica.native                 # download + cache the latest release binary
python -m pamica.native --version v0.2.1 # a specific release
python -m pamica.native --print         # print where it resolves to; do not download

Using a local build

Set PAMICA_NATIVE_BINARY to a locally built binary to bypass the resolver entirely (this is what the tests use, and the fallback on a platform with no prebuilt asset):

export PAMICA_NATIVE_BINARY=/path/to/amica15

Build one with native/build.sh (gfortran + LAPACK; the single-rank MPI shim means no MPI runtime is required). If no prebuilt binary matches your host and PAMICA_NATIVE_BINARY is unset, the resolver raises a clear, actionable error pointing at the build script.

Validation harness

validate_implementations.py can source the reference through this engine instead of the bundled macOS-only amica15mac fixture, so the real Fortran reference can be compared against pamica's backends on any platform (--backend torch, the default, or numpy, mlx, a comma-separated list, or all):

# Resolve/download the native binary (or honor PAMICA_NATIVE_BINARY):
uv run python validate_implementations.py --native-engine

# Or point at a specific binary, and compare every backend against it:
uv run python validate_implementations.py --fortran-binary /path/to/amica15 --backend all

pamica.native.engine.AMICANative

Run the native AMICA binary on data and expose the result as an AmicaOutput.

Parameters:

Name	Type	Description	Default
binary	path - like	Explicit binary path; otherwise resolved for the host (downloaded from the release on first use). Equivalent to setting PAMICA_NATIVE_BINARY.	None
version	str	Release tag to resolve the binary from when not given explicitly.	"latest"
threads	int	OMP_NUM_THREADS for the run (default: the binary's own default).	None
timeout	float	Seconds before the subprocess is killed (default: no timeout).	None
**params	object	Any Fortran input.param field (or a friendly alias: n_models, n_mix), overriding the defaults; e.g. max_iter, lrate, pdftype, do_newton.	{}

Notes

The defaults are pamica's shared ones since issue #354, the values the PyTorch, MLX and NumPy backends default to, under the binary's keywords (lrate 0.1, Newton off, max_iter 100, do_opt_block off, pdftype 0, ...). pamica settings the binary has no keyword for, such as keep_best and mineig_rel, are not written, so the binary returns its last iterate and uses the absolute mineig floor.

The binary's block_size counts one thread's share of a block, and it processes no block at all when max_threads * block_size exceeds the samples (an all-NaN fit, issue #292). Unless block_size is given, the engine therefore writes pamica's 8192-sample block, capped at the data's length, divided by max_threads (10 unless given): 819 for a long recording. A given block_size is written as is. Before issue #354 the defaults were the bundled sample_data/input.param's (lrate 0.05, Newton on from iteration 50, max_iter 2000, block_size 512); pass that file's settings as keywords to run it (see docs/api/native-backend.md).

Source code in pamica/native/engine.py

class AMICANative:
    """Run the native AMICA binary on data and expose the result as an
    ``AmicaOutput``.

    Parameters
    ----------
    binary : path-like, optional
        Explicit binary path; otherwise resolved for the host (downloaded from the
        release on first use). Equivalent to setting ``PAMICA_NATIVE_BINARY``.
    version : str, default "latest"
        Release tag to resolve the binary from when not given explicitly.
    threads : int, optional
        ``OMP_NUM_THREADS`` for the run (default: the binary's own default).
    timeout : float, optional
        Seconds before the subprocess is killed (default: no timeout).
    **params
        Any Fortran ``input.param`` field (or a friendly alias: ``n_models``,
        ``n_mix``), overriding the defaults; e.g. ``max_iter``, ``lrate``,
        ``pdftype``, ``do_newton``.

    Notes
    -----
    The defaults are pamica's shared ones since issue #354, the values the
    PyTorch, MLX and NumPy backends default to, under the binary's keywords
    (``lrate`` 0.1, Newton off, ``max_iter`` 100, ``do_opt_block`` off,
    ``pdftype`` 0, ...). pamica settings the binary has no keyword for, such
    as ``keep_best`` and ``mineig_rel``, are not written, so the binary
    returns its last iterate and uses the absolute ``mineig`` floor.

    The binary's ``block_size`` counts one thread's share of a block, and it
    processes no block at all when ``max_threads * block_size`` exceeds the
    samples (an all-NaN fit, issue #292). Unless ``block_size`` is given, the
    engine therefore writes pamica's 8192-sample block, capped at the data's
    length, divided by ``max_threads`` (10 unless given): 819 for a long
    recording. A given ``block_size`` is written as is. Before issue #354 the defaults were the bundled
    ``sample_data/input.param``'s (``lrate`` 0.05, Newton on from iteration
    50, ``max_iter`` 2000, ``block_size`` 512); pass that file's settings as
    keywords to run it (see ``docs/api/native-backend.md``).
    """

    def __init__(
        self,
        binary: Optional[Union[str, Path]] = None,
        *,
        version: str = "latest",
        threads: Optional[int] = None,
        timeout: Optional[float] = None,
        **params: object,
    ) -> None:
        # resolve() now: the subprocess runs with cwd set to a tempdir, so a
        # relative binary path would pass the existence check but fail to launch.
        self.binary = Path(binary).resolve() if binary is not None else None
        self.version = version
        self.threads = threads
        self.timeout = timeout
        self.params = params
        self.output_: Optional[AmicaOutput] = None

    def _resolve_binary(self) -> Path:
        if self.binary is not None:
            if not self.binary.exists():
                raise FileNotFoundError(f"binary not found: {self.binary}")
            return self.binary
        return resolver.resolve(self.version)

    def _input_params(
        self, n_channels: int, n_samples: int, params: dict[str, object]
    ) -> dict[str, object]:
        """The ``input.param`` settings :meth:`fit` writes for data of this
        shape: pamica's shared defaults, then the constructor's and ``fit``'s
        keywords (friendly aliases mapped), then the data's dimensions."""
        merged = _default_params()
        given = set()
        for src in (self.params, params):
            for key, value in src.items():
                merged[_ALIASES.get(key, key)] = value
                given.add(_ALIASES.get(key, key))
        if "block_size" not in given:
            # pamica's block_size counts the samples of one block; the
            # binary's counts one thread's share of it, and it runs
            # n_samples // (max_threads * block_size) blocks, so a block larger
            # than the data would leave it none and an all-NaN fit (issue
            # #292). The default is therefore pamica's block, at most the
            # data, split over the threads.
            # int(str(...)): a value passed as a file's text ("10") reads the same.
            threads = int(str(merged["max_threads"]))
            block = min(int(str(merged["block_size"])), n_samples)
            merged["block_size"] = max(1, block // threads)
        merged["data_dim"] = n_channels
        merged["field_dim"] = n_samples
        if not merged.get("pcakeep"):
            merged["pcakeep"] = n_channels  # default: keep all components
        # `files` must come first: amica15.f90 hard-stops if it parses other
        # keys before the data file. Dict insertion order preserves that.
        return {"files": "./data.fdt", "outdir": "./amicaout/", **merged}

    def fit(self, X: np.ndarray, **params: object) -> "AMICANative":
        """Run AMICA on ``X`` (shape ``(n_channels, n_samples)``) and store the
        result as ``self.output_`` (an ``AmicaOutput``)."""
        X = np.asarray(X)
        if X.ndim != 2:
            raise ValueError(f"X must be 2-D (n_channels, n_samples); got {X.shape}")
        n_channels, n_samples = X.shape
        param = self._input_params(n_channels, n_samples, params)

        binary = self._resolve_binary()

        with tempfile.TemporaryDirectory(prefix="amica_native_") as td:
            work = Path(td)
            # AMICA reads the data as raw byte_size floats in column-major order
            # (numpy_impl/data.py: reshape order="F"); write it that way.
            byte_size = param["byte_size"]
            dtype = (
                np.float32
                if isinstance(byte_size, int) and byte_size == 4
                else np.float64
            )
            X.astype(dtype).ravel(order="F").tofile(work / "data.fdt")

            outdir = work / "amicaout"
            outdir.mkdir()
            (work / "input.param").write_text(_render_param(param))

            env = None
            if self.threads is not None:
                import os

                env = {**os.environ, "OMP_NUM_THREADS": str(self.threads)}

            proc = subprocess.run(
                [str(binary), "input.param"],
                cwd=work,
                capture_output=True,
                text=True,
                timeout=self.timeout,
                env=env,
            )
            if proc.returncode != 0:
                raise RuntimeError(
                    f"native AMICA failed (exit {proc.returncode}).\n"
                    f"stdout tail:\n{proc.stdout[-2000:]}\n"
                    f"stderr tail:\n{proc.stderr[-2000:]}"
                )
            if not (outdir / "W").exists():
                raise RuntimeError(
                    "native AMICA produced no output (no 'W' file); stdout tail:\n"
                    f"{proc.stdout[-2000:]}"
                )
            # A collapsed fit writes NaN weights (and zero model probabilities);
            # loadmodout's pinv(W@S) would then fail with an opaque SVD error, so
            # detect it here and report it as the degenerate fit it is (cf. the
            # #50 degenerate-fit contract for the Python backends). An empty W
            # (truncated write) is caught too -- an all-NaN check vacuously passes
            # on a zero-length array.
            w_raw = np.fromfile(outdir / "W")
            if w_raw.size == 0 or not np.all(np.isfinite(w_raw)):
                raise RuntimeError(
                    "native AMICA produced a degenerate fit (non-finite weights); "
                    "the run did not converge. Try more iterations, more data, or "
                    "fewer models."
                )
            self.output_ = loadmodout(outdir)
        return self

    def transform(self, X: np.ndarray, model_idx: int = 0) -> np.ndarray:
        """Source activations for ``X`` from the fitted model (delegates to
        ``AmicaOutput.sources``)."""
        if self.output_ is None:
            raise RuntimeError("call fit() before transform().")
        return self.output_.sources(X, model_idx)

fit(X, **params)

Run AMICA on X (shape (n_channels, n_samples)) and store the result as self.output_ (an AmicaOutput).

Source code in pamica/native/engine.py

def fit(self, X: np.ndarray, **params: object) -> "AMICANative":
    """Run AMICA on ``X`` (shape ``(n_channels, n_samples)``) and store the
    result as ``self.output_`` (an ``AmicaOutput``)."""
    X = np.asarray(X)
    if X.ndim != 2:
        raise ValueError(f"X must be 2-D (n_channels, n_samples); got {X.shape}")
    n_channels, n_samples = X.shape
    param = self._input_params(n_channels, n_samples, params)

    binary = self._resolve_binary()

    with tempfile.TemporaryDirectory(prefix="amica_native_") as td:
        work = Path(td)
        # AMICA reads the data as raw byte_size floats in column-major order
        # (numpy_impl/data.py: reshape order="F"); write it that way.
        byte_size = param["byte_size"]
        dtype = (
            np.float32
            if isinstance(byte_size, int) and byte_size == 4
            else np.float64
        )
        X.astype(dtype).ravel(order="F").tofile(work / "data.fdt")

        outdir = work / "amicaout"
        outdir.mkdir()
        (work / "input.param").write_text(_render_param(param))

        env = None
        if self.threads is not None:
            import os

            env = {**os.environ, "OMP_NUM_THREADS": str(self.threads)}

        proc = subprocess.run(
            [str(binary), "input.param"],
            cwd=work,
            capture_output=True,
            text=True,
            timeout=self.timeout,
            env=env,
        )
        if proc.returncode != 0:
            raise RuntimeError(
                f"native AMICA failed (exit {proc.returncode}).\n"
                f"stdout tail:\n{proc.stdout[-2000:]}\n"
                f"stderr tail:\n{proc.stderr[-2000:]}"
            )
        if not (outdir / "W").exists():
            raise RuntimeError(
                "native AMICA produced no output (no 'W' file); stdout tail:\n"
                f"{proc.stdout[-2000:]}"
            )
        # A collapsed fit writes NaN weights (and zero model probabilities);
        # loadmodout's pinv(W@S) would then fail with an opaque SVD error, so
        # detect it here and report it as the degenerate fit it is (cf. the
        # #50 degenerate-fit contract for the Python backends). An empty W
        # (truncated write) is caught too -- an all-NaN check vacuously passes
        # on a zero-length array.
        w_raw = np.fromfile(outdir / "W")
        if w_raw.size == 0 or not np.all(np.isfinite(w_raw)):
            raise RuntimeError(
                "native AMICA produced a degenerate fit (non-finite weights); "
                "the run did not converge. Try more iterations, more data, or "
                "fewer models."
            )
        self.output_ = loadmodout(outdir)
    return self

transform(X, model_idx=0)

Source activations for X from the fitted model (delegates to AmicaOutput.sources).

Source code in pamica/native/engine.py

def transform(self, X: np.ndarray, model_idx: int = 0) -> np.ndarray:
    """Source activations for ``X`` from the fitted model (delegates to
    ``AmicaOutput.sources``)."""
    if self.output_ is None:
        raise RuntimeError("call fit() before transform().")
    return self.output_.sources(X, model_idx)