Phase 15: re-measure parity under the finished epic (#351)

.context/decisions/0003-best-iterate-safeguard.md

@@ -76,10 +76,59 @@ latter joined in Phase 3, issue #289) -- see `pamica/mlx_impl/core.py`'s
 `_snapshot_params`/`_restore_params` and
 `pamica/tests/mlx_tests/test_mlx_keepbest.py`.
 
+## Re-measured after epic #324 (2026-09-23, issue #351)
+
+Epic #324 changed the default trajectory of every backend (component-row
+`doscaling` #333, 1-based schedule gates #335, the reference's iteration order
+and unconditional `A`-freeze #339/#345, the normalized initial `A` #341,
+single-precision constants #344), so the figures above describe code that no
+longer ships. The ensemble was re-run with the #51 protocol
+(`.context/issue-351/keep_best_ensemble.py`: the sample EEG, `n_models=2`, the
+#51 keyword set, N = 20, pamica seeds 0-19), with two changes:
+the reference is the pinned v0.3.3 native binary, seeded 0-19 and
+single-threaded (the #51 runs used the clock-seeded `amica15mac`), and each
+fit runs 300 iterations once, with the 100- and 200-iteration values read from
+its trajectory (checked equal to separate 100-iteration fits, with `keep_best`
+on and off, on two seeds on both sides).
+
+| budget | Fortran mean (sd) | pamica return-last mean (sd) | pamica `keep_best` mean (sd) | sd ratio | pamica minus Fortran | Kolmogorov-Smirnov (KS) p | restores |
+|---:|---:|---:|---:|---:|---:|---:|---:|
+| 100 | -3.3550 (0.0029) | -3.3542 (0.0029) | -3.3542 (0.0029) | 1.0x | +0.0008 | 0.83 | 0 of 20 |
+| 200 | -3.3418 (0.0030) | -3.3416 (0.0022) | -3.3416 (0.0022) | 0.75x | +0.0002 | 0.57 | 0 of 20 |
+| 300 | -3.3392 (0.0027) | -3.3391 (0.0020) | -3.3391 (0.0020) | 0.75x | +0.0001 | 0.98 | 1 of 20 |
+
+- The late overshoots that motivated this decision do not occur in these
+  fits: the largest likelihood decrease in any of the 20 pamica trajectories
+  is 1.2e-5, and none of them used a natural-gradient fallback. At 100
+  iterations the return-last and `keep_best` distributions are identical, and
+  their spread equals the reference's (sd ratio 1.0x, where #51 measured 12.7x
+  and 2.0x). Phase 7 of the epic (#340) traced most of the old overshoots to
+  the column-rule `doscaling` (#333).
+- A restore fired in one of the 20 fits, and only at the 300-iteration budget: seed 3,
+  whose best iterate (iteration 299) beat its last by 4.2e-6.
+- The mean gap is gone at every budget: pamica's mean lies 8.1e-4 above the
+  reference's at 100 iterations, 2.1e-4 at 200 and 1.0e-4 at 300, where #51
+  measured -0.009 at 100 and attributed it to convergence speed. The
+  bundled tier's clock-seeded ensemble (`benchmarks/reproduce_table1.py`,
+  `AMICA` seeds 1-20) agrees: -3.3541 against -3.3543, KS p 0.83. Refitting
+  that ensemble's pamica half with the code before the epic's changes to the fit (e38aa11) against
+  the same 20 reference fits gives -3.3627 (sd 0.006, KS p 1e-5), close to
+  the #27/#51 figures. Seven of those 20 fits stopped early on `min_dll`,
+  whose check counted likelihood dips as small gains until #339 (their mean
+  -3.3679), and the 13 that ran the full 100 iterations average -3.3600, so
+  the old gap came partly from the early stops and partly from the update
+  rule of that code.
+- The decision stands: `keep_best` stays on by default. With these fits it
+  rarely acts, and when it does the gain is small, but a run that ends below
+  its peak (a likelihood decrease near the end, or a Newton fallback) still
+  returns its best iterate.
+
 ## Receipts
 
 - `pamica/torch_impl/amica_torch_ng.py` (`keep_best`, `_snapshot_params`/
   `_restore_params`, `final_ll_`, `_KEEP_BEST_TOL`).
 - `pamica/tests/torch_tests/test_ng_backend.py::test_keep_best_*`.
 - `.context/issue-51/ensemble_ll.py` (Fortran-vs-NG LL ensemble, real data).
+- `.context/issue-351/keep_best_ensemble.py` and `raw/keep_best/` (the
+  re-measurement above), `.context/issue-351/findings.md`.
 - Fortran schedule: `pamica/amica15.f90:1038-1058` (anneal-on-decrease).

.context/issue-27/multimodel_distributional_equivalence.md

@@ -1,5 +1,12 @@
 # Multi-model AMICA parity: distributional equivalence to Fortran (issue #27)
 
+> **Historical record.** These ensembles predate epic #324, which changed every
+> backend's default trajectory (issues #333, #335, #339, #341, #344, #345).
+> The re-measured study under the finished epic, with the same protocol against
+> the pinned v0.3.3 native binary, is in `.context/issue-351/findings.md`
+> (script `.context/issue-351/multimodel_ensemble.py`, which reuses this
+> directory's analysis and figure code).
+
 **Bottom line.** For multi-model AMICA (`n_models > 1`), the natural-gradient
 PyTorch backend (`AMICATorchNG`) is validated against the Fortran reference at
 the level that is actually well-posed: its **ensemble of solutions is

.context/issue-351/_reference_settings.py

@@ -0,0 +1,52 @@
+"""The reference settings of ``benchmarks/reproduce_table1.py``, shared by this
+directory's scripts so that every native-binary run passes its protocol's
+settings explicitly (a change of ``AMICANative``'s defaults cannot move them).
+
+``REFERENCE_SETTINGS`` is read from the script's source without importing it,
+so a script that imports an older pamica tree can use it too.
+"""
+
+from __future__ import annotations
+
+import ast
+import importlib.util
+import sys
+from pathlib import Path
+from typing import Any
+
+_TABLE1 = Path(__file__).resolve().parents[2] / "benchmarks/reproduce_table1.py"
+_NAME = "reproduce_table1_351"
+
+
+def _literal(name: str) -> dict[str, Any]:
+    for node in ast.parse(_TABLE1.read_text()).body:
+        if (
+            isinstance(node, ast.AnnAssign)
+            and isinstance(node.target, ast.Name)
+            and node.target.id == name
+            and node.value is not None
+        ):
+            return dict(ast.literal_eval(node.value))
+    raise LookupError(f"{name} not found in {_TABLE1}")
+
+
+REFERENCE_SETTINGS: dict[str, Any] = _literal("REFERENCE_SETTINGS")
+
+
+def _table1():
+    if _NAME in sys.modules:
+        return sys.modules[_NAME]
+    spec = importlib.util.spec_from_file_location(_NAME, _TABLE1)
+    assert spec is not None and spec.loader is not None
+    mod = importlib.util.module_from_spec(spec)
+    sys.modules[_NAME] = mod
+    spec.loader.exec_module(mod)
+    return mod
+
+
+def single_model_reference_kwargs(max_iter: int) -> dict[str, Any]:
+    return _table1().single_model_reference_kwargs(max_iter)
+
+
+def multimodel_reference_kwargs(max_iter: int) -> dict[str, Any]:
+    return _table1().multimodel_reference_kwargs(max_iter)

.context/issue-351/bundled_single_basins.py

@@ -0,0 +1,114 @@
+"""Basins of the bundled tier's single-model comparison (issue #351 item 2).
+
+``benchmarks/reproduce_table1.py --tier bundled`` compares five pamica fits
+(seeds 301-305) with five clock-seeded reference fits, one pair per seed, on
+the bundled sample (32 channels, k ~ 30), Newton off, 2000 iterations, and
+reports only the means. This script repeats that single-model protocol with
+the unmixing matrices kept, so every pair can be compared: the same pamica
+seeds (deterministic, so they reproduce the tier's pamica fits) and seeded
+reference runs (``seed`` 1..N, reproducible), with the tier's keyword sets.
+
+    uv run python .context/issue-351/bundled_single_basins.py OUT.npz [N_FORTRAN] [--threads T]
+"""
+
+from __future__ import annotations
+
+import argparse
+import importlib.util
+import itertools
+import json
+import sys
+from pathlib import Path
+
+import numpy as np
+
+HERE = Path(__file__).resolve().parent
+REPO = HERE.parents[1]
+
+
+def _table1():
+    spec = importlib.util.spec_from_file_location(
+        "reproduce_table1_basins", REPO / "benchmarks/reproduce_table1.py"
+    )
+    assert spec is not None and spec.loader is not None
+    mod = importlib.util.module_from_spec(spec)
+    sys.modules["reproduce_table1_basins"] = mod
+    spec.loader.exec_module(mod)
+    return mod
+
+
+def fortran_kwargs(t1, seed: int, threads: int) -> dict:
+    """``AMICANative`` keywords of the reference runs: the tier's single-model
+    settings, every one explicit, plus the seed."""
+    return dict(
+        threads=threads, max_threads=threads, timeout=3600, seed=seed,
+        **t1.single_model_reference_kwargs(2000),
+    )  # fmt: skip
+
+
+def main() -> None:
+    p = argparse.ArgumentParser()
+    p.add_argument("out", type=Path)
+    p.add_argument("n_fortran", nargs="?", type=int, default=10)
+    p.add_argument("--threads", type=int, default=4)
+    p.add_argument("--pamica-seeds", default="301,302,303,304,305")
+    a = p.parse_args()
+    t1 = _table1()
+    from pamica import AMICA, AMICANative
+
+    data, _ = t1.load_bundled_data()
+    binary = t1.resolve_binary(None, "v0.3.3")
+    runs: dict[str, dict] = {}
+    for s in range(1, a.n_fortran + 1):
+        eng = AMICANative(binary=binary, **fortran_kwargs(t1, s, a.threads))
+        eng.fit(data)
+        assert eng.output_ is not None
+        runs[f"fortran_seed{s}"] = {
+            "W": eng.output_.W[:, :, 0],
+            "ll": float(eng.output_.LL[-1]),
+        }
+        print(f"fortran seed {s}: LL={runs[f'fortran_seed{s}']['ll']:.6f}", flush=True)
+    for s in [int(x) for x in a.pamica_seeds.split(",") if x]:
+        model = AMICA(n_models=1, n_mix=3, device="cpu", verbose=False)
+        model.fit(
+            data, max_iter=2000, lrate=0.05, do_mean=True, do_sphere=True,
+            do_approx_sphere=True, do_newton=False, seed=s, block_size=512,
+            minlrate=1e-8, lratefact=0.5, maxdecs=3, newt_start=50, newt_ramp=10,
+            newtrate=1.0, rho0=1.5, minrho=1.0, maxrho=2.0, rholrate=0.05,
+            rholratefact=0.5, invsigmin=0.0, invsigmax=100.0, doscaling=True,
+            scalestep=1,
+        )  # fmt: skip
+        runs[f"pamica_seed{s}"] = {
+            "W": model.get_unmixing_matrix(0),
+            "ll": float(model.final_ll_),
+            "iters": len(model.ll_history_),
+            "stop": str(model.stop_reason_),
+        }
+        print(
+            f"pamica seed {s}: LL={model.final_ll_:.6f} iters={len(model.ll_history_)} "
+            f"stop={model.stop_reason_}",
+            flush=True,
+        )
+    names = list(runs)
+    np.savez(
+        a.out,
+        names=np.array(names),
+        W=np.stack([runs[k]["W"] for k in names]),
+        ll=np.array([runs[k]["ll"] for k in names]),
+    )
+    rep: dict = {"ll": {k: runs[k]["ll"] for k in names}, "pairs": {}}
+    for x, y in itertools.combinations(names, 2):
+        c = t1.xcorr(runs[x]["W"], runs[y]["W"])
+        rep["pairs"][f"{x} vs {y}"] = {
+            "mean_corr": float(c.mean()),
+            "min_corr": float(c.min()),
+            "amari": float(t1.amari_distance(runs[x]["W"], runs[y]["W"])),
+        }
+    for k in ("iters", "stop"):
+        rep[k] = {n: runs[n][k] for n in names if k in runs[n]}
+    a.out.with_suffix(".json").write_text(json.dumps(rep, indent=2))
+    print(json.dumps(rep["ll"], indent=2))
+
+
+if __name__ == "__main__":
+    main()

.context/issue-351/equivalence_check.py

@@ -0,0 +1,100 @@
+"""The ±0.05 equivalence margin of the #27 multi-model study, on saved ensembles
+(issue #351 item 3). No refitting.
+
+The #27 record states its acceptance criterion as "between not worse than
+within-Fortran; TOST within a margin" of ±0.05 on the mean Hungarian-matched
+correlation. Its first TOST was computed on the 190/400 pairwise values as if
+they were independent, which the record later retired as pseudoreplicated (each
+run appears in ~39 pairs) in favor of the run-level permutation test. For each
+ensemble given, this script reports, for the correlation and the Amari distance:
+
+* ``diff``: mean(between) - mean(within-Fortran), and whether it lies inside
+  ±0.05 (the documented descriptive check, stated for the correlation only);
+* ``pairwise_tost_p``: the retired pairwise TOST with the same margin, for
+  continuity with the old record (its p-value is not valid, as the record says);
+* ``run_bootstrap_90ci``: a run-level bootstrap 90% interval of ``diff``
+  (runs resampled with replacement within each implementation, 5000 draws),
+  which respects the shared-run dependence; an interval inside ±0.05 is the
+  run-level analogue of the TOST at the 5% level.
+
+    uv run python .context/issue-351/equivalence_check.py OUT.json NAME=PATH.npz [NAME=PATH.npz ...]
+"""
+
+from __future__ import annotations
+
+import importlib.util
+import json
+import sys
+from pathlib import Path
+
+import numpy as np
+from scipy import stats
+
+HERE = Path(__file__).resolve().parent
+REPO = HERE.parents[1]
+DELTA = 0.05
+
+
+def _load(path: Path, name: str):
+    spec = importlib.util.spec_from_file_location(name, path)
+    assert spec is not None and spec.loader is not None
+    mod = importlib.util.module_from_spec(spec)
+    sys.modules[name] = mod
+    spec.loader.exec_module(mod)
+    return mod
+
+
+def analyze(Fs: np.ndarray, Gs: np.ndarray, metric) -> dict:
+    allW = np.concatenate([Fs, Gs])
+    m, n = len(allW), len(Fs)
+    P = np.zeros((m, m))
+    for i in range(m):
+        for j in range(i + 1, m):
+            P[i, j] = P[j, i] = metric(allW[i], allW[j])
+    F = np.arange(n)
+    G = np.arange(n, m)
+    wF = P[np.ix_(F, F)][np.triu_indices(n, 1)]
+    bt = P[np.ix_(G, F)].ravel()
+    diff = float(bt.mean() - wF.mean())
+    se = np.sqrt(bt.var(ddof=1) / bt.size + wF.var(ddof=1) / wF.size)
+    p_tost = float(
+        max(stats.norm.sf((diff + DELTA) / se), stats.norm.cdf((diff - DELTA) / se))
+    )
+    rng = np.random.default_rng(0)
+    boot = []
+    for _ in range(5000):
+        f = rng.choice(F, n, replace=True)
+        g = rng.choice(G, len(G), replace=True)
+        # within-Fortran over distinct draws only (a run is not paired with itself)
+        sub = P[np.ix_(f, f)]
+        mask = f[:, None] != f[None, :]
+        w = sub[np.triu(mask, 1)].mean()
+        b = P[np.ix_(g, f)].mean()
+        boot.append(b - w)
+    lo, hi = np.percentile(boot, [5, 95])
+    return {
+        "diff_between_minus_within_fortran": diff,
+        "inside_pm_0.05": bool(abs(diff) < DELTA),
+        "pairwise_tost_p_retired": p_tost,
+        "run_bootstrap_90ci": [float(lo), float(hi)],
+        "run_bootstrap_90ci_inside_pm_0.05": bool(lo > -DELTA and hi < DELTA),
+    }
+
+
+def main() -> None:
+    out = Path(sys.argv[1])
+    ama27 = _load(REPO / ".context/issue-27/amari_distance.py", "ama27_eq")
+    rep = {}
+    for arg in sys.argv[2:]:
+        name, path = arg.split("=", 1)
+        d = np.load(path)
+        rep[name] = {
+            "corr": analyze(d["Fs"], d["Gs"], ama27.xcorr),
+            "amari": analyze(d["Fs"], d["Gs"], ama27.model_amari),
+        }
+        print(name, json.dumps(rep[name], indent=1), flush=True)
+    out.write_text(json.dumps(rep, indent=2))
+
+
+if __name__ == "__main__":
+    main()