Parity testing improvements - #39
Merged
Merged
Conversation
Overhaul the parity suite to catch real numerical divergence and close coverage gaps found while reviewing the existing tests. - Forward-model parity: use known fixture params, not MATLAB's recovered fit values, so forward parity is isolated from fit recovery. - Tolerances: add fit_recovery_strict (atol 1e-6) for well-conditioned fits; keep the loose profile only for ill-conditioned 2CXM. - Noisy-data parity: store deterministic noisy curves + MATLAB fits in the baseline and compare Python's fit of the identical curve, gated per parameter on identifiability (unstable params reported, never asserted). - Baseline drift guard: regenerate the baseline in CI and fail if the committed one no longer matches current MATLAB. - T1 end-to-end map parity: new MATLAB generator + committed VFA fixture + Python pipeline comparison on physiologically-plausible voxels. - Runtime parity: wire up the previously-dead --run-runtime-parity options to real tests that gate on numerical parity and report wall-clock ratio. - Mask-collapse: a required parity check with <2 valid voxels now fails instead of silently skipping; multi-model suites assert real coverage. - CI: add T1 map parity + a fast tofts/patlak sparse-ROI DCE check. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Metrics (§3): gate on RMSE + Corr; drop MAE/p95; add reported-only CI-normalized abs-diff (p95) and proportion-outside-CI for every Python-vs-MATLAB parameter check (both sides emit 95% CI maps). Regions + gating (§2): evaluate brain/GM/WM (pipeline fits a union ROI so GM/WM voxels are fit); gate only Tofts+Patlak Ktrans. Tofts-GM is reported-only — its Python-vs-MATLAB disagreement is non-identifiability, not a bug (Python SSE equal-or-better in 49/50 GM voxels). Standard suite is default-on (test_bbb_p19_region_parity, tofts+patlak); --parity-suite=allmodels adds ex_tofts/tissue_uptake/2cxm reported-only. T1 map parity (test_bids_t1_map_parity_nonlinear) is also default-on. Switches (§5, deprecate-don't-break): single --parity-suite selector and --parity-thresholds JSON; old --run-*/model flags kept as aliases. CI migration tracked in PORTING_STATUS.md. Consolidation: delete 4 overlapping parity tests + dead helper cluster; trim run_dce_parity.py and CI to one region step. Pipeline feature: ROI-only fit mode (stage_overrides.fit_voxels=0) — average-then-fit per ROI, skipping the per-voxel fit, matching MATLAB's method. Fast (whole-brain xls parity 8min -> ~3s) and less noise-biased for nonlinear models. Powers test_bbb_p19_roi_xls_parity. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The ROI-only mode made test_bbb_p19_roi_xls_parity fast, so it's now default-on like the other parity tests — leaving --run-parity with no consumers. Remove the --run-parity/--parity option + run_parity fixture, and the orphaned --run-full-parity/--full-parity (its only test was the deleted full-volume parity). Drop the stale --run-parity from the CI T1 step and update the READMEs + PORTING_STATUS. Also removes the parity-testing improvement plan doc (work complete). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Fit CIs (dce_models.py): replace placeholder "CI = point estimate" returns with genuine Jacobian-based 95% intervals (beta +/- t*sqrt(diag(MSE*inv(J'J))), the confint/nlparci equivalent) for tofts, ex_tofts, patlak, vp, tissue_uptake, fxr, and 2cxm. Derived params follow MATLAB propagation (tissue_uptake vp via the Tp CI; 2cxm Ktrans=E*Fp via the delta method). Matches analytic OLS to ~1e-8. Coverage metric (phantom_gt_helpers.py, run_phantom_gt_reliability.py): report per-region ground-truth-in-CI coverage (fraction of voxels where GT falls inside the fit's 95% CI; well-calibrated ~0.95) plus standardized error z=|GT-fit|/CI_halfwidth. This is the correct accuracy-under-noise signal and is scale-free where %GT blows up near zero GT. Full CPU sweep shows ex_tofts brain Ktrans is calibrated (~0.90-0.98) while tofts/patlak are systematically biased; sub-08 (perfect T1) still fails tofts, confirming model mismatch over T1 quality. Data docs: add BIDS participants.tsv/json, dataset_description.json, and a top-level README documenting sub-01..sub-09 (real-derived reference/variants vs synthetic GT phantoms), replacing implicit folder-name-only provenance. Status doc updated with the coverage findings and next steps. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Delete tests/data/ci_fixtures and fold its data into tests/data/BIDS_test as proper BIDS subjects, removing duplicated data: - ci downsample_x2_bids was byte-identical to sub-02downsample; keep the latter and repoint run_dce_benchmark.py at it. - bbb_p19_downsample_x3y3 -> sub-10bbbdownsample (DCE fit-parity fixture); ROIs become derivatives/.../desc-*_mask.nii, MATLAB baselines under derivatives/matlabref/. - vfa_small + tiny_settings_case -> sub-11tiny (VFA anat for T1-map parity + tiny DCE for settings sweeps). - Drop unused sub-03noisyhigh and sub-04noisylow. - Rename every subject-internal file from sub-01_* to its real subject label so the dataset is valid BIDS. Repoint consuming tests, dce_default.json/dceprep_default.json, and the generation scripts (now emit the BIDS layout). Dataset-level qualification now skips sessions lacking preprocessed inputs (the fit-only fixtures) instead of failing them. ROI-xls parity compares a canonical tissue token so BIDS mask names align with the frozen MATLAB reference. MATLAB DCE parity generator gains optional explicit-path overrides for the BIDS layout. Docs updated. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Align documentation with the parity-testing-improvements branch state: - README/walkthrough: drop the removed --dce-preferences CLI example line - dce_options: note use_dce_preferences defaults false (JSON-only) with no CLI flag for the preference file - PORTING_STATUS: refresh snapshot (177 passed), repoint the renamed sub-10bbbdownsample fixture, correct the parity-alias risk/checkpoint - COMPLETED: log the parity rework, ROI-only fit mode, default-on parity, CI drift guard, JSON-only config, banner/version, plus the newly pulled Jacobian fit CIs, phantom GT coverage metric, and fixture consolidation Also ignore local RUNNER_DATA/ scratch output. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Rework the OSIPI DCE reliability suite after finding the old "peer max" gate was near-circular: the imported peer pool includes the LEK/Edinburgh implementation that ROCKETSHIP ports (2cxm, tissue_uptake), so our fit reproduces LEK and our error tracked the peer maximum to ~4 significant figures. Gating (dual gate): - Hard-gate DCE reliability + fast-backend tests on OSIPI's official published acceptance tolerances (new osipi_official_tolerances.json, transcribed from OSIPI DCEmodels_data.py; new shared helper osipi_official_tolerances.py). - Demote the peer-error spread to a non-gating signal (summary + run_osipi_reliability). - Promote the 2cxm/2cum reliability + fast CPUfit checks from xfail to real passes (they were fixed but the markers were stale). Provenance (verified against the OSIPI testing-framework paper, van Houdt et al., MRM 2023, doi:10.1002/mrm.29826): - DRO ground truth is byte-identical to OSIPI @ 23d3714 (Manning et al., MRM 2021). - osipi_peer_error_summary.json is now reproducible in-repo: commit the upstream per-contributor result CSVs (dce_models_results/, dsc_models_results/) alongside the existing T1/SI2Conc dirs, and add generate_peer_error_summary.py, which recomputes the JSON from them to machine precision (--check verifies). Summary generator (generate_osipi_summary.py): fix the path bug and the accelerated-vs-validated fit mismatch; emit provenance, a dual-gate accuracy table, per-case ground-truth-vs-fit tables, and refreshed figures. Docs (README, peer_accuracy_summary.md, PORTING_STATUS) corrected to cite the source paper and note the OSIPI tolerances are deliberately wide validity checks, not accuracy bars. Tests: 181 passed, 9 skipped, 2 xfailed (--osipi-slow). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…M fixed) Verify the three non-MATLAB fitting backends against the OSIPI DROs and make the accelerated multi-compartment fits robust. Multi-backend verification (python / cpufit / gpufit): - Label the existing reliability fits as the python backend; add full-sweep cpufit/gpufit tests gated on OSIPI's official acceptance tolerances. - Per-backend accuracy report (osipi_summary.md) + figures; README reframed for ROCKETSHIP users (why the tests exist, what they show, how to reproduce). Patched accelerator + multi-start: - Rebuilt pyCpufit 1.4.1 (Gpufit dev 3db5b4d "Fix false CONVERGED on rejected step" + 607f127 global convergence) verified in use. - Backend-agnostic multi-start (dce_pipeline._accel_multistart_refine): re-fit only voxels that pin vp/Fp to a bound (or fail) from a few perturbed starts and keep the lowest chi-square, so it never degrades a good fit and applies identically to cpufit/gpufit. Resolves tissue_uptake (2CUM) on the OSIPI sweep (1 failing case -> 0); the cpufit/gpufit 2CUM sweep is promoted xfail -> passing. 2CXM remains precision/parameterization-limited on the float32 accelerated path (vp<->Fp degeneracy; ~6/24 cases). Initialization is not the lever (warm-start from linear Patlak makes it worse); the float64 python backend, which fits E=Ktrans/Fp, is the reference. Tracked in TODO.md and the osipi-verification reference doc. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Replace the suspect-gated vp/Fp-scaling multi-start with the random log-uniform coarse->refine strategy from the Gpufit bug harness (bug/experiments.py), which is more general (escapes arbitrary wrong basins, not just bound-pinned ones): - _accel_multistart_refine: per voxel, try the fixed start + 8 log-uniform draws within bounds, cheap coarse fit (30 iters) to pick the basin, one full refine (200 iters), keep the lowest chi-square. Never degrades a good base fit; varies only initial values so cpufit and gpufit share it. Config via prefs (accel_multistart[_starts/_coarse_iters/_seed]). - 2CUM still passes the OSIPI gate; 2CXM still misses the same ~6 low-flow (Fp=5) cases (matches the harness's 18/24) where vp is weakly identifiable. Correct the root-cause framing across doc/TODO/summary/tests: the residual 2CXM failures are the caller's Fp initial guess landing in a wrong basin (Fp_init 0.35 is 50-84x the true internal Fp; any Fp_init <= 0.05 recovers the good minimum), NOT float32 precision -- a DOUBLE_PRECISION cpufit build shows the same degenerate minima, so a double-precision build is not a fix. Warm-starting Ktrans+vp from a linear Patlak fit while leaving Fp high makes 2CXM worse; the lever is Fp. TODO: the planned real fix is to reparameterize the compiled 2CXM (and 2CUM) model to fit E=Ktrans/Fp in the Gpufit CPU/CUDA fork, mirroring the float64 python backend (added as the next-up item; not started). Regenerated osipi_summary.md + figures; updated the osipi-verification reference doc, TODO.md, and COMPLETED.md. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…Jac + O(N)) Switch the accelerated 2cxm/tissue_uptake fits to the extraction fraction E=Ktrans/Fp in the compiled Gpufit fork (Gpufit dev 3bf5734): removes the Ktrans=Fp pole and its corrupted numerical Jacobian, and replaces the O(N^2) convolution with an O(N) exponential recurrence + analytic Jacobians. cpufit is verified; the gpufit CUDA kernels carry the same math (host-verified) but await CUDA-hardware confirmation. Caller (python/dce_pipeline.py): - Map Ktrans/Fp prefs -> E init/bounds (_extraction_fraction_init_bounds, mirroring the float64 python reference _fit_2cxm_osipi_canonical); recover Ktrans=E*Fp on output. - Lower the Fp floor 1e-3 -> 1e-4/s (2cxm, tissue_uptake) so low-flow tissue (fp=5 mL/100mL/min ~= 8.3e-4/s) is representable -- the missing piece that lets the low-flow cases pass. Result: all 5 accelerated OSIPI cpufit sweeps pass, including all 24 2CXM cases (was ~6/24 xfail), ~3000x faster on 2cxm (12325 -> 4.1 ms/row). Rebuilt pyCpufit into .venv (dylib md5 4a56ad4f -> 0044e3df). Tests/docs: - Un-xfail cpufit 2CXM; gpufit 2CXM stays xfail(strict=False) pending hardware. - Add test_reparam_jacobian.py (analytic Jacobian vs central differences, L2<1e-6). - Remove the --osipi-slow gate so the now-fast 2CXM/2CUM sweeps + reliability fits run by default. - Update STATUS.md (folds in the deleted gpufit_2cxm_2cum_divergence.md), TODO.md, COMPLETED.md, PORTING_STATUS.md. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
run_DCE.yml still referenced the removed tests/data/ci_fixtures tree (bbb_p19_downsample_x3y3, vfa_small) and the pre-rename sub-01 VFA filename; repoint at the BIDS_test sub-10bbbdownsample/sub-11tiny fixtures and use generate_dce_tofts_parity_map's explicit BIDS path overrides instead of its flat subjectRoot layout. Also fix a latent bug in T1mapping_fit.m surfaced by the rename: its VFA-json discovery regex assumed purely numeric subject/session labels (sub-\d+), so it silently found zero files (and crashed on the undefined `tr`) for any real BIDS label like sub-01original. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Regenerated tofts+patlak reference maps with MATLAB's current find_end_ss auto-detect (steadyStateTime=-2, unchanged default) and force_cpu=1 to match GitHub Actions' GPU-less runners. The committed maps were last regenerated ~Feb 2026 and had drifted from what CI's live MATLAB run now produces, which is why parity_checks passed locally but failed in CI. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…t to piecewise_constant Adds a documented mechanism for fixed/predictable steady-state windows: a SteadyStateEndTimeIndex field in the AIF file's JSON sidecar (same discovery convention as the existing DCE metadata sidecar), checked between the explicit steady_state_end override and auto-detection. This replaces ad hoc stage_overrides pins as the recommended way to pin a value. Also fixes Python's implicit default auto-detect method from legacy_sobel to piecewise_constant. legacy_sobel ports MATLAB's dce_auto_aif.m (a different, AIF-ROI heuristic); piecewise_constant ports find_end_ss, which is what A_make_R1maps_func actually calls by default (steadyStateTime=-2, matching script_preferences.txt's injection_time=-2 default). Python's default never matched MATLAB's real default until now. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…y state _make_config backed test_bbb_p19_region_parity/roi_xls_parity with a hardcoded steady_state_start=1/end=2 window. CI's MATLAB baseline is regenerated live with MATLAB's real auto-detect (find_end_ss, steadyStateTime=-2), which no longer resolves to that window, causing a Python-vs-MATLAB numeric mismatch that looked like a correctness bug but was really a stale, mismatched window (tofts corr=0.878, patlak corr=-0.0075 in brain). Drop the pin and adopt the steady_state_auto_method=piecewise_constant + auto_find_injection=1 combo already proven for tofts in _make_tofts_post_8ef4988_config (added in dbab95b), extending it to the full region-parity and ROI-xls tests (tofts + patlak). Verified locally against a freshly regenerated MATLAB baseline: cpu-backend Ktrans corr now 0.98+ (tofts) and 0.999997 (patlak) across brain/gm/wm. _make_tofts_post_8ef4988_config simplifies to just its one remaining unique override (aif_biexp_timing_method) since the rest is now the shared default. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Zero callers (only _resolve_baseline_window is wired into the real pipeline); it still hardcoded the old [1,2]-style fallback, which would mislead anyone who edited it thinking it was live. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
check_baseline_drift.py only guards the synthetic-curve contract baseline (matlab_reference_v1.json) and never exercises A_make_R1maps_func/find_end_ss. The NIfTI-based sub-10bbbdownsample fixture (derivatives/matlabref/...) had no analogous guard: CI regenerated it fresh every run but never diffed the fresh output against the committed copy, so the committed maps could silently drift out from under test_bbb_p19_region_parity indefinitely -- which is exactly what happened (steady-state auto-detect moved between when the maps were last committed and when CI ran, and the only symptom was a confusing Python-vs- MATLAB numeric mismatch). check_matlabref_map_drift.py regenerates the same maps into a candidate dir and compares them against the committed reference (correlation + max abs diff, tight tolerance since it's MATLAB vs. itself). Wired into parity_checks before the existing generation step, so a real MATLAB-side drift now fails loudly and specifically instead of masquerading as a Python correctness bug. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…hanism Both hardcoded steady_state_start/end in stage_overrides. Replace with a committed SteadyStateEndTimeIndex sidecar next to each fixture's AIF file (same value as before, so numeric behavior is unchanged) -- exercising the same production code path a real user pinning a value would use, instead of a stage_overrides-only backdoor. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
_run_dce_for_phantom_session pinned steady_state_start/end via stage_overrides to align Stage-A's baseline with the phantom generator's known GT baseline_images. Copy the AIF file into the scratch output_dir and write a SteadyStateEndTimeIndex sidecar next to the copy instead, so the fixed value is set the same documented way a real user would pin one -- without ever writing into the committed fixture tree. Verified end-to-end against sub-05phantom: baseGT=4 baseFit=4 (exact match, same as the prior stage_overrides pin). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
CI's first run of this new guard false-failed: nonlinear voxel CI-bound estimates (ktrans/ve/vp confidence intervals) invert a near-singular Jacobian in a handful of poorly-conditioned voxels, and are known to be sensitive to MATLAB release/OS/parfor worker count even with the algorithm unchanged (dce_preferences.txt's tight voxel_MaxIter=voxel_MaxFunEvals=50 budget makes this worse). One such voxel blew max_abs_diff up to 2.1e10 on CI's R2022a runner vs. this machine's MATLAB, while correlation over the other ~7000 voxels barely moved (0.9744 worst case, all others >=0.999997). Drop the max-abs-diff gate (kept as a reported diagnostic only) and gate purely on correlation with a generous-but-meaningful floor (0.9 -- comfortably above the observed worst-case noise, far above the ~0/negative correlation a genuine algorithm regression produces, like the steady-state window bug this guard exists to catch). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
- Introduced shared machinery in `dce_fit_backends.py` for Stage-D fitting, allowing both CPU and accelerated backends to utilize a unified candidate assembly and multi-start optimization process. - Updated `model_patlak_fit` in `dce_models.py` to leverage the new backend, simplifying the fitting process and ensuring consistent behavior across platforms. - Enhanced `dce_pipeline.py` to route Patlak model fitting through the new backend, removing redundant code and improving maintainability. - Added tests to verify that the Patlak model now uses per-voxel linear seeding for initial parameters, aligning with CPU behavior and improving fitting accuracy. - Documented the motivation and goals for the Stage-D fit backend consolidation in a new markdown file.
- Updated `model_tissue_uptake_fit` and `model_2cxm_fit` to leverage the new Stage-D fitting functions from `dce_fit_backends`, simplifying the fitting process and improving maintainability. - Removed redundant code for parameter initialization and bounds handling, centralizing logic in the backend. - Enhanced `dce_pipeline` to call the new fitting functions directly, streamlining the model fitting workflow. - Updated tests to reflect changes in fitting logic and ensure compatibility with the new backend structure.
…r handling - Migrate individual model fitting functions to a shared batch processing architecture. - Introduce `_fit_stage_d_batch` to handle fitting for multiple models, reducing code duplication. - Update `_run_scipy_per_voxel` to catch exceptions per voxel, ensuring that one bad voxel does not affect others. - Remove deprecated `run_backend_fit` function and streamline backend selection in fitting functions. - Enhance error handling and logging for backend fitting attempts, providing clearer feedback on failures. - Update tests to reflect changes in fitting logic and ensure proper backend calls are made.
- tissue_uptake/2cxm python runners double-scaled hardcoded canonical bounds by rate_in_to_min because the merged settings dict (defaults + overrides) was passed as `prefs` to _merge_prefs_in_canonical_units with an empty defaults dict. Thread the caller's raw (pre-merge) overrides through FitInputs.raw_prefs so only genuine overrides get scaled, matching the pre-refactor behavior. - 2cxm's OSIPI-canonical CI computation used the densely-interpolated fit grid's point count as dof instead of the number of actually acquired time points, understating CI width. Rescale curve_fit's covariance to the real dof. - Restore _reject_algorithm_override enforcement, now unreachable from the real Stage-D dispatch path after the refactor moved off the dce_models wrapper functions. - Remove unreachable dead code after fit_2cxm_stage_d's real return. - Fix a TODO.md checkbox marked done whose own text says the work (CUDA hardware verification) hasn't happened yet.
The Pearson correlation this check used is a sum-of-products statistic, not robust to a single extreme outlier: one near-singular tofts ve confidence-interval voxel (a known non-identifiability edge case) landed on the numerically unstable side of MATLAB's confint() on GitHub's CI runner, blowing up to ~1e10 while every other one of ~6658 voxels matched almost exactly. That single voxel collapsed Pearson corr from ~1.0 to 0.21, failing the guard despite no real drift. Verified locally: MATLAB itself is deterministic (two local runs are byte-identical and match the committed baseline exactly); the CI runner consistently computes a different value for that one voxel, twice in a row. Spearman correlation stays ~0.999 in the presence of that outlier while still collapsing to ~0 for genuine algorithm drift (verified against real reference data and synthetic uncorrelated data).
Member
Author
|
Ready to merge |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Test plan
🤖 Generated with Claude Code