|
| 1 | +"""``AMICANative``: run the native AMICA Fortran binary as a pyAMICA backend. |
| 2 | +
|
| 3 | +The fourth run engine, alongside the NumPy, PyTorch (``AMICATorchNG``) and MLX |
| 4 | +backends. It writes the data and an ``input.param``, runs the dependency-free |
| 5 | +binary resolved by :mod:`pyAMICA.native.resolver`, and reads the result back |
| 6 | +through :func:`pyAMICA.numpy_impl.load.loadmodout` -- so its output is an |
| 7 | +``AmicaOutput`` with the same accessors (``.sources``, ``.W``, ``.A``, ...) as a |
| 8 | +loaded fit. This is the Fortran reference itself, so it is the parity oracle the |
| 9 | +Python backends are validated against. |
| 10 | +""" |
| 11 | + |
| 12 | +from __future__ import annotations |
| 13 | + |
| 14 | +import subprocess |
| 15 | +import tempfile |
| 16 | +from pathlib import Path |
| 17 | +from typing import Optional, Union |
| 18 | + |
| 19 | +import numpy as np |
| 20 | + |
| 21 | +from ..numpy_impl.load import AmicaOutput, loadmodout |
| 22 | +from . import resolver |
| 23 | + |
| 24 | +# Full default parameter set (mirrors sample_data/input.param) so the engine does |
| 25 | +# not depend on the sample data being installed. ``files``/``outdir``/``data_dim``/ |
| 26 | +# ``field_dim`` are set per fit; everything else is a tunable default. Keys are the |
| 27 | +# Fortran param names; friendly aliases are mapped in ``fit``. |
| 28 | +_DEFAULT_PARAMS: dict[str, object] = { |
| 29 | + "block_size": 512, |
| 30 | + "do_opt_block": 0, |
| 31 | + "blk_min": 256, |
| 32 | + "blk_step": 256, |
| 33 | + "blk_max": 1024, |
| 34 | + "num_models": 1, |
| 35 | + "max_threads": 10, |
| 36 | + "use_min_dll": 1, |
| 37 | + "min_dll": 1e-09, |
| 38 | + "use_grad_norm": 1, |
| 39 | + "min_grad_norm": 1e-07, |
| 40 | + "num_mix_comps": 3, |
| 41 | + "pdftype": 0, |
| 42 | + "max_iter": 2000, |
| 43 | + "num_samples": 1, |
| 44 | + "field_blocksize": 1, |
| 45 | + "do_history": 0, |
| 46 | + "histstep": 10, |
| 47 | + "share_comps": 0, |
| 48 | + "share_start": 100, |
| 49 | + "comp_thresh": 0.99, |
| 50 | + "share_iter": 100, |
| 51 | + "lrate": 0.05, |
| 52 | + "minlrate": 1e-08, |
| 53 | + "mineig": 1e-12, |
| 54 | + "lratefact": 0.5, |
| 55 | + "rholrate": 0.05, |
| 56 | + "rho0": 1.5, |
| 57 | + "minrho": 1.0, |
| 58 | + "maxrho": 2.0, |
| 59 | + "rholratefact": 0.5, |
| 60 | + "kurt_start": 3, |
| 61 | + "num_kurt": 5, |
| 62 | + "kurt_int": 1, |
| 63 | + "do_newton": 1, |
| 64 | + "newt_start": 50, |
| 65 | + "newt_ramp": 10, |
| 66 | + "newtrate": 1.0, |
| 67 | + "do_reject": 0, |
| 68 | + "numrej": 3, |
| 69 | + "rejsig": 3.0, |
| 70 | + "rejstart": 2, |
| 71 | + "rejint": 3, |
| 72 | + "writestep": 20, |
| 73 | + "write_nd": 0, |
| 74 | + "write_LLt": 1, |
| 75 | + "decwindow": 1, |
| 76 | + "max_decs": 3, |
| 77 | + "fix_init": 0, |
| 78 | + "update_A": 1, |
| 79 | + "update_c": 1, |
| 80 | + "update_gm": 1, |
| 81 | + "update_alpha": 1, |
| 82 | + "update_mu": 1, |
| 83 | + "update_beta": 1, |
| 84 | + "invsigmax": 100.0, |
| 85 | + "invsigmin": 0.0, |
| 86 | + "do_rho": 1, |
| 87 | + "load_rej": 0, |
| 88 | + "load_W": 0, |
| 89 | + "load_c": 0, |
| 90 | + "load_gm": 0, |
| 91 | + "load_alpha": 0, |
| 92 | + "load_mu": 0, |
| 93 | + "load_beta": 0, |
| 94 | + "load_rho": 0, |
| 95 | + "load_comp_list": 0, |
| 96 | + "do_mean": 1, |
| 97 | + "do_sphere": 1, |
| 98 | + "doPCA": 1, |
| 99 | + "pcakeep": 0, |
| 100 | + "pcadb": 30.0, |
| 101 | + "byte_size": 4, |
| 102 | + "doscaling": 1, |
| 103 | + "scalestep": 1, |
| 104 | +} |
| 105 | + |
| 106 | +# Friendly kwarg -> Fortran param-name aliases (match the Python backends' names). |
| 107 | +_ALIASES = {"n_models": "num_models", "n_mix": "num_mix_comps"} |
| 108 | + |
| 109 | + |
| 110 | +def _render_param(params: dict[str, object]) -> str: |
| 111 | + return "".join(f"{k} {_fmt(v)}\n" for k, v in params.items()) |
| 112 | + |
| 113 | + |
| 114 | +def _fmt(v: object) -> str: |
| 115 | + if isinstance(v, bool): |
| 116 | + return str(int(v)) |
| 117 | + if isinstance(v, float): |
| 118 | + return f"{v:.6e}" if (v != 0 and abs(v) < 1e-3) else f"{v:.6f}" |
| 119 | + return str(v) |
| 120 | + |
| 121 | + |
| 122 | +class AMICANative: |
| 123 | + """Run the native AMICA binary on data and expose the result as an |
| 124 | + ``AmicaOutput``. |
| 125 | +
|
| 126 | + Parameters |
| 127 | + ---------- |
| 128 | + binary : path-like, optional |
| 129 | + Explicit binary path; otherwise resolved for the host (downloaded from the |
| 130 | + release on first use). Equivalent to setting ``PAMICA_NATIVE_BINARY``. |
| 131 | + version : str, default "latest" |
| 132 | + Release tag to resolve the binary from when not given explicitly. |
| 133 | + threads : int, optional |
| 134 | + ``OMP_NUM_THREADS`` for the run (default: the binary's own default). |
| 135 | + timeout : float, optional |
| 136 | + Seconds before the subprocess is killed (default: no timeout). |
| 137 | + **params |
| 138 | + Any Fortran ``input.param`` field (or a friendly alias: ``n_models``, |
| 139 | + ``n_mix``), overriding the defaults; e.g. ``max_iter``, ``lrate``, |
| 140 | + ``pdftype``, ``do_newton``. |
| 141 | + """ |
| 142 | + |
| 143 | + def __init__( |
| 144 | + self, |
| 145 | + binary: Optional[Union[str, Path]] = None, |
| 146 | + *, |
| 147 | + version: str = "latest", |
| 148 | + threads: Optional[int] = None, |
| 149 | + timeout: Optional[float] = None, |
| 150 | + **params: object, |
| 151 | + ) -> None: |
| 152 | + # resolve() now: the subprocess runs with cwd set to a tempdir, so a |
| 153 | + # relative binary path would pass the existence check but fail to launch. |
| 154 | + self.binary = Path(binary).resolve() if binary is not None else None |
| 155 | + self.version = version |
| 156 | + self.threads = threads |
| 157 | + self.timeout = timeout |
| 158 | + self.params = params |
| 159 | + self.output_: Optional[AmicaOutput] = None |
| 160 | + |
| 161 | + def _resolve_binary(self) -> Path: |
| 162 | + if self.binary is not None: |
| 163 | + if not self.binary.exists(): |
| 164 | + raise FileNotFoundError(f"binary not found: {self.binary}") |
| 165 | + return self.binary |
| 166 | + return resolver.resolve(self.version) |
| 167 | + |
| 168 | + def fit(self, X: np.ndarray, **params: object) -> "AMICANative": |
| 169 | + """Run AMICA on ``X`` (shape ``(n_channels, n_samples)``) and store the |
| 170 | + result as ``self.output_`` (an ``AmicaOutput``).""" |
| 171 | + X = np.asarray(X) |
| 172 | + if X.ndim != 2: |
| 173 | + raise ValueError(f"X must be 2-D (n_channels, n_samples); got {X.shape}") |
| 174 | + n_channels, n_samples = X.shape |
| 175 | + |
| 176 | + merged = dict(_DEFAULT_PARAMS) |
| 177 | + for src in (self.params, params): |
| 178 | + for key, value in src.items(): |
| 179 | + merged[_ALIASES.get(key, key)] = value |
| 180 | + merged["data_dim"] = n_channels |
| 181 | + merged["field_dim"] = n_samples |
| 182 | + if not merged.get("pcakeep"): |
| 183 | + merged["pcakeep"] = n_channels # default: keep all components |
| 184 | + |
| 185 | + binary = self._resolve_binary() |
| 186 | + |
| 187 | + with tempfile.TemporaryDirectory(prefix="amica_native_") as td: |
| 188 | + work = Path(td) |
| 189 | + # AMICA reads the data as raw byte_size floats in column-major order |
| 190 | + # (numpy_impl/data.py: reshape order="F"); write it that way. |
| 191 | + byte_size = merged["byte_size"] |
| 192 | + dtype = ( |
| 193 | + np.float32 |
| 194 | + if isinstance(byte_size, int) and byte_size == 4 |
| 195 | + else np.float64 |
| 196 | + ) |
| 197 | + X.astype(dtype).ravel(order="F").tofile(work / "data.fdt") |
| 198 | + |
| 199 | + outdir = work / "amicaout" |
| 200 | + outdir.mkdir() |
| 201 | + # `files` must come first: amica15.f90 hard-stops if it parses other |
| 202 | + # keys before the data file. Dict insertion order preserves that. |
| 203 | + param = {"files": "./data.fdt", "outdir": "./amicaout/", **merged} |
| 204 | + (work / "input.param").write_text(_render_param(param)) |
| 205 | + |
| 206 | + env = None |
| 207 | + if self.threads is not None: |
| 208 | + import os |
| 209 | + |
| 210 | + env = {**os.environ, "OMP_NUM_THREADS": str(self.threads)} |
| 211 | + |
| 212 | + proc = subprocess.run( |
| 213 | + [str(binary), "input.param"], |
| 214 | + cwd=work, |
| 215 | + capture_output=True, |
| 216 | + text=True, |
| 217 | + timeout=self.timeout, |
| 218 | + env=env, |
| 219 | + ) |
| 220 | + if proc.returncode != 0: |
| 221 | + raise RuntimeError( |
| 222 | + f"native AMICA failed (exit {proc.returncode}).\n" |
| 223 | + f"stdout tail:\n{proc.stdout[-2000:]}\n" |
| 224 | + f"stderr tail:\n{proc.stderr[-2000:]}" |
| 225 | + ) |
| 226 | + if not (outdir / "W").exists(): |
| 227 | + raise RuntimeError( |
| 228 | + "native AMICA produced no output (no 'W' file); stdout tail:\n" |
| 229 | + f"{proc.stdout[-2000:]}" |
| 230 | + ) |
| 231 | + # A collapsed fit writes NaN weights (and zero model probabilities); |
| 232 | + # loadmodout's pinv(W@S) would then fail with an opaque SVD error, so |
| 233 | + # detect it here and report it as the degenerate fit it is (cf. the |
| 234 | + # #50 degenerate-fit contract for the Python backends). An empty W |
| 235 | + # (truncated write) is caught too -- an all-NaN check vacuously passes |
| 236 | + # on a zero-length array. |
| 237 | + w_raw = np.fromfile(outdir / "W") |
| 238 | + if w_raw.size == 0 or not np.all(np.isfinite(w_raw)): |
| 239 | + raise RuntimeError( |
| 240 | + "native AMICA produced a degenerate fit (non-finite weights); " |
| 241 | + "the run did not converge. Try more iterations, more data, or " |
| 242 | + "fewer models." |
| 243 | + ) |
| 244 | + self.output_ = loadmodout(outdir) |
| 245 | + return self |
| 246 | + |
| 247 | + def transform(self, X: np.ndarray, model_idx: int = 0) -> np.ndarray: |
| 248 | + """Source activations for ``X`` from the fitted model (delegates to |
| 249 | + ``AmicaOutput.sources``).""" |
| 250 | + if self.output_ is None: |
| 251 | + raise RuntimeError("call fit() before transform().") |
| 252 | + return self.output_.sources(X, model_idx) |
0 commit comments