neuromechanist · Jul 7, 2026 · Jul 7, 2026
diff --git a/.context/feature_parity.md b/.context/feature_parity.md
@@ -0,0 +1,97 @@
+# Feature Parity Report: Fortran vs PyTorch AMICA
+
+Compares the Fortran reference (`amica15.f90`, binary `amica15mac`) with the current PyTorch
+backend `AMICATorchNG` (the sole PyTorch backend; the earlier Adam/autograd `AMICATorch` and
+`AMICATorchV2` were removed in #32). The NumPy reference (`AMICA_NumPy`) is retained as an oracle
+and CLI. Living status is in `../AGENTS.md` (Known Issues); ADRs are in `decisions/`.
+
+**Headline:** single-model parity is done (LL ~ -3.40 vs Fortran -3.4018; Hungarian component
+correlation ~0.997 with Newton positive-definite, 0 fallbacks). Adaptive PDF (#26), multi-model
+(#27, distributional equivalence), best-iterate (#51), and the degenerate-fit contract (#50) are
+done. The one unimplemented Fortran feature is component sharing; the one unverified success
+criterion is runtime.
+
+## Core AMICA features
+
+| Feature | Fortran | AMICATorchNG | Status | Notes |
+|---------|---------|--------------|--------|-------|
+| Multiple models | yes | yes | Complete | `n_models`; per-model bias `c` update (#27) |
+| Mixture of Gaussians | yes | yes | Complete | `num_mix` components |
+| Generalized Gaussian PDF | yes | yes | Complete | default `pdftype=0`, byte-for-byte unchanged |
+| Natural gradient | yes | yes | Complete | A-update parity fix (#24) |
+| Learning-rate adaptation | yes | yes | Complete | Fortran-matched lrate schedule |
+| Convergence criteria | yes | yes | Complete | dLL and grad-norm checks |
+| Best-iterate return | yes | yes | Complete | `keep_best` / `final_ll_` (#51) |
+
+## Optimization methods
+
+| Feature | Fortran | AMICATorchNG | Status | Notes |
+|---------|---------|--------------|--------|-------|
+| Natural gradient | yes | yes | Complete | exact-EM mixture updates |
+| Newton method | yes | yes | Complete | positive-definite, 0 fallbacks on sample data |
+| Newton ramping | yes | yes | Complete | Fortran-style ramp |
+
+## Data preprocessing
+
+| Feature | Fortran | AMICATorchNG | Status | Notes |
+|---------|---------|--------------|--------|-------|
+| Mean removal | yes | yes | Complete | `do_mean` |
+| Sphering (whitening) | yes | yes | Complete | symmetric-ZCA (#24) |
+| PCA dimension reduction | yes | yes | Complete | `pcakeep` |
+
+## Advanced features
+
+| Feature | Fortran | AMICATorchNG | Status | Notes |
+|---------|---------|--------------|--------|-------|
+| Outlier rejection | yes | yes | Complete | `do_reject` |
+| Adaptive PDF selection | yes | yes | Complete | 5 `pdftype` families + ext-Infomax switcher (#26) |
+| Component sharing | yes | **no** | **Missing** | `share_comps` reassignment not ported (torch `comp_list` is trivial per-model layout) |
+| History tracking | yes | yes | Complete | `ll_history`, grad norms |
+
+## Input/output
+
+| Feature | Fortran | AMICATorchNG | Status | Notes |
+|---------|---------|--------------|--------|-------|
+| Binary data loading | yes | yes | Complete | `.fdt` files |
+| Parameter files | yes | yes | Complete | JSON |
+| Checkpoint save/load | yes | yes | Complete | NG state_dict persistence (#36); untested (#15) |
+| Result loading | yes | yes | Complete | via `loadmodout()` |
+
+## Numerical stability
+
+| Feature | Fortran | AMICATorchNG | Status | Notes |
+|---------|---------|--------------|--------|-------|
+| Log-space computation | yes | yes | Complete | for stability |
+| Eigenvalue thresholds | yes | yes | Complete | `min_eig` |
+| NaN/Inf handling | yes | yes | Complete | degenerate-fit contract refuses NaN output (#50) |
+| Min/max bounds regression tests | yes | yes | **Untested** | mincond/minlog/maxdble/mineig present, no regression tests |
+
+## Performance and hardware
+
+| Feature | Fortran | AMICATorchNG | Notes |
+|---------|---------|--------------|-------|
+| OpenMP parallelization | yes | n/a | PyTorch uses tensor ops / device parallelism |
+| GPU support | no | yes | CUDA / MPS / CPU automatic (float64 parity runs on CPU/CUDA) |
+| Runtime vs Fortran | baseline | **unmeasured** | the 2-3x criterion has never been benchmarked |
+
+## Validation status (issue #24)
+
+| Metric | Fortran | AMICATorchNG | Match |
+|--------|---------|--------------|-------|
+| Final LL (single model) | -3.4018 | ~ -3.40 | yes (Jacobian-normalized) |
+| Component correlation | - | ~0.997 | yes (Hungarian-matched, clears >0.95) |
+| Newton | posdef | posdef | yes (0 fallbacks on sample data) |
+| Multi-model per-block sufficient stats | - | bit-exact | yes |
+
+Multi-model full-fit partition cross-correlation (~0.64 single-run) is intrinsic estimator spread,
+not a defect: Fortran agrees with itself at ~0.63, and the NG partition ensemble is statistically
+indistinguishable from Fortran's run-to-run distribution (#27, `.context/issue-27/`).
+
+## Remaining roadmap
+
+- [ ] **Performance benchmark** vs the Fortran binary (CPU/CUDA/MPS); verify or revise the 2-3x
+      runtime criterion.
+- [ ] **Component sharing** (`share_comps`, `share_start`/`share_int`, `load_comp_list`,
+      `comp_used`) ported to `AMICATorchNG` and validated on real data.
+- [ ] **Test hardening:** single-channel / single-sample edge cases, numerical-stability
+      regression tests, and `save`/`load` + `plot_components` coverage (#15).
diff --git a/.context/migration_plan.md b/.context/migration_plan.md
@@ -0,0 +1,40 @@
+# Migration Plan: NumPy to PyTorch (COMPLETE)
+
+**Status: complete.** This file is kept as a historical record. The NumPy-to-PyTorch migration
+landed via the natural-gradient EM backend (`AMICATorchNG`); the PyTorch path is the default and
+reaches Fortran parity. The NumPy implementation is retained as an oracle and CLI, not as a
+fallback. Current status lives in `../AGENTS.md` (Known Issues) and `progress_summary.md`.
+
+## What the migration delivered
+
+- **Default backend is PyTorch.** `from pyAMICA import AMICA` wraps `AMICATorchNG`. The legacy
+  reference is exposed separately as `AMICA_NumPy` and is not selected via a `backend=` flag; there
+  is no `legacy/` directory and no runtime backend switch.
+- **Feature parity.** Natural gradient + Newton (positive-definite, with Fortran-style ramping),
+  exact-EM mixture updates, symmetric-ZCA sphere, Jacobian LL, outlier rejection (`do_reject`), all
+  five `pdftype` density families (#26), and multi-model support with the per-model bias `c` update
+  (#27). Single-model results match the Fortran fixed point (LL ~ -3.40, component correlation
+  ~0.997).
+- **Device support.** Automatic CUDA / MPS / CPU selection. The NG backend computes in float64 for
+  Fortran parity, which MPS cannot represent, so parity runs use CPU or CUDA; the `AMICA` wrapper
+  falls back to CPU automatically when a device is not pinned. Run with
+  `PYTORCH_ENABLE_MPS_FALLBACK=1` for ops MPS does not support.
+- **Environment.** Migrated from the conda `torch-312` env to UV; the stack is declared in
+  `pyproject.toml` with `uv.lock` pinned and `.python-version` set to 3.12.
+- **Validation.** `validate_implementations.py` runs both backends on real sample EEG plus the
+  Fortran binary and matches components via the Hungarian algorithm. No mocks or synthetic data.
+
+## Not carried over from the original plan
+
+- **`backend='numpy'` runtime switch / `legacy/` dir / gradual dual-backend transition.** Dropped.
+  The two implementations are separate public classes, not a runtime toggle.
+- **Adam / L-BFGS / autograd "bonus" optimizers.** The basic `AMICATorch` and `AMICATorchV2`
+  (Adam/autograd) backends were removed in #32; `AMICATorchNG` (natural-gradient EM) is the only
+  PyTorch backend.
+
+## Still open (tracked outside this plan)
+
+- Performance benchmark vs Fortran (the "2-3x runtime" criterion is unmeasured).
+- Component sharing (`share_comps`), the last unimplemented Fortran feature.
+- Edge-case + numerical-stability regression tests; `save`/`load` and `plot_components` coverage
+  (issue #15).
diff --git a/.context/plan.md b/.context/plan.md
@@ -8,13 +8,15 @@
 ## Development Tasks
 <!-- Status markers: [ ] pending, [~] in progress, [x] complete -->
 
-### Priority 1: Parity blockers
-- [~] Fix likelihood sign/scaling vs Fortran (~13x factor; GG normalization corrected, needs revalidation)
-- [~] Stabilize Newton optimization (NaN at Newton-start iter; clipping + Fortran-matched ramp added, needs revalidation)
-- [ ] Add numerical-stability bounds from Fortran (`mincond=1e-15`, `minlog=-1500`, `maxdble=1e32`, `mineig=1e-15`)
-- [ ] Add NaN/Inf checks and epsilon to all divisions (notably `dmu / dalpha`)
-- [ ] Ensure identical initialization vs Fortran (seed, sphering/whitening, starting mixing matrix)
-- [ ] Raise component correlation with Fortran (~0.46-0.9, run-dependent) to >0.95
+### Priority 1: Parity blockers — DONE (#24)
+- [x] Fix likelihood sign/scaling vs Fortran (the ~13x factor was a pre-parity basic-backend
+      artifact; NG uses Jacobian LL, single-model LL ~ -3.40 vs Fortran -3.4018)
+- [x] Stabilize Newton optimization (ported from NumPy, positive-definite, 0 fallbacks on sample data)
+- [x] Add numerical-stability bounds from Fortran (`mincond`, `minlog`, `maxdble`, `mineig` present;
+      regression tests still owed, see Priority 3)
+- [x] Add NaN/Inf checks (degenerate-fit contract #50 refuses NaN output instead of suppressing it)
+- [x] Ensure identical initialization vs Fortran (symmetric-ZCA sphere, seed, starting mixing matrix)
+- [x] Raise component correlation with Fortran to >0.95 (now ~0.997, Hungarian-matched)
 
 ### Priority 2: Missing core features
 - [x] Port outlier rejection (`do_reject` / `reject_data`) to the torch backend (done in `AMICATorchNG`)
@@ -47,23 +49,26 @@
       suite collects cleanly and passes (35 passed, 6 xfailed for documented parity/algorithm
       issues, 0 errors); see `.context/phase1_baseline.md`
 - [x] Integration tests comparing against Fortran outputs (`tests/torch_tests/`) - Phase 1, issue #10
-- [ ] Numerical-stability regression tests
+- [ ] Numerical-stability regression tests (mincond/minlog/maxdble/mineig)
 - [ ] Edge cases (single channel, single sample)
+- [ ] `AMICA.save`/`load` and `plot_components` coverage (issue #15)
+- [ ] Performance benchmark: NG runtime vs Fortran binary (verify the 2-3x criterion)
 
 ### Infrastructure / migration
 - [x] Migrate environment from conda `torch-312` to UV; declare the PyTorch stack in `pyproject.toml`
       (Phase 1, issue #10: `torch`, `pytest`, `pytest-cov` added to `pyproject.toml`/`uv.lock`,
       `.python-version` pinned to 3.12 for torch/MPS compatibility)
-- [ ] Set up CI (see `.rules/ci_cd.md`)
+- [x] Set up CI (see `.rules/ci_cd.md`): ruff lint/format, pytest (excluding slow/Fortran-binary
+      parity), build + clean-env import matrix on Python 3.12/3.13; typos check. Green on `main`.
 - [x] Fix legacy test function signatures (`load_data_file` extra parameter) - Phase 1, issue #10
 
 ## Success Criteria
-- [ ] Component correlation > 0.95 with Fortran
-- [ ] LL convergence within 1% of Fortran
-- [ ] No NaN/Inf during optimization
-- [ ] Runtime within 2-3x of Fortran
-- [ ] Real (non-mock) test suite green
+- [x] Component correlation > 0.95 with Fortran (~0.997, single model)
+- [x] LL convergence within 1% of Fortran (LL ~ -3.40 vs -3.4018)
+- [x] No NaN/Inf during optimization (degenerate-fit contract #50)
+- [ ] Runtime within 2-3x of Fortran (unmeasured; benchmark owed)
+- [x] Real (non-mock) test suite green
 
 ## Notes
 - Correctness is defined by parity with the Fortran binary, not by convergence alone.
-- Detailed feature status: `../FEATURE_PARITY.md`; migration timeline: `../MIGRATION_PLAN.md`.
+- Detailed feature status: `feature_parity.md`; migration record: `migration_plan.md`.
diff --git a/.context/progress_summary.md b/.context/progress_summary.md
@@ -0,0 +1,68 @@
+# Progress Summary
+
+Internal status snapshot for the PyTorch AMICA effort. The living source of truth for parity
+status is the **Known Issues** section of `../AGENTS.md`; this file summarizes what has landed and
+what remains as of the v0.1.0 preparation.
+
+## Delivered
+
+### Natural-gradient EM backend at Fortran parity (epic #9 / issue #24)
+- `AMICATorchNG` (`pyAMICA/torch_impl/core.py`) is the sole PyTorch backend; the `AMICA`
+  scikit-learn-style wrapper wraps it directly.
+- Root cause of the historical parity gap: the natural-gradient A-update was transposed and
+  multiplied on the wrong side. Fix plus exact-EM mixture updates, the digamma rho update, the
+  symmetric-ZCA sphere, the output transpose, and the Jacobian log-likelihood bring single-model
+  results to the Fortran fixed point.
+- Newton is ported from the NumPy reference and stays positive-definite (0 fallbacks on the
+  sample data).
+- **Metrics:** single-model LL ~ -3.40 (Fortran -3.4018); Hungarian-matched component correlation
+  ~0.997, clearing the >0.95 gate.
+
+### Adaptive PDF selection (issue #26)
+- All five `amica15.f90` `pdftype` density families ported to `AMICATorchNG`: 0 generalized
+  Gaussian (default, byte-for-byte unchanged), 2 Gaussian, 3 logistic, 4 sub-Gaussian cosh+, and
+  the extended-Infomax kurtosis switcher (`pdftype=1`).
+- Fixed families are bit-exact vs the literal Fortran `z0`/`fp` (~1e-15) and converge within
+  ~0.005 LL of the binary. The dynamic switcher is dead code in the binary (no bit-exact oracle),
+  so it is validated by real-data LL. See ADR 0002.
+
+### Multi-model AMICA (issue #27)
+- Validated by **distributional equivalence**: multi-model AMICA is not partition-identifiable, so
+  the NG-vs-Fortran partition cross-correlation distribution is statistically indistinguishable
+  from Fortran's own run-to-run spread (Mann-Whitney p=0.97, TOST within +/-0.05). Per-block
+  sufficient stats are bit-exact vs Fortran.
+- Per-model bias `c` update (`update_c`) ported to both backends, guarded to a no-op for
+  `n_models=1` so single-model parity stays bit-exact. See `.context/issue-27/`.
+
+### Best-iterate safeguard (issue #51)
+- `AMICATorchNG.fit` returns the highest-LL iterate (`keep_best`, default on; `final_ll_` reports
+  the returned iterate's LL, `ll_history` keeps the true trajectory), not the last iterate under
+  the non-monotone lrate schedule. Cuts multi-model LL sd from 12.7x to 2.0x Fortran's at a
+  matched 100-iter budget. Single-model #24 parity stays bit-exact (monotone => no restore). See
+  ADR 0003.
+
+### Degenerate-fit contract (issue #50)
+- The `AMICA` wrapper no longer treats a degenerate fit (`stop_reason` nan_ll / singular_ll) as
+  usable. `fit` sets `is_fitted_` only on a converged fit and exposes `converged_` /
+  `stop_reason_`; `transform` / `get_mixing_matrix` / `get_unmixing_matrix` / `save` raise a clear
+  degenerate error instead of returning NaN sources.
+
+### Structure and infrastructure
+- Module rename into `numpy_impl/` and `torch_impl/` with topic-based names (issue #34); the public
+  import surface (`from pyAMICA import AMICA, AMICA_NumPy, AMICATorchNG`) is stable.
+- UV is the canonical environment (`pyproject.toml` + `uv.lock`); the legacy conda env is retired.
+- CI is live and green on `main`: ruff lint/format, pytest (excluding slow/Fortran-binary parity),
+  and a build + clean-env import matrix on Python 3.12 and 3.13. Typos check green.
+- Validation harness (`validate_implementations.py`) runs both backends and matches components via
+  the Hungarian algorithm, on real sample EEG plus the Fortran binary (`amica15mac`).
+
+## Remaining before / around v0.1.0
+
+- **Performance benchmark:** the "runtime within 2-3x of Fortran" success criterion has never been
+  measured. Needs a repeatable NG-vs-Fortran benchmark across CPU/CUDA/MPS.
+- **Component sharing:** the last unimplemented Fortran feature. Torch `comp_list` is only the
+  trivial per-model block layout; `share_comps` reassignment (`share_start` / `share_int`,
+  `load_comp_list`, `comp_used`) is absent. NumPy has partial plumbing only.
+- **Test hardening:** no single-channel / single-sample edge tests, no numerical-stability
+  regression tests (mincond / minlog / maxdble / mineig bounds); `AMICA.save`/`load` and
+  `plot_components` remain untested (issue #15).
diff --git a/AGENTS.md b/AGENTS.md
@@ -58,7 +58,8 @@ computes in float64 for Fortran parity, which MPS cannot represent, so parity ru
 - Newton and exact-EM updates are implemented in `AMICATorchNG` and the legacy NumPy `numpy_impl/core.py`
   (both Fortran-faithful). Adaptive PDF (#26) is DONE (all five `pdftype` families + ext-Infomax
   switcher); full multi-model matching (#27) is validated by distributional equivalence.
-- See `FEATURE_PARITY.md`, `MIGRATION_PLAN.md`, and `PROGRESS_SUMMARY.md` for detailed roadmaps.
+- See `.context/feature_parity.md`, `.context/migration_plan.md`, and `.context/progress_summary.md`
+  for detailed roadmaps.
 
 ## Known Issues (parity blockers)
 **Single-model parity: DONE (#24).** The natural-gradient A-update transpose fix (plus exact-EM
@@ -149,12 +150,12 @@ byte-for-byte unchanged. See `.context/decisions/` and `tests/torch_tests/test_n
 - `.context/research.md` - Fortran-vs-Python parity analysis (data structures, subroutines, divergence causes)
 - `.context/ideas.md` - PyTorch design decisions and library options
 - `.context/scratch_history.md` - Debugging notes, failed attempts, lessons
+- `.context/feature_parity.md` - Feature comparison and implementation roadmap
+- `.context/migration_plan.md` - NumPy to PyTorch migration record (complete)
+- `.context/progress_summary.md` - Achievements and validation metrics snapshot
 - `.context/decisions/` - Architecture Decision Records (copy `0000-template.md` to start one)
 
 ## Project Docs (top-level)
-- `FEATURE_PARITY.md` - Feature comparison and implementation roadmap
-- `MIGRATION_PLAN.md` - NumPy to PyTorch migration timeline
-- `PROGRESS_SUMMARY.md` - Achievements and validation metrics snapshot
 - `README.md` - Overview and quick start
 
 ---