diff --git a/.context/issue-144-parity-data-adequacy/bundled_sample_newton0.py b/.context/issue-144-parity-data-adequacy/bundled_sample_newton0.py new file mode 100644 index 0000000..72b4819 --- /dev/null +++ b/.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() diff --git a/.context/issue-144-parity-data-adequacy/run_5seed_newton0.sh b/.context/issue-144-parity-data-adequacy/run_5seed_newton0.sh new file mode 100644 index 0000000..54aff08 --- /dev/null +++ b/.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 ===" diff --git a/.context/issue-144-parity-data-adequacy/run_fortran_only.py b/.context/issue-144-parity-data-adequacy/run_fortran_only.py new file mode 100644 index 0000000..3924da8 --- /dev/null +++ b/.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 [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(" 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() diff --git a/.context/issue-144-parity-data-adequacy/run_ng_only.py b/.context/issue-144-parity-data-adequacy/run_ng_only.py new file mode 100644 index 0000000..448af6c --- /dev/null +++ b/.context/issue-144-parity-data-adequacy/run_ng_only.py @@ -0,0 +1,72 @@ +"""NG (CUDA) phase + correlation for one seed, picking up a Fortran W already +written by run_fortran_only.py to the same out_dir. Run these sequentially +(one GPU) after the Fortran-only jobs (which can run concurrently on CPU). + + uv run python run_ng_only.py [do_newton] +""" + +import sys +from pathlib import Path + +import numpy as np +from scipy.optimize import linear_sum_assignment + +HERE = Path(__file__).resolve().parent +REPO = HERE.parents[1] +sys.path.insert(0, str(REPO)) +from pyAMICA.torch_impl import AMICATorchNG # 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 main(): + npy_path = Path(sys.argv[1]) + seed = int(sys.argv[2]) + max_iter = int(sys.argv[3]) + out_dir = Path(sys.argv[4]) + do_newton = bool(int(sys.argv[5])) if len(sys.argv) > 5 else True + + data = np.load(npy_path).astype(np.float64) + nw, field = data.shape + + W_fortran = np.fromfile(out_dir / "fortran_output/W", dtype=np.float64).reshape( + nw, nw, order="F" + ) + fort_ll = float((out_dir / "fortran_ll.txt").read_text()) + + print( + f"seed {seed}: running AMICATorchNG on CUDA (do_newton={do_newton})...", + flush=True, + ) + m = AMICATorchNG( + n_channels=nw, n_models=1, n_mix=3, block_size=512, lrate=0.05, + minlrate=1e-8, lratefact=0.5, maxdecs=3, do_newton=do_newton, + 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, seed=seed, device="cuda", + ) # fmt: skip + m.fit(data, max_iter=max_iter, verbose=False) + if m.stop_reason in AMICATorchNG._DEGENERATE_STOP_REASONS: + raise RuntimeError( + f"seed {seed}: NG fit ended degenerate (stop_reason={m.stop_reason!r}); " + "refusing to report a correlation against a singular/NaN W" + ) + W_ng = m.get_unmixing_matrix(0) + + corrs = xcorr(W_fortran, W_ng) + n_above_95 = int((corrs > 0.95).sum()) + print( + f"seed {seed}: mean_corr={corrs.mean():.4f} min={corrs.min():.4f} " + f"n_above_0.95={n_above_95}/{nw} fortran_LL={fort_ll:.4f} ng_LL={m.final_ll_:.4f}", + flush=True, + ) + + +if __name__ == "__main__": + main() diff --git a/.context/issue-144-parity-data-adequacy/test_ds002718_32ch.py b/.context/issue-144-parity-data-adequacy/test_ds002718_32ch.py new file mode 100644 index 0000000..534ea9c --- /dev/null +++ b/.context/issue-144-parity-data-adequacy/test_ds002718_32ch.py @@ -0,0 +1,153 @@ +"""Single-model Fortran-vs-NG parity on a data-adequate real recording (k>=60 +per issue #90's documented threshold), using the bundled amica15mac binary +(also used for Table 1's bundled-sample Amari row; Table 1's headline +correlation instead uses the Linux amica15_linux build via +run_5seed_newton0.sh) and AMICATorchNG, instead of the under-determined +bundled 32ch/30504-frame EEGLAB tutorial recording (k=29.8). + +Data: OpenNeuro ds002718 sub-002 (Wakeman-Henson faces), first 32 of the first +70 EEG channels, first N frames -- real data, not committed (not bundled; see +benchmarks/README_dimsweep.md for the download recipe this mirrors). + + uv run python .context/issue-144-parity-data-adequacy/test_ds002718_32ch.py [n_seeds] [max_iter] +""" + +import os +import shutil +import subprocess +import sys +import tempfile +from pathlib import Path + +import numpy as np +from scipy.optimize import linear_sum_assignment + +HERE = Path(__file__).resolve().parent +REPO = HERE.parents[1] +sys.path.insert(0, str(REPO)) +from pyAMICA.torch_impl import AMICATorchNG # noqa: E402 + +BIN = REPO / "pyAMICA/sample_data/amica15mac" +INPUT_PARAM = REPO / "pyAMICA/sample_data/input.param" + + +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 main(): + npy_path = Path(sys.argv[1]) + n_seeds = int(sys.argv[2]) if len(sys.argv) > 2 else 1 + max_iter = int(sys.argv[3]) if len(sys.argv) > 3 else 300 + + data = np.load(npy_path).astype(np.float64) + nw, field = data.shape + k = field / nw**2 + print(f"data: {nw} channels x {field} frames, k={k:.1f}, max_iter={max_iter}") + + work = Path(tempfile.mkdtemp(prefix="ds002718_parity_")) + fdt_path = work / "data.fdt" + data.astype(np.float32).T.tofile(fdt_path) + + results = [] + for seed in range(n_seeds): + d = work / f"run_{seed}" + (d / "fortran_output").mkdir(parents=True, exist_ok=True) + shutil.copy(fdt_path, d / "data.fdt") + 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("pcakeep"): + # AMICATorchNG has no PCA source reduction (n_sources == + # n_channels always), so pcakeep must track nw or Fortran's + # W comes back reduced-rank (the template's literal 32 is + # only correct by coincidence at nw=32). + lines.append(f"pcakeep {nw}") + elif ln.startswith("use_min_dll"): + # Force the full max_iter budget instead of Fortran's own + # early-stopping: otherwise Fortran can converge and stop + # well short of max_iter while AMICATorchNG (no equivalent) + # keeps optimizing past that point, drifting weakly- + # determined components to a different, still-valid optimum. + lines.append("use_min_dll 0") + elif ln.startswith("use_grad_norm"): + lines.append("use_grad_norm 0") + else: + lines.append(ln) + (d / "input.param").write_text("\n".join(lines) + "\n") + + orig = os.getcwd() + os.chdir(d) + try: + r = subprocess.run( + [str(BIN), "input.param"], capture_output=True, text=True, timeout=1800 + ) + finally: + os.chdir(orig) + if r.returncode != 0: + print(f"seed {seed}: Fortran failed: {r.stderr[-400:]}") + continue + W_fortran = np.fromfile(d / "fortran_output/W", dtype=np.float64).reshape( + nw, nw, order="F" + ) + 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") + fort_ll = float("nan") + + m = AMICATorchNG( + n_channels=nw, n_models=1, n_mix=3, block_size=512, lrate=0.05, + minlrate=1e-8, lratefact=0.5, maxdecs=3, do_newton=True, + 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, seed=seed, device="cpu", + ) # fmt: skip + m.fit(data, max_iter=max_iter, verbose=False) + if m.stop_reason in AMICATorchNG._DEGENERATE_STOP_REASONS: + print( + f"seed {seed}: NG fit ended degenerate (stop_reason={m.stop_reason!r}), skipping" + ) + continue + W_ng = m.get_unmixing_matrix(0) + + corrs = xcorr(W_fortran, W_ng) + print( + f"seed {seed}: mean_corr={corrs.mean():.4f} min={corrs.min():.4f} " + f"fortran_LL={fort_ll:.4f} ng_LL={m.final_ll_:.4f}" + ) + results.append(corrs.mean()) + + if results: + print( + f"\n{nw}ch/{field}fr (k={k:.1f}), n={len(results)}: " + f"mean={np.mean(results):.4f} range={np.min(results):.4f}-{np.max(results):.4f}" + ) + + if len(results) < n_seeds: + print(f"\n{n_seeds - len(results)}/{n_seeds} seed(s) failed or were degenerate") + return 1 + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/.context/issue-144-parity-data-adequacy/test_ds002718_hallu.py b/.context/issue-144-parity-data-adequacy/test_ds002718_hallu.py new file mode 100644 index 0000000..e2333ce --- /dev/null +++ b/.context/issue-144-parity-data-adequacy/test_ds002718_hallu.py @@ -0,0 +1,181 @@ +"""Single-model Fortran-vs-NG parity on the FULL ds002718 sub-002 recording +(70 channels, 747750 frames, k=152.6 -- the exact configuration issue #90's +k-sweep already validated at ~0.98 mean |corr|, native-fortran-f64 vs +torch-cuda-f64 specifically at 0.995), run on hallu: native Linux Fortran +build (24 threads) + PyTorch CUDA, across multiple independent seeds. + + uv run python .context/issue-144-parity-data-adequacy/test_ds002718_hallu.py [n_seeds] [max_iter] [threads] [seed_offset] +""" + +import os +import subprocess +import sys +import tempfile +from pathlib import Path + +import numpy as np +from scipy.optimize import linear_sum_assignment + +HERE = Path(__file__).resolve().parent +REPO = HERE.parents[1] +sys.path.insert(0, str(REPO)) +from pyAMICA.torch_impl import AMICATorchNG # noqa: E402 + +BIN = REPO / "pyAMICA/sample_data/amica15_linux" +INPUT_PARAM = REPO / "pyAMICA/sample_data/input.param" + + +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 write_fdt(data: np.ndarray, path: Path) -> None: + """(n_channels, n_samples) -> amica's raw float32 .fdt, channel-fastest + (column-major) order -- byte-identical to EEGLAB's own .fdt layout.""" + path.write_bytes(np.ascontiguousarray(data).astype(" 2 else 1 + max_iter = int(sys.argv[3]) if len(sys.argv) > 3 else 2000 + threads = int(sys.argv[4]) if len(sys.argv) > 4 else 24 + seed_offset = int(sys.argv[5]) if len(sys.argv) > 5 else 0 + + data = np.load(npy_path).astype(np.float64) + nw, field = data.shape + k = field / nw**2 + print( + f"data: {nw} channels x {field} frames, k={k:.1f}, max_iter={max_iter}, " + f"threads={threads}", + flush=True, + ) + + work = Path(tempfile.mkdtemp(prefix="ds002718_full_parity_")) + fdt_path = work / "data.fdt" + write_fdt(data, fdt_path) + + results = [] + for seed in range(seed_offset, seed_offset + n_seeds): + d = work / f"run_{seed}" + (d / "fortran_output").mkdir(parents=True, exist_ok=True) + os.symlink(fdt_path, d / "data.fdt") + 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"): + # AMICATorchNG has no PCA source reduction (n_sources == + # n_channels always), so pcakeep must track nw or Fortran's + # W comes back reduced-rank (the template's literal 32 is + # only correct by coincidence at nw=32). + 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): + # otherwise Fortran can converge and stop well short of + # max_iter 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") + else: + lines.append(ln) + (d / "input.param").write_text("\n".join(lines) + "\n") + + orig = os.getcwd() + os.chdir(d) + env = {**os.environ, "OMP_NUM_THREADS": str(threads)} + try: + print( + f"seed {seed}: running Fortran (native, {threads} threads)...", + flush=True, + ) + r = subprocess.run( + [str(BIN), "input.param"], capture_output=True, text=True, + timeout=1800, env=env, + ) # fmt: skip + finally: + os.chdir(orig) + if r.returncode != 0: + print(f"seed {seed}: Fortran failed: {r.stderr[-500:]}", flush=True) + continue + W_fortran = np.fromfile(d / "fortran_output/W", dtype=np.float64).reshape( + nw, nw, order="F" + ) + 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") + + print(f"seed {seed}: running AMICATorchNG on CUDA...", flush=True) + m = AMICATorchNG( + n_channels=nw, n_models=1, n_mix=3, block_size=512, lrate=0.05, + minlrate=1e-8, lratefact=0.5, maxdecs=3, do_newton=True, + 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, seed=seed, device="cuda", + ) # fmt: skip + m.fit(data, max_iter=max_iter, verbose=False) + if m.stop_reason in AMICATorchNG._DEGENERATE_STOP_REASONS: + print( + f"seed {seed}: NG fit ended degenerate (stop_reason={m.stop_reason!r}), skipping", + flush=True, + ) + continue + W_ng = m.get_unmixing_matrix(0) + + corrs = xcorr(W_fortran, W_ng) + n_above_95 = int((corrs > 0.95).sum()) + print( + f"seed {seed}: mean_corr={corrs.mean():.4f} min={corrs.min():.4f} " + f"n_above_0.95={n_above_95}/{nw} " + f"fortran_LL={fort_ll:.4f} ng_LL={m.final_ll_:.4f}", + flush=True, + ) + results.append(corrs.mean()) + + if results: + print( + f"\n{nw}ch/{field}fr (k={k:.1f}), n={len(results)}: " + f"mean={np.mean(results):.4f} range={np.min(results):.4f}-{np.max(results):.4f}", + flush=True, + ) + + if len(results) < n_seeds: + print( + f"\n{n_seeds - len(results)}/{n_seeds} seed(s) failed or were degenerate", + flush=True, + ) + return 1 + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/.context/issue-27/multimodel_ensemble.py b/.context/issue-27/multimodel_ensemble.py index 9027aae..dd5f8c1 100644 --- a/.context/issue-27/multimodel_ensemble.py +++ b/.context/issue-27/multimodel_ensemble.py @@ -7,6 +7,7 @@ next to this script. Real sample data + Fortran binary only (NO MOCK). uv run python .context/issue-27/multimodel_ensemble.py [N] + uv run python .context/issue-27/multimodel_ensemble.py --from-cache # reuse ensemble.npz, no re-fitting The Fortran binary (`sample_data/amica15mac`) is x86_64 and runs under Rosetta on Apple silicon. Absolute cross-corr magnitudes depend on config/seed; the @@ -161,10 +162,16 @@ def gap(mask): # within-group-A minus A-vs-rest (large => between worse) def figure(within_F, within_G, between, F_ll, G_ll, diff, p_perm, ks, out): + # Sized to its ACTUAL print footprint, not a big on-screen canvas: paper.md + # embeds this at width=100% of a ~5.36in single-column page (measured from the + # compiled paper.pdf), so figsize is set to that width directly -- LaTeX then + # displays it near 1:1 instead of shrinking a much larger canvas down to fit, + # which previously collapsed every font to an unreadable ~3pt in the printed + # PDF even though it looked fine on screen (figure-qa print-scale finding). plt.rcParams.update( - {"font.size": 11, "axes.spines.top": False, "axes.spines.right": False} + {"font.size": 7, "axes.spines.top": False, "axes.spines.right": False} ) - fig, (axA, axB) = plt.subplots(1, 2, figsize=(11, 4.3)) + fig, (axA, axB) = plt.subplots(1, 2, figsize=(5.36, 4.2)) bins = np.linspace(0.5, 1.0, 26) for arr, col, lab in [ (within_F, C_FORT, "within-Fortran"), @@ -176,21 +183,13 @@ def figure(within_F, within_G, between, F_ll, G_ll, diff, p_perm, ks, out): arr, bins=bins, density=True, histtype="step", color=col, lw=2, label=lab ) axA.axvline(arr.mean(), color=col, ls="--", lw=1.2) - axA.set_xlabel("stacked 2x32 Hungarian cross-correlation") + axA.set_xlabel("Hungarian-matched cross-correlation\n(stacked 2x32 components)") axA.set_ylabel("density") - axA.set_title("A Partition-agreement distributions", loc="left", fontweight="bold") - axA.legend(frameon=False, fontsize=9, loc="upper right") - axA.text( - 0.02, - 0.97, - f"means: within-F {within_F.mean():.3f} | within-pyA {within_G.mean():.3f} | " - f"between {between.mean():.3f}\n" - f"means differ by {abs(diff):.3f} (within a 0.05 margin)\n" - f"run-level permutation (between not worse): p={p_perm:.2f}", - transform=axA.transAxes, - va="top", - fontsize=8.5, - bbox=dict(boxstyle="round", fc="white", ec="0.7", alpha=0.9), + axA.set_title( + "A Partition-agreement distributions", + loc="left", + fontweight="bold", + fontsize=8, ) bins_ll = np.linspace( min(F_ll.min(), G_ll.min()) - 0.005, max(F_ll.max(), G_ll.max()) + 0.005, 24 @@ -201,56 +200,116 @@ def figure(within_F, within_G, between, F_ll, G_ll, diff, p_perm, ks, out): arr, bins=bins_ll, density=True, histtype="step", color=col, lw=2, label=lab ) axB.axvline(arr.mean(), color=col, ls="--", lw=1.2) - axB.set_xlabel("final log-likelihood (per sample-channel)") + axB.set_xlabel("final log-likelihood\n(mean per sample-channel)") axB.set_ylabel("density") - axB.set_title("B Likelihood distributions", loc="left", fontweight="bold") - axB.legend(frameon=False, fontsize=9, loc="upper left") - axB.text( - 0.98, - 0.97, - f"Fortran {F_ll.mean():.4f} (sd {F_ll.std():.3f})\n" - f"pyAMICA {G_ll.mean():.4f} (sd {G_ll.std():.3f})\nKS p={ks:.0e}", - transform=axB.transAxes, + axB.set_title( + "B Likelihood distributions", loc="left", fontweight="bold", fontsize=8 + ) + + # Both the legend and the stats box previously sat inside the axes and ended up + # overlapping the histogram bars (and each other) no matter where they were + # anchored -- with 3 distributions filling most of the plotted range, there was + # no empty pocket big enough to hold either. Reserve a fixed bottom margin for + # both instead and place them with figure-fraction (not axes-fraction) + # coordinates: axes-fraction anchors turned out to depend on the final axes + # height that tight_layout picks, which isn't known in advance and caused the + # legend/text to collide with the xlabel above it. get_position() gives each + # axes' true horizontal center after the layout below is fixed, so this keeps + # each panel's legend/text under its own histogram, not bleeding into the + # other panel. + fig.subplots_adjust(top=0.80, bottom=0.42, left=0.11, right=0.97, wspace=0.45) + cx_a = sum(axA.get_position().intervalx) / 2 + cx_b = sum(axB.get_position().intervalx) / 2 + + handles_a, labels_a = axA.get_legend_handles_labels() + fig.legend( + handles_a, labels_a, frameon=False, fontsize=6, + loc="upper center", bbox_to_anchor=(cx_a, 0.31), + ) # fmt: skip + fig.text( + cx_a, + 0.20, + f"mean corr. (run pairs)\n" + f"within-Fortran: {within_F.mean():.3f} (n={len(within_F)})\n" + f"within-pyAMICA: {within_G.mean():.3f} (n={len(within_G)})\n" + f"between: {between.mean():.3f} (n={len(between)})\n" + f"diff: {diff:+.3f} (margin +/-0.05)\n" + f"perm. p={p_perm:.2f}", + ha="center", va="top", - ha="right", - fontsize=8.5, + fontsize=5.5, + bbox=dict(boxstyle="round", fc="white", ec="0.7", alpha=0.95), + ) + + handles_b, labels_b = axB.get_legend_handles_labels() + fig.legend( + handles_b, labels_b, frameon=False, fontsize=6, + loc="upper center", bbox_to_anchor=(cx_b, 0.31), + ) # fmt: skip + fig.text( + cx_b, + 0.23, + f"mean final LL\n" + f"Fortran: {F_ll.mean():.4f} (sd {F_ll.std():.3f})\n" + f"pyAMICA: {G_ll.mean():.4f} (sd {G_ll.std():.3f})\n" + f"gap: {abs(F_ll.mean() - G_ll.mean()):.3f} (100-iter budget)\n" + f"KS p={ks:.0e}", + ha="center", + va="top", + fontsize=5.5, bbox=dict(boxstyle="round", fc="white", ec="0.7", alpha=0.9), ) + + fig.text( + 0.5, + 0.90, + "Dashed vertical lines mark each distribution's mean.", + ha="center", + fontsize=6.5, + style="italic", + ) fig.suptitle( "Multi-model AMICA (n_models=2): pyAMICA vs Fortran ensembles, real sample EEG", fontweight="bold", - y=1.02, + fontsize=8.5, + y=0.97, ) - fig.tight_layout() fig.savefig( - out / "multimodel_ensemble_distributions.png", bbox_inches="tight", dpi=200 + out / "multimodel_ensemble_distributions.png", bbox_inches="tight", dpi=300 ) fig.savefig(out / "multimodel_ensemble_distributions.pdf", bbox_inches="tight") def main(): - n = int(sys.argv[1]) if len(sys.argv) > 1 else 20 - data = load_data() - work = Path(tempfile.mkdtemp(prefix="amica_ensemble_")) - print(f"scratch: {work}") - Fs, F_ll, Gs, G_ll = [], [], [], [] - for k in range(n): - print(f"Fortran {k + 1}/{n}", flush=True) - W, ll = run_fortran(work, str(k)) - Fs.append(W) - F_ll.append(ll) - for k in range(n): - print(f"NG {k + 1}/{n}", flush=True) - W, ll = run_ng(data, seed=k) - Gs.append(W) - G_ll.append(ll) - Fs, Gs = np.array(Fs), np.array(Gs) - F_ll, G_ll = np.array(F_ll), np.array(G_ll) - - # Persist the raw ensemble so the figure/tests can be regenerated without - # re-running the 40 fits (the earlier run's data was lost, forcing this rerun). - np.savez(HERE / "ensemble.npz", Fs=Fs, Gs=Gs, F_ll=F_ll, G_ll=G_ll) + if "--from-cache" in sys.argv: + # Reuse the persisted ensemble (e.g. to regenerate the figure after a + # labeling/legend change) instead of re-running 40 real fits. + d = np.load(HERE / "ensemble.npz") + Fs, Gs, F_ll, G_ll = d["Fs"], d["Gs"], d["F_ll"], d["G_ll"] + else: + n = int(sys.argv[1]) if len(sys.argv) > 1 else 20 + data = load_data() + work = Path(tempfile.mkdtemp(prefix="amica_ensemble_")) + print(f"scratch: {work}") + Fs, F_ll, Gs, G_ll = [], [], [], [] + for k in range(n): + print(f"Fortran {k + 1}/{n}", flush=True) + W, ll = run_fortran(work, str(k)) + Fs.append(W) + F_ll.append(ll) + for k in range(n): + print(f"NG {k + 1}/{n}", flush=True) + W, ll = run_ng(data, seed=k) + Gs.append(W) + G_ll.append(ll) + Fs, Gs = np.array(Fs), np.array(Gs) + F_ll, G_ll = np.array(F_ll), np.array(G_ll) + + # Persist the raw ensemble so the figure/tests can be regenerated without + # re-running the 40 fits (the earlier run's data was lost, forcing this rerun). + np.savez(HERE / "ensemble.npz", Fs=Fs, Gs=Gs, F_ll=F_ll, G_ll=G_ll) + n = len(Fs) within_F = pairwise(Fs, Fs, True) within_G = pairwise(Gs, Gs, True) between = pairwise(Gs, Fs, False) diff --git a/.context/issue-27/multimodel_ensemble_distributions.pdf b/.context/issue-27/multimodel_ensemble_distributions.pdf index 9a2de52..9088e1d 100644 Binary files a/.context/issue-27/multimodel_ensemble_distributions.pdf and b/.context/issue-27/multimodel_ensemble_distributions.pdf differ diff --git a/.context/issue-27/multimodel_ensemble_distributions.png b/.context/issue-27/multimodel_ensemble_distributions.png index ead9146..98f41ca 100644 Binary files a/.context/issue-27/multimodel_ensemble_distributions.png and b/.context/issue-27/multimodel_ensemble_distributions.png differ diff --git a/.context/paper-applications-draft.md b/.context/paper-applications-draft.md new file mode 100644 index 0000000..b6d6eb3 --- /dev/null +++ b/.context/paper-applications-draft.md @@ -0,0 +1,60 @@ +# Draft: "Applications" section (for a fuller paper, not the JOSS 1750-word draft) + +Contributed by Scott Makeig (2026-07-14), suggested to go before or after "State of the Field" +in `paper.md`. Held here rather than merged into the JOSS paper: the JOSS draft is word-budget- +constrained and this section is scoped for a longer paper/preprint. Revisit when drafting that. + +--- + +## Applications + +The ever-advancing compute speed of desktop computing can facilitate routine applications of +advanced brain data modeling features available in AMICA: + +1. **Source separation.** We recently have shown (Gwen et al., 202?) that the default max + iterations (2000) proposed in AMICA is a convenience value -- AMICA source separation continues + to increase slowly as the number of training iterations is increased. At 25 ms per iteration, a + 2000-iteration decomposition (e.g., of the 70-channel example EEG datasets used in this paper) + should require less than an hour to compute, in many cases making available a larger compute + horizon within which to further optimize source separation. + + An efficient measure of source separation performance is Mutual Information Reduction (MIR) + introduced by Palmer in (Delorme et al., 2012). + + Scott's strong recommendation: add MIR as a built-in pyAMICA option, applicable at decomposition + end and/or optionally at specified decomposition waypoints. See the MIR/PMI port epic (tracked + as a GitHub issue) for the implementation side of this. + +2. **Brain dynamic instability** is a hallmark of human brain dynamics, both normal and + pathologic -- yet source and source network instability is not yet commonly measured in M/EEG + studies. In its multi-model mode, AMICA separates its training data into domains fit to + different ICA models that compete for data points during training; this separation of the + training data into source-model domains has been shown to be powerful for brain state + monitoring, during sleep, quiet rest, and active task performance (Hsu et al., 201?). Accurate, + data-driven segregation of unlabeled datasets into as many as 20 rest and active emotion + imagination periods has also been demonstrated (Hsu et al., 201?). A plug-in EEGLAB toolbox for + evaluating and plotting multi-model AMICA solutions is available (Ozgur...). + +3. **Source and source network stability.** Artoni et al. (191?) have demonstrated and + contributed RELICA, an EEGLAB plug-in for estimating the stability of sources returned by ICA. + Akalin Acar & Makeig (202?) showed that the pattern of scalp projection variance of clusters of + near-identical sources returned across bootstrap training data decompositions can reveal the + nature of biological source network instability. Again, at 25 ms per training iteration, + twenty-five 2000-iteration bootstrap decompositions of a given dataset can be computed in less + than a day, allowing assays of source and source network instability as well as stability. + +All these possibilities, once beyond the reach of routine desktop computing, are now practical to +apply using current desktop hardware. Using pyAMICA in still more powerful compute environments can +only increase the depth of detail and statistical power of studies of brain dynamics in complex or +even real-life protocols -- either exploratory or confirmatory. For applications in supercomputer +environments, however, Palmer's FORTRAN version (AMICA 5.??) customized for use at the San Diego +Supercomputer Center (freely available via the Neuroscience Gateway (nsgportal.org; ????, 201?)) +might prove still more efficient. + +--- + +Notes for whoever drafts the fuller paper: +- All the "???"/"201?" citations above are as Scott sent them; need real citations before use. +- The 25 ms/iteration figure matches the MLX benchmark (`.context/issue-77/benchmark_findings.md`), + not every backend -- qualify per-backend if reused. +- Ties into [[amica-parity-epic-status]] and the postAmicaUtility MIR/PMI port epic. diff --git a/.github/workflows/draft-pdf.yml b/.github/workflows/draft-pdf.yml index 0d2813c..bcc232f 100644 --- a/.github/workflows/draft-pdf.yml +++ b/.github/workflows/draft-pdf.yml @@ -1,8 +1,11 @@ name: Draft PDF # Builds paper.md into a JOSS-formatted PDF (Pandoc + inara, via the official -# openjournals action) and uploads it as an artifact. This verifies the paper -# compiles under JOSS's toolchain before submission and gives a per-push preview. +# openjournals action). Not a JOSS submission step -- this just keeps the +# committed paper.pdf in the repo root in sync with paper.md/paper.bib for +# local review. On a direct push it commits the rebuilt PDF back (skipped for +# pull_request runs, since a PR checkout is a detached merge ref, not the head +# branch, so pushing there would go to the wrong place or fail outright). on: push: paths: @@ -16,6 +19,9 @@ on: - .github/workflows/draft-pdf.yml workflow_dispatch: +permissions: + contents: write + jobs: paper: runs-on: ubuntu-latest @@ -34,3 +40,15 @@ jobs: name: paper # The action writes the compiled PDF next to the paper source. path: paper.pdf + - name: Commit and push rebuilt PDF + if: github.event_name == 'push' + run: | + git config user.name "github-actions[bot]" + git config user.email "github-actions[bot]@users.noreply.github.com" + git add paper.pdf + if ! git diff --cached --quiet; then + git commit -m "Rebuild paper.pdf [skip ci]" + git push + else + echo "paper.pdf unchanged, nothing to commit" + fi diff --git a/.gitignore b/.gitignore index af17e45..c436b66 100644 --- a/.gitignore +++ b/.gitignore @@ -26,6 +26,7 @@ htmlcov/ # Temp output from NG e2e test (Fortran run artifacts) pyAMICA/tests/torch_tests/_ng_e2e_tmp/ +pyAMICA/tests/torch_tests/_ng_e2e_json_tmp/ # Benchmark run artifacts benchmarks/results/ diff --git a/README.md b/README.md index cf21701..f5d6a55 100644 --- a/README.md +++ b/README.md @@ -12,8 +12,8 @@ EEG/EMG blind source separation and is a drop-in replacement for EEGLAB's AMICA: single-model output is byte-identical to the Fortran reference and loads directly in EEGLAB. -Single-model results match the Fortran reference (log-likelihood ~ -3.40, -Hungarian-matched component correlation ~ 0.997); see the +Single-model results match the Fortran reference (Hungarian-matched component correlation ~ 0.998 +on well-determined data, Newton disabled); see the [documentation](https://eeglab.org/pyAMICA/) for validation details and the backend-selection guide. diff --git a/docs/assets/figures/multimodel-ensemble.png b/docs/assets/figures/multimodel-ensemble.png index ead9146..98f41ca 100644 Binary files a/docs/assets/figures/multimodel-ensemble.png and b/docs/assets/figures/multimodel-ensemble.png differ diff --git a/docs/guides/validation.md b/docs/guides/validation.md index 434cfeb..90a2ee6 100644 --- a/docs/guides/validation.md +++ b/docs/guides/validation.md @@ -3,9 +3,15 @@ **Correctness in pyAMICA is defined as parity with the reference Fortran binary, not merely as convergence.** A run is correct when it reproduces the Fortran output within numerical tolerance. This page collects the full verification evidence: bit-exact score functions, single-model parity, -multi-model distributional equivalence, cross-platform device and precision invariance, +multi-model distributional similarity, cross-platform device and precision invariance, the EEGLAB drop-in round-trip, and the remaining validated behaviors. -Every result uses the bundled real sample EEG and the reference Fortran binary; none uses synthetic data. +Every result uses real EEG and the reference Fortran binary; none uses synthetic data. +The multi-model and score-function checks use only the bundled EEGLAB tutorial sample (32 channels, +30504 samples at 128 Hz, ~238 s); the single-model headline additionally uses an external recording +(OpenNeuro ds002718 sub-002, 70 channels) at $k=\text{frames}/\text{channels}^2\approx153$, well past +the ~60 threshold where cross-backend agreement plateaus (below), since the bundled sample sits at +that threshold's boundary ($k\approx30$). Extending to a multi-subject, multi-dataset validation is +planned future work, not yet done. Throughout, IC abbreviates independent component and LL log-likelihood. ## Validation at a glance @@ -14,8 +20,9 @@ Throughout, IC abbreviates independent component and LL log-likelihood. |---|---|---| | Source-density score and log-density (non-GG families) | vs the literal `amica15.f90` expressions | bit-exact ($<10^{-12}$) | | Per-block sufficient statistics and one M-step | vs Fortran | bit-exact ($\sim\!10^{-15}$) | -| Single-model solution | log-likelihood, component correlation, Amari distance vs Fortran | LL within ~0.005 of $-3.4018$; correlation 0.997; Amari 0.006 | -| Multi-model solution | distributional equivalence over 20-run ensembles | indistinguishable from Fortran's own run-to-run spread ($p = 0.96$) | +| Single-model solution (`do_newton=0`, $k\approx153$) | log-likelihood, component correlation vs Fortran | LL within ~0.0005 of $-3.6993$; correlation 0.998 | +| Single-model solution (`do_newton=0`, bundled, $k\approx30$) | Amari distance vs Fortran | 0.006 | +| Multi-model solution | distributional similarity over 20-run ensembles | indistinguishable from Fortran's own run-to-run spread ($p = 0.96$) | | Device and precision invariance | same independent components across CPU/CUDA/MPS/MLX, float32/float64, Linux/macOS | identical (1.000) across all eight torch/MLX combinations | | Cross-backend log-likelihood | converged LL across every backend | agree to ~3 significant digits (max pairwise ~0.003) | | EEGLAB output | `write_amica_output` round-trip through `loadmodout15` | single-model bytes are an exact serialization; loads with correct layout | @@ -32,14 +39,21 @@ a standard unmixing-matrix comparison metric (Amari, Cichocki & Yang, 1996) that ## Single-model parity -On real sample EEG the natural-gradient backend reaches Fortran's solution: +With Newton disabled (`do_newton=0`), the natural-gradient backend reaches Fortran's solution, +averaged over 5 seeds at a matched 2000-iteration budget. The headline correlation is measured on an +external, unambiguously well-determined recording (OpenNeuro ds002718, 70 channels, $k\approx153$) so +it is not sensitive to whether a particular small dataset happens to be well-conditioned; the bundled +32-channel sample ($k\approx30$, at the project's own data-adequacy boundary) gives a consistent Amari +distance: -- Log-likelihood ~ -3.40 (Fortran ~ -3.4018). -- Hungarian-matched component correlation ~0.997, clearing the >0.95 gate. -- Amari distance ~0.006. +- Log-likelihood ~ -3.6993 ($k\approx153$; Fortran ~ -3.6993, gap ~0.0003). +- Hungarian-matched component correlation ~0.998 ($k\approx153$; Fortran-vs-Fortran self-consistency + over the same 5 seeds: ~0.999), clearing the >0.95 gate. On the bundled sample ($k\approx30$) both + numbers are consistent: ~0.998 pyAMICA-vs-Fortran, ~0.998 Fortran-vs-Fortran. +- Amari distance ~0.006 (bundled sample; Fortran-vs-Fortran: ~0.005). The fixed source-density families are bit-exact against the literal Fortran score/derivative expressions (~1e-15), -and the backend converges to the binary's solution within ~0.005 log-likelihood. +and the backend converges to the binary's solution within ~0.005 log-likelihood on either dataset. ### Source-density families are bit-exact @@ -66,10 +80,10 @@ sub-Gaussian (code 4) densities on a kurtosis schedule; its dynamic switch has n log-likelihood instead. Each fixed family converges within ~0.005 LL of the binary at a matched Newton budget. See `pyAMICA/tests/torch_tests/test_ng_pdf_families.py` and ADR 0002. -## Multi-model equivalence +## Multi-model distributional similarity Multi-model AMICA is not partition-identifiable, so exact partition parity with Fortran is the wrong acceptance bar. -The right test is whether the two implementations sample the same distribution over solutions. Running an ensemble of `N = 20` fits per implementation on the real sample EEG (`n_models = 2`, 3 mixture components, 100 iterations, matched schedule), +The right test is whether the two implementations sample a similar distribution over solutions. Running an ensemble of `N = 20` fits per implementation on the bundled sample EEG (`n_models = 2`, 3 mixture components, 100 iterations, matched schedule), the pyAMICA-vs-Fortran partition cross-correlation distribution overlaps Fortran's own run-to-run distribution: | Distribution (pairwise Hungarian-matched \|corr\|) | Mean | SD | Range | @@ -374,14 +388,35 @@ Tests live under `pyAMICA/tests/`: `torch_tests/test_ng_backend.py`, `torch_test ## Reproducing these results -Everything on this page runs from the bundled sample data and the reference binary, with no external download: +The multi-model and score-function checks run from the bundled sample data and the reference binary, +with no external download: ```bash uv run python validate_implementations.py # single- and multi-model parity report uv run pytest # the full parity/behavior test suite ``` +The bundled-sample single-model numbers (`do_newton=0`, ~0.998 correlation / 0.006 Amari, $k\approx30$) +are the 5-seed, 2000-iteration sweep produced by +`uv run python .context/issue-144-parity-data-adequacy/bundled_sample_newton0.py 5 2000`. +The headline correlation/LL ($k\approx153$) additionally needs the external ds002718 recording +(OpenNeuro, sub-002, manual download) and is produced by +`.context/issue-144-parity-data-adequacy/run_5seed_newton0.sh`; that script's paths, thread count, and +Fortran binary are specific to the workstation it was run on, so it is not yet a portable one-command +reproduction (tracked as future work). `validate_implementations.py`'s own CLI defaults (a single seed, +100 iterations, `do_newton` read from `sample_params.json` where it is `true`) do not reproduce either +row directly. + The multi-model ensemble and Amari detail regenerate from saved fits (no re-fitting) with `uv run python .context/issue-27/amari_distance.py`. The cross-platform benchmark and equivalence figures are produced by `benchmarks/benchmark_decompose.py` (and the sweep scripts alongside it); the underlying findings are in `.context/issue-84/` and `.context/issue-90/`. + +### Parameter files + +`sample_data/sample_params.json` is the JSON parameter file used above (loaded via +`AMICA.from_params_file`); its keys mostly reuse Fortran's `.param` names (`lrate`, `do_newton`, +`rho0`, `block_size`, `max_iter`, `num_models`, ...), but not all of them match one-to-one +(for example `num_mix` here vs `num_mix_comps` in Fortran's `input.param`), and pyAMICA does not +yet parse the literal Fortran `.param` text format. A native `.param` reader, so the same file +drives both implementations, is tracked as a future issue. diff --git a/paper.md b/paper.md index 65c1fb5..af650a2 100644 --- a/paper.md +++ b/paper.md @@ -32,34 +32,35 @@ bibliography: paper.bib Independent Component Analysis (ICA) is a standard method for separating electroencephalography (EEG) and electromyography (EMG) recordings into maximally independent sources, which isolates brain, muscle, and artifact activity for downstream analysis. Adaptive Mixture ICA (AMICA) generalizes single-model ICA to a mixture of ICA models with adaptive source densities, and produces the most dipolar (and thus most physiologically interpretable) component decompositions of EEG among the widely used algorithms benchmarked by @delorme2012independent. -Its reference implementation, however, is a Fortran program parallelized with the Message Passing Interface (MPI) and distributed as a compiled binary driven from MATLAB/EEGLAB, +Its reference implementation is a Fortran program parallelized with the Message Passing Interface (MPI) and distributed as a compiled binary driven from MATLAB/EEGLAB, which is difficult to install, runs only on the CPU, and is not usable from a Python scientific workflow. `pyAMICA` is a Python implementation of AMICA that reproduces the reference Fortran results within numerical tolerance while running on the CPU, NVIDIA GPUs (CUDA), and Apple GPUs (Apple's MLX array framework [@mlx2023]). -It is built on PyTorch [@paszke2019pytorch], NumPy [@harris2020array], -and SciPy [@virtanen2020scipy], and exposes a scikit-learn-style estimator. +It is a complete NumPy/PyTorch reimplementation, not a wrapper around the Fortran binary, +built on PyTorch [@paszke2019pytorch], NumPy [@harris2020array], +and SciPy [@virtanen2020scipy] (Python $\geq$3.12, PyTorch $\geq$2.12; tested on Apple Silicon/ARM64 and `x86_64`/CUDA Linux), and exposes a scikit-learn-style estimator under a BSD-3-Clause license. In double precision it reproduces the reference score-function algebra to machine precision; it also runs in single precision (float32), which the CPU-only binary does not offer and which is numerically faithful (agreeing with double precision to four to five significant digits) yet required to use Apple GPUs, which have no float64 and host the fastest backend (MLX). `pyAMICA` writes output in the binary format that EEGLAB's AMICA loader reads: a single-model output file is byte-identical in layout to a native AMICA file and needs no manual re-interpretation, -and multi-model output round-trips through the same loader. Correctness is defined as parity with the Fortran reference for the single-model case and, because multi-model AMICA is not partition-identifiable, as distributional equivalence for the multi-model case; +and multi-model output round-trips through the same loader. Correctness is defined as parity with the Fortran reference for the single-model case and, because multi-model AMICA is not partition-identifiable, as a similar distribution of solutions for the multi-model case; both are validated on real EEG against the reference binary. +A minimal example is a scikit-learn-style call, `AMICA(n_models=1, n_mix=3).fit(X)` on a `(channels, samples)` array (full workflow in the README). The software is at (archived at doi:10.5281/zenodo.21312148). # Statement of need AMICA yields components that are well suited to equivalent-dipole source localization and to automated classification of independent components [@piontonachini2019iclabel], and it separates EEG more effectively than most alternatives [@delorme2012independent; @palmer2012amica]. -The reference implementation is the original Fortran code: it depends on an MPI toolchain and precompiled binaries that are awkward to build across platforms, -it cannot use a GPU, and it is invoked as an external process from MATLAB rather than called from Python. As neuroimaging analysis has moved toward Python, +The reference Fortran implementation shares these limitations and is invoked as an external process from MATLAB rather than called from Python. As neuroimaging analysis has moved toward Python, for example MNE-Python [@gramfort2013meg], an AMICA that runs natively in Python and on a GPU, and that is validated to reproduce the Fortran reference numerically, is needed for modern pipelines and for studying the algorithm in an open codebase. General-purpose Python ICA implementations do not fill this gap. `scikit-learn` and `MNE-Python` provide FastICA [@hyvarinen2000independent] and Infomax [@bell1995information; @lee1999independent], and Picard [@ablin2018faster] provides a faster maximum-likelihood ICA, but none implement AMICA's mixture of models, adaptive generalized-Gaussian source densities, or Newton updates, -and so they do not reproduce AMICA decompositions. `pyAMICA` targets researchers who need AMICA specifically: +and so they do not reproduce AMICA decompositions. `pyAMICA` targets researchers who need AMICA: EEG/EMG analysts who want AMICA-quality decompositions inside a Python pipeline, users of GPU hardware who want faster runs than the CPU-only binary, and methodologists who need a transparent reference implementation to inspect and build on. @@ -70,42 +71,48 @@ exact-EM mixture updates, a positive-definite Newton step [@palmer2008newton], symmetric zero-phase-component-analysis (ZCA) sphering, the five source-density families of the reference (generalized Gaussian, Gaussian, logistic, sub-Gaussian, and the extended-Infomax kurtosis switcher), a mixture of ICA models, and component sharing across models. -`pyAMICA`'s conformity with the reference binary is measured on real sample EEG (Table 1) with two complementary metrics: -Hungarian-matched component correlation, and the Amari distance [@amari1996new], a standard unmixing-matrix comparison metric that needs no assignment step since it is permutation- and scale-invariant by construction. -Both agree closely for the single model, and the source-density score functions and per-block sufficient statistics are bit-exact against the literal Fortran expressions. +`pyAMICA`'s conformity with the reference binary is measured with two complementary metrics: Hungarian-matched component correlation, +and the Amari distance [@amari1996new], a standard unmixing-matrix comparison metric that needs no assignment step since it is permutation- and scale-invariant by construction. +Both implementations were run for AMICA's usual 2000 iterations with Newton off (`do_newton=0`) and otherwise-default parameters +(two separate configuration files drive them, JSON for `pyAMICA` and Fortran's own text format, with a transcribed subset of settings; a shared-format reader is planned). +The single-model headline (Table 1) uses an external recording (OpenNeuro ds002718, 70 channels, $k\approx153$), well past the ~60 threshold where cross-backend agreement plateaus (documentation); +the bundled 32-channel sample used below sits at that threshold's boundary ($k\approx30$) and gives a consistent Amari distance. +Score functions and per-block sufficient statistics are exact to floating-point resolution against the literal Fortran expressions on the bundled sample. A mixture of ICA models is not partition-identifiable, so exact partition parity is the wrong bar for the multi-model case; -it is instead assessed by distributional equivalence across ensembles of 20 runs each (\autoref{fig:ensemble}). -Both metrics again agree: a run-level permutation test, which permutes the 40 runs as intact units to respect the dependence among pairwise values, -finds no evidence that cross-implementation agreement is worse than Fortran's own run-to-run agreement. -The single-run values are therefore intrinsic estimator spread rather than a shortfall. -Equivalence is claimed for the partition structure; the multi-model log-likelihood distributions still differ slightly at a matched 100-iteration budget -(`pyAMICA` reaches Fortran's mean with about twice as many iterations), so full-likelihood equivalence is not yet claimed. Per-run detail for both metrics is in the documentation. - -| Regime | Metric | Result (correlation / Amari distance) | -|---|---|---| -| Single-model | Log-likelihood gap to Fortran ($-3.4018$) | within ~0.005 | -| Single-model | Conformity with Fortran | 0.997 / 0.006 | -| Single-model | Score functions and sufficient statistics | bit-exact ($\sim\!10^{-15}$) | -| Multi-model | Single run (`pyAMICA`-Fortran; Fortran-Fortran) | 0.65; 0.64 (sd 0.05) / 0.163; 0.174 (sd 0.02) | -| Multi-model | Ensemble equivalence, 20 runs each (permutation $p$) | $p=0.96$ / $p>0.999$ | - -: Single-model parity and multi-model distributional equivalence of `pyAMICA` -with the Fortran reference on the bundled sample EEG. +it is instead assessed by whether the two implementations sample a similar distribution of solutions, across ensembles of 20 bundled-sample runs each (\autoref{fig:ensemble}). +A run-level permutation test, which permutes the 40 runs as intact units to respect the dependence among pairwise values, finds no evidence that cross-implementation +agreement is worse than Fortran's own run-to-run agreement, so single-run values reflect intrinsic estimator spread rather than a shortfall. +The multi-model log-likelihood distributions still differ slightly at a matched 100-iteration budget +(`pyAMICA` reaches Fortran's mean with about twice as many iterations), so full-likelihood similarity is not yet claimed. Per-run detail for both metrics is in the documentation. -![Multi-model solution-ensemble partition-correlation distributions for 20 `pyAMICA` and 20 Fortran fits of the sample EEG. -The within-Fortran, within-`pyAMICA`, and between-implementation distributions overlap, -so the single-run correlation reflects the estimator's intrinsic run-to-run spread rather than a gap to the reference.\label{fig:ensemble}](docs/assets/figures/multimodel-ensemble.png){ width=75% } +| Regime | Metric | Result (mean; correlation / Amari distance) | +|---|---|---| +| Single-model | Log-likelihood gap to Fortran (mean per-sample-channel LL, $-3.6993$, $k\approx153$) | within ~0.0005 | +| Single-model | Conformity with Fortran | 0.998 ($k\approx153$) / 0.006 (bundled, $k\approx30$) | +| Single-model | Score functions (non-default families) and sufficient statistics | exact to floating-point resolution ($\sim\!10^{-15}$) | +| Multi-model | Single-run magnitude (`pyAMICA`-Fortran; Fortran-Fortran) | 0.65; 0.64 (sd 0.05) / 0.163; 0.174 (sd 0.02) | +| Multi-model | Ensemble similarity, 20 runs each: mean difference, between $-$ within-Fortran (permutation $p$) | $+0.011$ ($p=0.96$) / $-0.011$ ($p>0.999$) | + +: Single-model parity (external ds002718, $k\approx153$; Amari and score-function checks on the bundled sample) and +multi-model distributional similarity (bundled sample) of `pyAMICA` with the Fortran reference. All values are means +(sd shown where relevant) over matched components (single-model) or 190 within- and 400 cross-implementation pairs (multi-model). +Full methodology and reproduction steps are in the documentation. + +![Multi-model solution-ensemble partition-correlation distributions (panel A) and log-likelihood distributions (panel B) for 20 `pyAMICA` and 20 Fortran fits of the sample EEG; dashed lines mark each distribution's mean. +The within-Fortran, within-`pyAMICA`, and between-implementation correlation distributions overlap, +so the single-run correlation reflects the estimator's intrinsic run-to-run spread rather than a gap to the reference. +Panel B's apparent separation is a 0.009 log-likelihood gap on a ~0.035 axis.\label{fig:ensemble}](docs/assets/figures/multimodel-ensemble.png){ width=100% } All backends converge to the same single-model log-likelihood on real EEG (maximum pairwise difference ~0.003). -On Apple Silicon the MLX backend is the fastest option and is roughly flat with channel count (15-25 ms per iteration from 16 to 70 channels; see the documentation), +On Apple Silicon the MLX backend is the fastest option and is flat with channel count (15-25 ms per iteration from 16 to 70 channels; see the documentation), about eight times faster than double-precision multithreaded CPU; PyTorch-MPS is not a win (at or worse than the CPU). On NVIDIA hardware double-precision CUDA is the reproducible path and is overhead-bound at EEG scale, so single precision gives it little additional speedup (Table 2). Native Fortran itself scales with CPU cores, unlike the CPU backends above: -with enough cores pinned it is competitive with, or faster than, the GPU it is compared against on each machine (Table 2), -though only by dedicating most cores of a much larger, hotter host than a laptop GPU. -A data-size sweep further shows cross-backend component equivalence rising with frames per channel and plateauing near 0.98 once the decomposition is well-determined, -where two independent double-precision implementations (native Fortran and PyTorch-CUDA) agree at 0.995; -single-precision runs are seven-significant-digit, so double precision remains the default for parity. +with enough cores pinned it beats the CUDA GPU on the workstation (Table 2), though only by dedicating most cores of a larger, hotter host than a laptop GPU; +it does not catch Apple's MLX on laptop hardware. +A data-size sweep further shows cross-backend component agreement rising with frames per channel and plateauing near 0.98 once the decomposition is well-determined, +where two independent double-precision implementations (native Fortran and PyTorch-CUDA) agree at a mean of 0.995; +single-precision runs agree with double precision to four to five significant digits, so double precision remains the default for parity. | Backend (device) | Precision | ms / iteration | |---|---|---:| @@ -121,15 +128,15 @@ single-precision runs are seven-significant-digit, so double precision remains t `pdftype`=0, `block_size`=512; warm, minimum of repeated runs). CPU, MPS, and MLX on Apple Silicon; CUDA on a separate NVIDIA RTX 4090; CUDA float32 is comparable (~36 ms). -The two native-Fortran rows are from a separate core-count sweep (documentation) on the same two machines, at the core count where each backend's throughput levels off; +The two native-Fortran rows are from a separate core-count sweep (documentation) on the same two machines, at each backend's throughput plateau; the other CPU rows above use the platform default thread count, so they are not core-matched to Fortran. -Unlike the correctness comparison, this benchmark uses external data (OpenNeuro ds002718) and specific GPU hardware. +Unlike the correctness comparison, this benchmark uses external data (OpenNeuro ds002718, one subject so far) and specific GPU hardware. -The correctness harness compares `pyAMICA` against Fortran with two metrics, -Hungarian-matched component correlation and Amari distance, -and uses only the bundled real sample EEG and Fortran binary, with no external download and no synthetic data. +The correctness harness compares `pyAMICA` against Fortran with two metrics, Hungarian-matched component correlation and Amari distance, and never uses synthetic data; +the multi-model and score-function checks need no external download (bundled sample only). The full per-channel and multi-model performance tables, the per-run Amari-distance detail, -and the data-size sweep are in the documentation (). +and the data-size sweep, along with the step-by-step commands to reproduce every number here, +are in the documentation (). # State of the field @@ -137,8 +144,8 @@ and the data-size sweep are in the documentation ( 0.95 + assert amari_distance(fortran["W"], ng_results["W"]) < 0.05 + + @pytest.mark.slow @pytest.mark.skipif(not DATA_FILE.exists(), reason="sample data missing") def test_full_fit_parity_numpy_vs_ng(tmp_path): diff --git a/pyproject.toml b/pyproject.toml index beb6715..16e752f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -104,3 +104,15 @@ exclude_lines = [ # analysis scripts that ARE tested (loaded via importlib from the test suite) # were fixed on their own merits rather than excluded. (issue #122) exclude = [".context/**"] + +[tool.ty.rules] +# Optional-extra lazy imports (mne: benchmarks/benchmark_decompose.py; mlx: +# pyAMICA/mlx_impl/core.py, benchmarks/benchmark_dimsweep.py) carry +# `# ty: ignore[unresolved-import]` comments, since CI's base env (no extras, +# matching the other jobs) can't resolve either. Whether such a comment is +# "used" or "unused" then depends on whether the local env happens to have the +# extra installed, so this rule would otherwise flip between a false pass and +# a false "unused-ignore-comment" warning depending on the developer's own +# environment rather than the code. Disabled project-wide rather than per-line +# since there's no per-environment suppression for the checker itself. +unused-ignore-comment = "ignore" diff --git a/validate_implementations.py b/validate_implementations.py index e4f6b53..cfc8c5b 100644 --- a/validate_implementations.py +++ b/validate_implementations.py @@ -80,28 +80,6 @@ def load_sample_data() -> Tuple[np.ndarray, Dict]: return data, params -def create_fortran_params(params: Dict, output_dir: Path, seed: int) -> Path: - """Create parameter file for Fortran AMICA.""" - param_file = output_dir / "params.txt" - - # Create Fortran-style parameter file - with open(param_file, "w") as f: - f.write(f"datafile '{output_dir}/data.fdt'\n") - f.write(f"num_chans {params['data_dim']}\n") - f.write(f"num_frames {params['field_dim'][0]}\n") - f.write(f"num_models {params.get('num_models', 1)}\n") - f.write(f"num_mix {params.get('num_mix', 3)}\n") - f.write(f"max_iter {params.get('max_iter', 100)}\n") - f.write(f"lrate {params.get('lrate', 0.05)}\n") - f.write(f"do_mean {1 if params.get('do_mean', True) else 0}\n") - f.write(f"do_sphere {1 if params.get('do_sphere', True) else 0}\n") - f.write(f"do_newton {1 if params.get('do_newton', False) else 0}\n") - f.write(f"outdir '{output_dir}/fortran_output'\n") - f.write(f"seed {seed}\n") - - return param_file - - def run_fortran_amica( data: np.ndarray, params: Dict, output_dir: Path, seed: int ) -> Optional[Dict]: