Data-size (k-factor) sweep: cross-backend IC equivalence vs frames (#90)

.context/issue-84/phase3_decomp_equivalence.md

@@ -47,8 +47,10 @@ When the decomposition is well-determined (high k) **every backend, including Fo
 identical ICs** (0.997 at k=574). At k=30 -- the rule-of-thumb minimum -- only the strongest
 ~20/70 components are reproducible; the rest are under-determined and different inits settle
 into different, equally valid local optima (AMICA is non-convex). This validates both backend
-correctness *and* the data-adequacy rule. Pushing 70ch to k=152 with the full recording is the
-follow-up test in #90.
+correctness *and* the data-adequacy rule. Pushing 70ch to k=152 with the full recording (#90,
+`.context/issue-90/`) confirms the trend saturates: equivalence climbs to ~0.98 at k~=60 and
+plateaus (the two f64 implementations reach 0.995), it does not reach 1.0 because of float32
+precision + non-convex local-optimum spread on the weakest components.
 
 ## Result 3: total decompose time (2000 iters, 147k frames, 70ch)
 
@@ -85,5 +87,6 @@ MPS, on Apple.
   montage). Tracked informally; refine later.
 - Channel subsets use the first N electrodes (spatially clustered) -> #91 (use distributed
   subsets so reduced-channel maps are whole-head).
-- Data-size/k sweep at 70ch (full recording, k=152) -> #90.
+- Data-size/k sweep at 70ch (full recording, k=152) -> #90 DONE
+  (`.context/issue-90/ksweep_findings.md`): equivalence saturates at ~0.98 past k~=60.
 - Output-format/convention parity for a true EEGLAB drop-in -> #92.

.context/issue-90/ksweep_findings.md

@@ -0,0 +1,97 @@
+# Issue #90: data-size (k-factor) sweep at 70 channels
+
+Follow-up to Phase 3 (#87). Phase 3 swept k along the **channel** axis (frames fixed at 147k,
+so more channels = lower k) and found cross-backend IC equivalence tracks the data-adequacy
+factor **k = frames / channels^2** (EEGLAB rule of thumb: k >= 20-30). #90 sweeps the other
+axis: hold **channels = 70** and grow **frames** using the full ds002718 sub-002 recording
+(747,750 frames, ~50 min @ 250 Hz), so k rises from 15 to 152.
+
+Hypothesis: as k grows the decomposition becomes well-determined, so cross-backend equivalence
+should rise toward ~1.0.
+
+## Run
+`benchmarks/benchmark_decompose.py --frames 73500,147000,294000,490000,747750 --channels 70
+--iters 2000` on hallu (Linux x86_64, RTX 4090, 32 cores). Backends: `native-fortran-f64`,
+`torch-cuda-f64`, `torch-cuda-f32` (the fast, well-behaved backends -- Phase 3 already proved
+every torch/MLX device/precision recovers an identical decomposition, so a 3-backend cross-check
+here is sufficient and the large-frame runs are ~5x longer per backend than the k=30 runs).
+
+## Result: equivalence climbs steeply, then saturates at ~0.98 (not 1.0)
+
+Mean pairwise Hungarian-matched |corr| across the 3 backends
+(`benchmarks/figures/ksweep_70ch.png`):
+
+| frames  |     k | mean \|corr\| | min pair | pairwise comps > 0.95 |
+|--------:|------:|--------------:|---------:|----------------------:|
+|  73,500 |    15 |        0.9105 |   0.8656 |                 55.2% |
+| 147,000 |    30 |        0.9288 |   0.8931 |                 56.7% |
+| 294,000 |    60 |    **0.9818** |   0.9717 |                 90.0% |
+| 490,000 |   100 |        0.9826 |   0.9790 |                 94.8% |
+| 747,750 | 152.6 |        0.9817 |   0.9712 |                 92.4% |
+
+The hypothesis holds **directionally** with an important nuance: equivalence rises sharply from
+0.91 (k=15) to 0.98 (k=60) and then **plateaus at ~0.982**, it does not keep climbing to 1.0.
+The transition is complete by **k ~= 60** (frames ~= 3 * channels^2); the sweep samples no point
+between k=30 (0.929) and k=60 (0.982), so the exact knee within that gap is unresolved (it is an
+upper bound, not a pinpoint). Below it the decomposition is under-determined and backends settle
+into different (equally valid) local optima on the weak components; above it the
+strong-and-medium components lock and only the residual precision / init spread remains.
+
+## Why the plateau sits at 0.98, not 1.0
+
+The per-config equivalence matrix at k=152 (`benchmarks/figures/ksweep_70ch_matrix.png`) splits
+the residual cleanly by precision:
+
+| pair | \|corr\| @ k=152 |
+|---|---|
+| native-fortran-f64 vs torch-cuda-**f64** | **0.995** |
+| native-fortran-f64 vs torch-cuda-f32     | 0.971 |
+| torch-cuda-f64 vs torch-cuda-f32         | 0.979 |
+
+- **The two double-precision implementations agree at 0.995** -- an independent Fortran binary
+  and the PyTorch-CUDA backend recover essentially the same ICs once k is large. This is the
+  strongest drop-in-replacement evidence in the epic: cross-*implementation*, not just
+  cross-device.
+- The residual mean-gap below 1.0 is dominated by **float32**, for two compounding reasons:
+  (1) f32 rounding accumulates over the larger data and 2000 iters; (2) at k=152 the f32 run
+  hit the natural-gradient **lrate floor (1e-12) at iter 1735 and stopped early**, landing on a
+  marginally different optimum than the full-2000-iter f64 runs. Both are valid converged
+  solutions; this is why Phase 3's "f32 == f64 at 1.000" (k=30, both ran the full 2000) relaxes
+  to ~0.98 on the full recording. It is a convergence/precision effect, not a backend defect.
+
+## Practical guideline (the useful deliverable)
+
+For **backend-reproducible** ICs at 70ch, budget **k ~= 60 (frames ~= 3 * channels^2)**, roughly
+2x the EEGLAB k >= 20-30 rule. The classic rule guarantees only that *the strongest* components
+are stable; k ~= 60 is where the *bulk* of the decomposition (90%+ of pairwise component
+correlations > 0.95) becomes reproducible across independent AMICA implementations. Past that,
+more data buys convergence robustness but not more cross-backend agreement -- the ceiling is set
+by float precision and the non-convex local-optimum spread, not by data adequacy.
+
+## Per-config decompose time (2000 iters, 70ch, hallu)
+
+Time scales ~linearly with frames, as expected; native-fortran-f64 is fastest per config
+(torch-cuda ~1.4-1.5x), and torch-cuda-f32 undercuts f64 (and stops early at k=152).
+
+| frames  |     k | fortran-f64 | cuda-f64 | cuda-f32 |
+|--------:|------:|------------:|---------:|---------:|
+|  73,500 |    15 |      122.6s |   185.0s |   175.0s |
+| 147,000 |    30 |      248.1s |   366.4s |   344.7s |
+| 294,000 |    60 |      502.3s |   730.8s |   679.3s |
+| 490,000 |   100 |      848.2s |  1216.8s |  1125.7s |
+| 747,750 | 152.6 |     1303.1s |  1856.2s |  1496.5s |
+
+(Per-sample normalized final LL falls monotonically with more data, -3.285 at k=15 to -3.699 at
+k=152; that is data heterogeneity -- a longer recording spans more brain states the single model
+must explain -- not a worse fit.)
+
+## Figures
+- `benchmarks/figures/ksweep_70ch.png` -- mean cross-backend |corr| vs k (the headline: steep
+  climb to the k~=60 knee, then a flat ~0.982 plateau below the 1.0 line).
+- `benchmarks/figures/ksweep_70ch_matrix.png` -- per-config equivalence matrix at k=152,
+  showing the f64-vs-f64 = 0.995 cell and the f32 rows at ~0.97-0.98.
+
+## Note on the harness
+`--compare` writes the k-sweep line plot to `--figure` and the per-config matrix to a derived
+`{stem}_matrix{suffix}` path, so the two no longer overwrite each other (the first pass of this
+run clobbered the line plot with the matrix; fixed in `_compare`).

benchmarks/benchmark_decompose.py

@@ -233,6 +233,81 @@ def _variance_order(W, data, channels):
     return np.argsort(-projvar)
 
 
+def _k_sweep(runs, figure=None):
+    """Data-size hypothesis test (#90): at fixed channels, mean cross-backend
+    Hungarian-matched |corr| vs k = frames/ch^2. Equivalence should rise toward
+    1.0 as k grows (the decomposition becomes well-determined)."""
+    nc = runs[0]["channels"]
+    by_f: dict = {}
+    for r in runs:
+        if r["frames"] <= 0:  # legacy npz without a recorded frame count -> skip
+            continue
+        by_f.setdefault(r["frames"], []).append(r)
+    rows, lines = [], []
+    for nf in sorted(by_f):
+        g = by_f[nf]
+        if len(g) < 2:
+            continue
+        pair_means, all_c = [], []
+        for i in range(len(g)):
+            for j in range(i + 1, len(g)):
+                c = _match_correlation(g[i]["W"], g[j]["W"])
+                pair_means.append(c.mean())
+                all_c.append(c)
+        allc = np.concatenate(all_c)
+        k = nf / nc**2
+        mean_eq = float(np.mean(pair_means))
+        rows.append((nf, k, mean_eq))
+        lines.append(
+            f"  {nf:8d}  {k:6.1f}    {mean_eq:.4f}   {min(pair_means):.4f}"
+            f"     {(allc > 0.95).mean() * 100:5.1f}%   {len(g)}"
+        )
+    if not rows:  # every frame count had <2 backends -- nothing to compare
+        print(
+            f"\n  data-size (k) sweep @ {nc} channels: need >=2 backends at a shared "
+            "frame count -- nothing to compare"
+        )
+        return
+    print(f"\n=== data-size (k) sweep @ {nc} channels ===")
+    print("  frames        k   mean|corr|   min    comps>0.95   n")
+    for ln in lines:
+        print(ln)
+    if figure:
+        _plot_ksweep(rows, nc, figure)
+
+
+def _plot_ksweep(rows, channels, path):
+    import matplotlib
+
+    matplotlib.use("Agg")
+    import matplotlib.pyplot as plt
+
+    ks = [r[1] for r in rows]
+    eq = [r[2] for r in rows]
+    fig, ax = plt.subplots(figsize=(6, 4))
+    ax.plot(ks, eq, "o-", color="#2a7")
+    for nf, k, e in rows:
+        ax.annotate(
+            f"{nf // 1000}k",
+            (k, e),
+            fontsize=7,
+            xytext=(0, 6),
+            textcoords="offset points",
+            ha="center",
+        )
+    ax.axhline(1.0, ls="--", lw=0.8, color="0.6")
+    ax.set_xlabel("k = frames / channels^2")
+    ax.set_ylabel("mean cross-backend |correlation|")
+    ax.set_title(
+        f"AMICA cross-backend IC equivalence vs data size @ {channels} channels\n"
+        "(more data -> better-determined decomposition -> higher equivalence)"
+    )
+    ax.grid(alpha=0.3)
+    fig.tight_layout()
+    fig.savefig(path, dpi=150)
+    print(f"wrote {path}")
+
+
 def _compare(dirs, figure=None, montage=None, topo_figure=None, n_topo=6, data=None):
     full = np.load(data).astype(np.float64) if data else None
     runs = []
@@ -243,6 +318,7 @@ def _compare(dirs, figure=None, montage=None, topo_figure=None, n_topo=6, data=N
                 {
                     "label": str(z["label"]),
                     "channels": int(z["channels"]),
+                    "frames": int(z["frames"]) if "frames" in z else -1,
                     "time": float(z["time"]),
                     "final_ll": float(z["final_ll"]),
                     "W": z["W"],
@@ -253,12 +329,39 @@ def _compare(dirs, figure=None, montage=None, topo_figure=None, n_topo=6, data=N
         print("no result npz found")
         return
 
-    # equivalence is per channel-count (W shapes must match); use the largest
-    by_ch: dict = {}
+    # data-size (k) sweep (#90): at the top channel count, if multiple frame counts
+    # were run, report mean cross-backend equivalence vs k = frames/ch^2. When a
+    # k-sweep runs it owns --figure (the headline plot); the per-config matrix below
+    # then goes to a derived "{stem}_matrix{suffix}" path so neither overwrites the
+    # other.
+    top_ch = max(r["channels"] for r in runs)
+    top_runs = [r for r in runs if r["channels"] == top_ch]
+    # legacy npz (pre-#90) carry no frame count (frames == -1); they can't join a
+    # k-sweep and won't group with per-frame npz in the matrix below, so warn rather
+    # than silently group them wrong when the two formats are merged in one --compare.
+    if any(r["frames"] < 0 for r in top_runs) and any(
+        r["frames"] > 0 for r in top_runs
+    ):
+        print(
+            f"  warning: {top_ch}ch mixes legacy npz (no recorded frame count) with "
+            "per-frame npz; legacy runs are excluded from the k-sweep and grouped "
+            "separately in the equivalence matrix"
+        )
+    real_frames = {r["frames"] for r in top_runs if r["frames"] > 0}
+    matrix_figure = figure
+    if len(real_frames) > 1:
+        _k_sweep(top_runs, figure)
+        if figure:
+            p = Path(figure)
+            matrix_figure = str(p.with_name(f"{p.stem}_matrix{p.suffix}"))
+
+    # single-config equivalence matrix + topomaps: use the largest (channels, frames)
+    by_cf: dict = {}
     for r in runs:
-        by_ch.setdefault(r["channels"], []).append(r)
-    target = max(by_ch)
-    group = by_ch[target]
+        by_cf.setdefault((r["channels"], r["frames"]), []).append(r)
+    target_cf = max(by_cf)
+    target = target_cf[0]
+    group = by_cf[target_cf]
     labels = [r["label"] for r in group]
     n = len(group)
 
@@ -289,10 +392,14 @@ def _compare(dirs, figure=None, montage=None, topo_figure=None, n_topo=6, data=N
         f"(min {offdiag.min():.4f}) -- ~1.0 => same ICs recovered"
     )
 
-    if figure:
-        _plot(labels, M, target, figure)
+    if matrix_figure:
+        _plot(labels, M, target, matrix_figure)
     if topo_figure and montage:
-        _plot_topomaps(group, montage, target, topo_figure, n_topo, full)
+        # de-sphere must recompute the sphere from the SAME frames the group was fit
+        # on (#90 --frames can truncate); pass the matching slice, not the full array.
+        nf = target_cf[1]
+        topo_data = full[:, :nf] if (full is not None and nf > 0) else full
+        _plot_topomaps(group, montage, target, topo_figure, n_topo, topo_data)
 
 
 def _load_info(montage_tsv, n_channels):
@@ -449,6 +556,12 @@ def main() -> int:
     ap = argparse.ArgumentParser(description=__doc__)
     ap.add_argument("--data", help="path to the (channels, frames) real-EEG .npy")
     ap.add_argument("--channels", default="70")
+    ap.add_argument(
+        "--frames",
+        default=None,
+        help="data-size sweep (#90): comma-separated frame counts (first N samples). "
+        "Default: all frames. At fixed channels this sweeps k = frames/ch^2.",
+    )
     ap.add_argument("--iters", type=int, default=2000)
     ap.add_argument("--threads", type=int, default=None)
     ap.add_argument("--backends", default="auto")
@@ -494,33 +607,42 @@ def main() -> int:
     tag = _platform_tag()
     print(f"platform {tag} | channels {channels} | iters {args.iters} | {backends}")
 
+    # clamp to the available frames and dedupe up front so two requested counts that
+    # collapse to the same length don't re-fit identical data and clobber each other's npz.
+    frame_list = (
+        sorted({min(int(f), full.shape[1]) for f in args.frames.split(",")})
+        if args.frames
+        else [full.shape[1]]
+    )
     for nc in channels:
-        data = np.ascontiguousarray(full[:nc, :])
-        k = data.shape[1] / nc**2
-        print(f"\n{nc}ch, {data.shape[1]} frames (k={k:.1f})")
-        for b in backends:
-            label = f"{b}@{tag}"
-            try:
-                t0 = perf_counter()
-                res = _fit(b, data, args.iters, args.fortran_bin, args.threads)
-                fname = out / f"{b}_{nc}ch_{tag}.npz"
-                np.savez_compressed(
-                    fname,
-                    label=label,
-                    backend=b,
-                    channels=nc,
-                    iters=args.iters,
-                    time=res["time"],
-                    final_ll=res["final_ll"],
-                    W=res["W"],
-                    A=res["A"],
-                )
-                print(
-                    f"  {b:24s} {res['time']:8.1f} s   LL={res['final_ll']:+.5f}"
-                    f"   [{perf_counter() - t0:.0f}s wall]  -> {fname.name}"
-                )
-            except Exception as exc:  # noqa: BLE001 - report and continue
-                print(f"  {b:24s} FAILED: {type(exc).__name__}: {str(exc)[:80]}")
+        for nf in frame_list:
+            data = np.ascontiguousarray(full[:nc, :nf])
+            k = nf / nc**2
+            print(f"\n{nc}ch, {nf} frames (k={k:.1f})")
+            for b in backends:
+                label = f"{b}@{tag}"
+                try:
+                    t0 = perf_counter()
+                    res = _fit(b, data, args.iters, args.fortran_bin, args.threads)
+                    fname = out / f"{b}_{nc}ch_{nf}f_{tag}.npz"
+                    np.savez_compressed(
+                        fname,
+                        label=label,
+                        backend=b,
+                        channels=nc,
+                        frames=nf,
+                        iters=args.iters,
+                        time=res["time"],
+                        final_ll=res["final_ll"],
+                        W=res["W"],
+                        A=res["A"],
+                    )
+                    print(
+                        f"  {b:24s} {res['time']:8.1f} s   LL={res['final_ll']:+.5f}"
+                        f"   [{perf_counter() - t0:.0f}s wall]  -> {fname.name}"
+                    )
+                except Exception as exc:  # noqa: BLE001 - report and continue
+                    print(f"  {b:24s} FAILED: {type(exc).__name__}: {str(exc)[:80]}")
     return 0