Add per-model bias c update for multi-model AMICA (#27) (#49)

* Add per-model bias c update for multi-model AMICA

Port Fortran's update_c (amica17.f90:1423-1429/1899-1901) into AMICATorchNG
and the NumPy oracle: c[i,h] = sum_t v_h*x / sum_t v_h, the per-model
responsibility-weighted data-space mean. The E-step now centers each model's
data before unmixing (b = W(x - c)); transform() does the same. Replaces the
old gradient-style dc = sum(g) accumulator, which was accumulated but never
applied (c was frozen at 0).

Guarded to a no-op for n_models=1: with v==1 the update collapses to the
(zero) mean of mean-removed data, and skipping it keeps single-model parity
bit-exact (issue #24).

Tests: single-model c stays exactly zero after fit; multi-model c equals the
responsibility-weighted data mean and the two models center differently;
existing NG<->NumPy sufficient-stat parity holds with dc renamed to dc_numer.

Controlled 2-model A/B vs the Fortran binary (same config/seed, c toggled):
cross-corr 0.631 -> 0.642 (+0.011), LL unchanged. The c omission was a minor
contributor; the dominant multi-model gap is intrinsic partition ambiguity
(see .context/issue-27/multimodel_c_update.md). Issue #27.

* Address PR review: guard dead-model c, add tests

Review findings (pr-review-toolkit, Sonnet):

- silent-failure: the new c = dc_numer/dgm division could be 0/0 = NaN for a
  dead model (dgm[h]==0). Unlike log(gm[h])=-inf (which softmax tolerates), a
  NaN c poisons the next iteration's cross-model softmax for every model. Added
  a containment guard in both backends: a zero-responsibility model keeps its
  prior c, mirroring the existing mu/beta/rho non-finite guards. This also
  resolves the NumPy restart-preserves-NaN-c concern.

- tests: added multi-model coverage the change opened up but the first commit
  left unexercised -- NumPy backend c update on real data, NG<->NumPy finalized
  c parity, transform() with nonzero c (verified vs W(x-c) by hand), the
  dead-model containment guard, multi-model dc_numer blocking invariance, and
  do_reject + multi-model c finiteness. Strengthened the Newton multi-model test
  to assert finite c and full iteration count.

- comments: fixed two pre-existing docstrings this change made stale
  (transform()/get_weights() said "X^T @ W", now "(X-c)^T @ W"); corrected the
  Fortran citation for wc = W@c (:2178 get_unmixing_matrices); clarified update_c
  is a flag not a routine; AGENTS.md "LL unchanged" -> "LL comparable".

All fast suites green (48 passed). Single-model paths unchanged (guard is
n_models>1). Issue #27.

* Add multi-model distributional-equivalence validation

Records the parity confirmation for multi-model AMICA (issue #27). Because
mixture-of-ICA is not partition-identifiable, exact partition parity with
Fortran is the wrong acceptance bar; the right test is whether the two
implementations sample the same distribution over solutions.

On an N=20-each ensemble (real sample EEG, n_models=2, 100 iters), the
NG-vs-Fortran partition cross-corr distribution is statistically equivalent to
Fortran's own run-to-run distribution (Mann-Whitney p=0.97; TOST equivalent
within +/-0.05; within-Fortran/within-NG/between all ~0.63-0.64). The single-run
~0.64 cross-corr earlier read as a shortfall is just intrinsic estimator spread:
Fortran agrees with itself at 0.63.

Adds:
- .context/issue-27/multimodel_distributional_equivalence.md (method, results,
  acceptance criteria)
- multimodel_ensemble.py (reproduction harness) + the figure (PNG/PDF)
- research.md / AGENTS.md pointers; AGENTS.md #27 now reads VALIDATED

Open residual tracked as #51: NG's LL distribution is ~0.02 lower and more
variable than Fortran's (optimizer quality, not a correctness bug -- one M-step
is bit-exact). Issue #27.

Author: neuromechanist

.context/issue-27/multimodel_c_update.md

@@ -0,0 +1,61 @@
+# Issue #27: per-model bias `c` update + 2-model partition re-measurement
+
+## What was done
+Ported the update gated by Fortran's `update_c` flag (numerator accumulation at
+amica17.f90:1423-1429, division at :1899-1901) into both `AMICATorchNG` and the
+legacy NumPy `pyAMICA.py`:
+
+- `c[i,h] = dc_numer[i,h] / dc_denom[i,h]`, with
+  `dc_numer[i,h] = sum_t v_h(t) * x(i,t)` (sphered-data space) and
+  `dc_denom[i,h] = sum_t v_h(t) = dgm[h]`.
+- `c` is the per-model, per-channel responsibility-weighted **data-space** mean.
+  The E-step centers each model's data before unmixing: `b = W(x - c)`
+  (Fortran subtracts `wc = W c`). The Python E-step now computes
+  `b = (X - c).T @ W` (data-space subtraction, equivalent, keeps `c`'s
+  semantics identical to Fortran).
+- Guarded to a no-op for `n_models = 1`: with `v == 1` the update collapses to
+  the (zero) mean of mean-removed data, and skipping it keeps single-model
+  parity **bit-exact** (issue #24). Without the guard a ~1e-13 float-sum
+  residual would perturb the otherwise machine-exact single-model trajectory.
+- Replaced the old `dc = sum(g)` accumulator (gradient-style bias that was
+  accumulated but never applied — `c` was frozen at 0) with the data-space
+  `dc_numer`. `transform()` in both backends now unmixes as `W(x - c)`.
+- Dead-model containment: a model with zero total responsibility (`dgm[h]==0`)
+  would give `0/0`; the guard keeps that model's PRIOR `c` instead of writing a
+  NaN. A NaN `c` would poison the next iteration's cross-model `softmax` for
+  every model (unlike `log(gm[h])=-inf`, which `softmax` tolerates), so this
+  mirrors the existing mu/beta/rho non-finite guards in the same method.
+
+## Controlled re-measurement (real sample EEG, 2 models, 100 iters)
+`scratchpad/measure_multimodel_xcorr.py`: Fortran binary (num_models=2) run to
+convergence, then `AMICATorchNG` (n_models=2) run twice with an otherwise
+identical config/seed — once with the `c` update ON, once forced OFF (c
+re-zeroed each iteration to reproduce the pre-fix trajectory). Metric: stacked
+`2*NW`-component Hungarian |correlation| of the per-model unmixing rows vs
+Fortran (same orientation as the validated single-model harness).
+
+| run              | final LL | mean cross-corr | min cross-corr |
+|------------------|---------:|----------------:|---------------:|
+| NG, c OFF (pre)  |  -3.3754 |          0.6306 |         0.2647 |
+| NG, c ON  (fix)  |  -3.3762 |          0.6415 |         0.2729 |
+| Fortran          |  -3.3596 |            1.000 |          1.000 |
+
+**delta mean cross-corr = +0.011** (LL comparable, essentially unchanged).
+
+## Conclusion
+The omitted `c` update was a genuine, fixable contributor to the multi-model
+gap, but a **minor** one: with all else held fixed it lifts the 2-model
+partition cross-correlation by only ~0.011. The dominant residual gap is
+**intrinsic partition ambiguity** (mixture-of-ICA has many near-degenerate
+partitions; NG is self-consistent, cross-corr 1.0 across block sizes). The
+`>0.95` target is not reachable via the `c` fix alone. This matches issue #27's
+prior expectation; the fix is retained because it is Fortran-faithful and
+correct (unit-validated: `c` equals the responsibility-weighted data mean and
+the two models center differently), and it removes the last known
+non-ambiguity discrepancy in the multi-model M-step.
+
+Absolute cross-corr magnitudes depend on run configuration and random init
+(Fortran uses its own RNG; NG uses seed=0), so the ~0.63 baseline here is not
+identical to the ~0.77 figure previously noted under a different setup. The
+**controlled A/B delta** (+0.011, same config/seed, only `c` toggled) is the
+config-independent result.

.context/issue-27/multimodel_distributional_equivalence.md

@@ -0,0 +1,102 @@
+# Multi-model AMICA parity: distributional equivalence to Fortran (issue #27)
+
+**Bottom line.** For multi-model AMICA (`n_models > 1`), the natural-gradient
+PyTorch backend (`AMICATorchNG`) is validated against the Fortran reference at
+the level that is actually well-posed: its **ensemble of solutions is
+statistically indistinguishable from Fortran's**. A single run's partition
+cross-correlation (~0.64) is *not* a parity defect; it is the intrinsic spread
+of the estimator, and Fortran exhibits the same spread against itself.
+
+## Why a single cross-correlation number is the wrong test
+
+Single-model ICA is (essentially) identifiable: up to permutation/sign/scale
+there is one solution, so two correct implementations converge to it (we see
+~0.997 component correlation vs Fortran). Multi-model AMICA is a **mixture of
+ICA models with a soft partition** (responsibilities `v_h(t)`), which is **not
+partition-identifiable**: many near-degenerate partitions achieve essentially
+the same likelihood, and EM converges to whichever basin the random
+initialization sits in. Two correct runs, even of the *same* implementation from
+different seeds, are therefore *not expected* to produce the same partition.
+
+So "match Fortran's partition" is not a well-posed acceptance criterion. The
+well-posed question is the one a statistician asks of a stochastic estimator:
+**do the two implementations sample the same distribution over solutions?**
+
+## Method (real sample EEG, NO MOCK)
+
+Ran `N = 20` fits per implementation on the sample EEG (`n_models = 2`, 3 mixture
+components, 100 iterations, matched schedule: `lrate 0.05`, Newton from iter 50,
+`newtrate 1.0`). Both implementations are genuinely stochastic run-to-run
+(Fortran reseeds its init RNG from entropy; NG varies with `seed`). For each run
+we stored the stacked `2*32` unmixing components and the final log-likelihood.
+
+Metric: **Hungarian-matched mean |correlation|** of the stacked components (the
+same matching used for single-model parity, which quotients out component
+permutation, sign, scale, and model-label switching). From the 20+20 runs we
+formed three distributions of pairwise agreement:
+
+- `within-Fortran` — all Fortran-Fortran pairs (Fortran's intrinsic spread), n=190
+- `within-NG` — all NG-NG pairs (NG's intrinsic spread), n=190
+- `between` — all NG-Fortran pairs (cross-implementation agreement), n=400
+
+## Results
+
+![Ensemble distributions](multimodel_ensemble_distributions.png)
+
+| distribution | mean cross-corr | sd | range |
+|---|---:|---:|---|
+| within-Fortran | 0.6339 | 0.042 | [0.567, 0.772] |
+| within-NG | 0.6438 | 0.046 | [0.537, 0.798] |
+| between (NG-Fortran) | 0.6381 | 0.047 | [0.525, 0.938] |
+
+- **Mann-Whitney** (one-sided, H1: `between < within-Fortran`): **p = 0.973** —
+  no evidence cross-implementation agreement is *worse* than Fortran's own
+  run-to-run agreement (if anything it is marginally higher).
+- **TOST equivalence** of the mean cross-corrs within ±0.05: **p ≈ 1e-32 →
+  EQUIVALENT** (difference +0.0042).
+
+The three distributions lie on top of each other. **The partition behavior of
+`AMICATorchNG` is statistically equivalent to Fortran's.** The ~0.64 single-run
+cross-corr that earlier looked like a shortfall is fully explained: Fortran
+agrees with *itself* at 0.634.
+
+### One residual: the likelihood distribution (tracked as #51)
+
+| | mean LL (per sample-channel) | sd |
+|---|---:|---:|
+| Fortran | -3.3545 | 0.003 |
+| NG | -3.3738 | 0.040 |
+
+KS p ≈ 1e-5. NG's LL is ~0.019 lower on average and ~13x more variable (one seed
+converged to ~-3.55, a stuck run, driving most of the variance; Fortran reaches
+nearly the same LL every run despite different partitions — confirming those
+partitions are near-degenerate in likelihood).
+
+This is an **optimizer-quality** signal, not a model-correctness bug: the
+per-block sufficient statistics and one M-step are bit-exact vs Fortran
+(~1e-15), so the equations are right. The inflated variance (not a uniform mean
+shift) points to occasional convergence to slightly worse local optima —
+likely the iteration cap or a schedule mismatch. Tracked as issue **#51**.
+
+## Multi-model acceptance criteria (the definition of done)
+
+This replaces the aspirational `>0.95` single-run cross-corr in #27's title,
+which asks the algorithm to be more identifiable than it mathematically is:
+
+1. **Algebra:** per-block sufficient statistics and one M-step bit-exact vs
+   Fortran (~1e-15) and vs the NumPy oracle (~1e-8). *(Held — test suite.)*
+2. **Partition:** the NG solution ensemble is statistically equivalent to
+   Fortran's (`between` not worse than `within-Fortran`; TOST within a margin).
+   *(Held — this document.)*
+3. **Likelihood:** NG's LL distribution equivalent to Fortran's on the sample
+   ensemble. *(Open — issue #51; small residual, optimizer tuning.)*
+4. **Self-consistency:** fixed seed reproduces exactly; cross-corr 1.0 across
+   block sizes. *(Held — `test_blocking_invariance*`.)*
+
+## Reproduction
+
+`.context/issue-27/multimodel_ensemble.py` runs the ensemble (Fortran binary +
+NG) and writes `ensemble.npz`; `plot_ensemble.py` renders the figure. Both use
+the real sample data and the macOS Fortran binary (x86_64, runs under Rosetta).
+Absolute cross-corr magnitudes depend on config/seed; the **controlled
+within-vs-between comparison** is the config-independent result.