Address JOSS reviewer feedback on paper and Figure 1

.context/issue-144-parity-data-adequacy/bundled_sample_newton0.py

@@ -0,0 +1,210 @@
+"""do_newton=0 parity + Fortran-vs-Fortran self-consistency on the bundled
+32-channel sample EEG. This is the source of Table 1's bundled-sample Amari
+distance row and the score-function/multi-model checks, which paper.md still
+describes as needing no external download; Table 1's single-model headline
+correlation/LL instead comes from the external ds002718 recording (see
+run_5seed_newton0.sh), since the bundled sample sits at the project's own
+data-adequacy boundary (k~30).
+
+Runs N seeds of: Fortran (amica15mac, do_newton overridden to 0) and
+AMICATorchNG (do_newton=False), then reports:
+  - mean/min Hungarian-matched correlation + Amari distance, NG vs Fortran
+  - mean/min Hungarian-matched correlation + Amari distance, Fortran vs Fortran
+
+    uv run python bundled_sample_newton0.py [n_seeds] [max_iter] [out_dir]
+"""
+
+import itertools
+import json
+import shutil
+import subprocess
+import sys
+import tempfile
+from pathlib import Path
+
+import numpy as np
+import torch
+from scipy.optimize import linear_sum_assignment
+
+HERE = Path(__file__).resolve().parent
+REPO = HERE.parents[1]
+sys.path.insert(0, str(REPO))
+from pyAMICA import AMICA  # noqa: E402
+from pyAMICA.torch_impl.utils import load_eeglab_data  # noqa: E402
+
+
+def xcorr(Wa, Wb):
+    na = Wa / (np.linalg.norm(Wa, axis=1, keepdims=True) + 1e-12)
+    nb = Wb / (np.linalg.norm(Wb, axis=1, keepdims=True) + 1e-12)
+    corr = np.abs(na @ nb.T)
+    r, c = linear_sum_assignment(1 - corr)
+    return corr[r, c]
+
+
+def amari_index(gain):
+    n = gain.shape[0]
+    abs_gain = np.abs(gain)
+    row_max = abs_gain.max(axis=1)
+    col_max = abs_gain.max(axis=0)
+    if np.any(row_max == 0) or np.any(col_max == 0):
+        raise ValueError("amari_index: a row or column is all-zero (degenerate W)")
+    row_term = (abs_gain.sum(axis=1) / row_max - 1).sum()
+    col_term = (abs_gain.sum(axis=0) / col_max - 1).sum()
+    return (row_term + col_term) / (2 * n * (n - 1))
+
+
+def amari_distance(Wa, Wb):
+    forward = amari_index(Wa @ np.linalg.pinv(Wb))
+    backward = amari_index(Wb @ np.linalg.pinv(Wa))
+    return float((forward + backward) / 2)
+
+
+def run_fortran(data_dim, max_iter, seed, work_dir):
+    work_dir.mkdir(parents=True, exist_ok=True)
+    shutil.copy(
+        REPO / "pyAMICA/sample_data/eeglab_data.fdt", work_dir / "eeglab_data.fdt"
+    )
+    template = (REPO / "pyAMICA/sample_data/input.param").read_text().splitlines()
+    lines = []
+    for line in template:
+        if line.startswith("files"):
+            lines.append("files ./eeglab_data.fdt")
+        elif line.startswith("outdir"):
+            lines.append("outdir ./fortran_output/")
+        elif line.startswith("max_iter"):
+            lines.append(f"max_iter {max_iter}")
+        elif line.startswith("do_newton"):
+            lines.append("do_newton 0")
+        else:
+            lines.append(line)
+    (work_dir / "input.param").write_text("\n".join(lines) + "\n")
+    (work_dir / "fortran_output").mkdir(exist_ok=True)
+
+    result = subprocess.run(
+        [str(REPO / "pyAMICA/sample_data/amica15mac"), "input.param"],
+        cwd=work_dir,
+        capture_output=True,
+        text=True,
+        timeout=600,
+    )
+    if result.returncode != 0:
+        raise RuntimeError(f"Fortran failed (seed {seed}): {result.stderr}")
+
+    W = np.fromfile(work_dir / "fortran_output/W", dtype=np.float64).reshape(
+        data_dim, data_dim, order="F"
+    )
+    ll = None
+    out_txt = work_dir / "fortran_output/out.txt"
+    if out_txt.exists():
+        for line in out_txt.read_text().splitlines():
+            if line.strip().startswith("iter"):
+                parts = line.split()
+                if "LL" in parts:
+                    ll = float(parts[parts.index("LL") + 2])
+    if ll is None:
+        print(f"seed {seed}: WARNING could not parse LL from Fortran out.txt")
+    return W, ll
+
+
+def main():
+    n_seeds = int(sys.argv[1]) if len(sys.argv) > 1 else 5
+    max_iter = int(sys.argv[2]) if len(sys.argv) > 2 else 2000
+    out_dir = (
+        Path(sys.argv[3])
+        if len(sys.argv) > 3
+        else Path(tempfile.gettempdir()) / "pyamica_newton0_bundled"
+    )
+    seeds = list(range(301, 301 + n_seeds))
+
+    with open(REPO / "pyAMICA/sample_data/sample_params.json") as f:
+        params = json.load(f)
+    data_dim = params["data_dim"]
+    field_dim = params["field_dim"][0]
+
+    data = load_eeglab_data(
+        str(REPO / "pyAMICA/sample_data/eeglab_data.fdt"),
+        data_dim=data_dim,
+        field_dim=field_dim,
+        dtype=np.float32,
+    ).astype(np.float64)
+
+    tmp_root = out_dir
+    tmp_root.mkdir(parents=True, exist_ok=True)
+
+    fortran_Ws = {}
+    ng_corrs = []
+    ng_amaris = []
+    for seed in seeds:
+        work_dir = tmp_root / f"seed{seed}"
+        W_f, ll_f = run_fortran(data_dim, max_iter, seed, work_dir)
+        fortran_Ws[seed] = W_f
+        print(f"seed {seed}: Fortran done, LL={ll_f}", flush=True)
+
+        np.random.seed(seed)
+        torch.manual_seed(seed)
+        model = AMICA(n_models=1, n_mix=params.get("num_mix", 3), verbose=False)
+        model.fit(
+            data,
+            max_iter=max_iter,
+            lrate=params.get("lrate", 0.05),
+            do_mean=params.get("do_mean", True),
+            do_sphere=params.get("do_sphere", True),
+            do_approx_sphere=params.get("do_approx_sphere", True),
+            do_newton=False,
+            seed=seed,
+        )
+        W_ng = model.get_unmixing_matrix(0)
+
+        corrs = xcorr(W_f, W_ng)
+        amari = amari_distance(W_f, W_ng)
+        ng_corrs.append(corrs)
+        ng_amaris.append(amari)
+        print(
+            f"seed {seed}: NG vs Fortran mean_corr={corrs.mean():.4f} "
+            f"min={corrs.min():.4f} amari={amari:.4f} ng_ll={model.final_ll_:.4f}",
+            flush=True,
+        )
+
+    all_means = np.array([c.mean() for c in ng_corrs])
+    all_mins = np.array([c.min() for c in ng_corrs])
+    print(
+        f"\nNG-vs-Fortran ({n_seeds} seeds): mean_corr={all_means.mean():.4f} "
+        f"(sd {all_means.std():.4f}), min_corr overall={all_mins.min():.4f}, "
+        f"mean_amari={np.mean(ng_amaris):.4f}",
+        flush=True,
+    )
+
+    ff_means = []
+    ff_mins = []
+    ff_amaris = []
+    for a, b in itertools.combinations(seeds, 2):
+        corrs = xcorr(fortran_Ws[a], fortran_Ws[b])
+        amari = amari_distance(fortran_Ws[a], fortran_Ws[b])
+        ff_means.append(corrs.mean())
+        ff_mins.append(corrs.min())
+        ff_amaris.append(amari)
+        print(
+            f"fortran seed {a} vs {b}: mean_corr={corrs.mean():.4f} "
+            f"min={corrs.min():.4f} amari={amari:.4f}",
+            flush=True,
+        )
+
+    if not ff_means:
+        print(
+            "\nFortran-vs-Fortran: skipped (need at least 2 seeds for a pair)",
+            flush=True,
+        )
+        return
+
+    ff_means = np.array(ff_means)
+    ff_mins = np.array(ff_mins)
+    print(
+        f"\nFortran-vs-Fortran ({n_seeds} seeds, {len(ff_means)} pairs): "
+        f"mean_corr={ff_means.mean():.4f} (sd {ff_means.std():.4f}), "
+        f"min_corr overall={ff_mins.min():.4f}, mean_amari={np.mean(ff_amaris):.4f}",
+        flush=True,
+    )
+
+
+if __name__ == "__main__":
+    main()

.context/issue-144-parity-data-adequacy/run_5seed_newton0.sh

@@ -0,0 +1,27 @@
+#!/bin/bash
+# 5-seed Fortran-vs-NG comparison with Newton disabled (do_newton=0), forced
+# full 2000-iteration budget (use_min_dll/use_grad_norm off), on the full
+# ds002718 sub-002 recording (70ch, 747750 frames, k=152.6).
+#
+# Pipelined: each seed's Fortran phase (CPU, 24 threads) runs strictly
+# sequentially (splitting threads across concurrent Fortran jobs was
+# catastrophically slow in an earlier attempt), but as soon as a seed's
+# Fortran finishes, its NG/CUDA phase (GPU) is launched in the background and
+# the NEXT seed's Fortran starts immediately -- CPU and GPU work overlap
+# since they don't compete for the same resource.
+set -e
+cd ~/pyAMICA-issue144
+NPY=benchmarks/data/ds002718_sub-002_eeg70_full.npy
+SCRIPTS=.context/issue-144-parity-data-adequacy
+
+for s in 201 202 203 204 205; do
+  echo "=== seed $s: Fortran (24 threads, do_newton=0) ==="
+  uv run python -u $SCRIPTS/run_fortran_only.py "$NPY" "$s" 2000 24 /tmp/newton0_seed$s 0
+  echo "=== seed $s: launching NG/CUDA in background, moving to next seed's Fortran ==="
+  uv run python -u $SCRIPTS/run_ng_only.py "$NPY" "$s" 2000 /tmp/newton0_seed$s 0 \
+    > /tmp/newton0_seed${s}_ng.log 2>&1 &
+done
+
+echo "=== all Fortran phases launched/done, waiting on remaining background NG jobs ==="
+wait
+echo "=== ALL DONE ==="

.context/issue-144-parity-data-adequacy/run_fortran_only.py

@@ -0,0 +1,111 @@
+"""Fortran-only phase for one seed (CPU-bound, safe to run concurrently with
+another seed's CUDA phase or other Fortran-only jobs at reduced thread counts).
+Writes W + LL to a fixed per-seed directory for a later, separate NG phase to
+pick up.
+
+    uv run python run_fortran_only.py <npy> <seed> <max_iter> <threads> <out_dir> [do_newton]
+"""
+
+import os
+import subprocess
+import sys
+from pathlib import Path
+
+import numpy as np
+
+HERE = Path(__file__).resolve().parent
+REPO = HERE.parents[1]
+BIN = REPO / "pyAMICA/sample_data/amica15_linux"
+INPUT_PARAM = REPO / "pyAMICA/sample_data/input.param"
+
+
+def write_fdt(data: np.ndarray, path: Path) -> None:
+    path.write_bytes(np.ascontiguousarray(data).astype("<f4").tobytes(order="F"))
+
+
+def main():
+    npy_path = Path(sys.argv[1])
+    seed = int(sys.argv[2])
+    max_iter = int(sys.argv[3])
+    threads = int(sys.argv[4])
+    out_dir = Path(sys.argv[5])
+    do_newton = int(sys.argv[6]) if len(sys.argv) > 6 else 1
+    out_dir.mkdir(parents=True, exist_ok=True)
+
+    data = np.load(npy_path).astype(np.float64)
+    nw, field = data.shape
+
+    fdt_path = out_dir / "data.fdt"
+    if not fdt_path.exists():
+        write_fdt(data, fdt_path)
+
+    (out_dir / "fortran_output").mkdir(parents=True, exist_ok=True)
+    lines = []
+    for ln in INPUT_PARAM.read_text().splitlines():
+        if ln.startswith("files"):
+            lines.append("files ./data.fdt")
+        elif ln.startswith("outdir"):
+            lines.append("outdir ./fortran_output/")
+        elif ln.startswith("data_dim"):
+            lines.append(f"data_dim {nw}")
+        elif ln.startswith("field_dim"):
+            lines.append(f"field_dim {field}")
+        elif ln.startswith("max_iter"):
+            lines.append(f"max_iter {max_iter}")
+        elif ln.startswith("max_threads"):
+            lines.append(f"max_threads {threads}")
+        elif ln.startswith("pcakeep"):
+            lines.append(f"pcakeep {nw}")
+        elif ln.startswith("use_min_dll"):
+            # Force the full max_iter budget instead of Fortran's own
+            # early-stopping (matches benchmark_dimsweep.py's convention for
+            # fixed-length, matched runs): otherwise Fortran can converge and
+            # stop well short of 2000 while AMICATorchNG (no early-stopping
+            # equivalent) keeps optimizing past that point, letting weakly-
+            # determined components drift to a different, still-valid optimum
+            # -- an asymmetry, not real disagreement.
+            lines.append("use_min_dll 0")
+        elif ln.startswith("use_grad_norm"):
+            lines.append("use_grad_norm 0")
+        elif ln.startswith("do_newton"):
+            lines.append(f"do_newton {do_newton}")
+        else:
+            lines.append(ln)
+    (out_dir / "input.param").write_text("\n".join(lines) + "\n")
+
+    orig = os.getcwd()
+    os.chdir(out_dir)
+    env = {**os.environ, "OMP_NUM_THREADS": str(threads)}
+    print(
+        f"seed {seed}: starting Fortran ({threads} threads, do_newton={do_newton})...",
+        flush=True,
+    )
+    try:
+        r = subprocess.run(
+            [str(BIN), "input.param"], capture_output=True, text=True,
+            timeout=3600, env=env,
+        )  # fmt: skip
+    finally:
+        os.chdir(orig)
+    if r.returncode != 0:
+        print(f"seed {seed}: FAILED: {r.stderr[-500:]}", flush=True)
+        sys.exit(1)
+    fort_ll = next(
+        (
+            float(ln.split("LL =")[1].split()[0])
+            for ln in reversed(r.stdout.splitlines())
+            if "LL =" in ln
+        ),
+        None,
+    )
+    if fort_ll is None:
+        print(
+            f"seed {seed}: WARNING could not parse LL from Fortran stdout", flush=True
+        )
+        fort_ll = float("nan")
+    (out_dir / "fortran_ll.txt").write_text(str(fort_ll))
+    print(f"seed {seed}: Fortran done, LL={fort_ll:.4f}", flush=True)
+
+
+if __name__ == "__main__":
+    main()