Paper benchmark re-run: timing + distortion figures (post-merge code)

[mvdoc]

WIP / draft. Paper benchmark re-run + figure pipeline. Benchmark-only (no package behavior changes). FreeSurfer 6 results are reused, never re-run.

The package default change (projection → continuity-only) is split into its own PR: #72. This PR's benchmark uses continuity-only explicitly, so it is independent of #72.

Why the diff is large (~9.5k lines)

dev has never contained the benchmark/ directory (it lived on benchmark-autoresearch). This PR introduces the whole benchmark suite to dev for the first time — the pre-existing autoresearch tooling plus the new paper scripts below.

What this PR introduces

Timing benchmark

• time_cores.py — flatten-only core-scaling (1/8/16 physical cores via CPU affinity; run-local truegeo + k-ring caches). --reproject re-projects continuity-only into <run>/patches.
• time_e2e.py — end-to-end (FS-free continuity-only projection + prep + flatten) @16 cores, both hemispheres per subject.

Distortion benchmark (true-geodesic, Ju 2005)

• truedist.py — the quality yardstick: libigl heat-method geodesics from sampled sources, scored as relative error |d2d − d_geo|/d_geo at each map's distance-optimal scale, reported local (≤30 mm) and global (all pairs). The k-ring energy metric is kept only as a diagnostic.
• group_flatten.py — distortion-at-scale runner: 50 subjects / 100 hemispheres from Narratives, 4 hemispheres concurrent × 8 CPUs each, resumable; both configs.
• fs6_recompare.py — re-verifies the FreeSurfer-6-vs-pyflatten comparison by scoring FS's {hemi}.flat and a fresh pyflatten flat on the same patch against the same reference.
• rebuild_truegeo.py / recompute_truedist.py — rebuild references / rewrite metric columns from saved flats without re-flattening.

Figures

• fig_paper.py — runtime, core-scaling, distortion, speed/accuracy, and per-hemisphere paired FS6-vs-pyflatten panels.
• fig_group_raincloud.py — paired LH/RH and config-paired group-distortion rainclouds (median-based).
• fig_flatmaps.py — per-hemisphere 3-method Fischl-1999 flatmap + histogram.
• fig_pipeline.py / fig_pipeline_detail.py — pipeline schematics.

Critical metric correction (this branch)

The true-geodesic reference was originally computed on the inflated surface (what patch.3d stores), where the heat solver ill-conditions and emits spurious near-zero/negative geodesics — blowing the relative metric to 200–350% on a handful of hemispheres. Fixed to compute on the fiducial surface (the surface the energy preserves), with a geodesic ≥ ½·chord sanitizer for residual solver failures. All prior true-geodesic numbers were regenerated on the corrected reference. The FS-vs-pyflatten verdict is unchanged by the fix (same reference scores both methods).

Robustness hardening (code review)

truegeo provenance stamping + fiducial-only cache reuse (a pre-fix inflated .npz can no longer be silently reused); clear errors on empty/misaligned references; fs6_recompare now delegates to the shared truedist metric instead of a divergent reimplementation (verified byte-identical); degenerate-input guard on the raincloud KDE; per-row resilience in the CSV recompute.

Method

• Two shipped variants, both Tutte flip-free init: robust_fast (default; k7/n6, line-search 7, ~3× faster) and tutte_default (k7/n12, quality-first).
• Paired hemispheres; 16 physical cores; CPU-only throughout.

Headline results (corrected fiducial metric)

• Group (100 hemispheres): robust_fast local median 14.3% / global 12.2%; tutte_default 14.2% / 12.1% — configs ~tied, LH/RH symmetric, max local ~22% (genuine).
• FS6 vs pyflatten (40 hemispheres, same patch): pyflatten better local 38/40 (14.85 vs 15.21% mean), FreeSurfer better global 38/40 (12.18 vs 12.45%) — within ~0.3 pp (the local-vs-global trade-off of Ju 2005).
• Speed: ~9.5 min (robust_fast) / ~14.8 min (tutte_default) per brain @16c vs FreeSurfer 6 ~181 min (≈12–19×); sublinear core scaling (latency-bound optimizer).
• Default choice: robust_fast Pareto-dominates (parity on distortion, ~3× faster than tutte_default).

Artifacts (not in repo)

/data2/projects/autoflatten/paper_bench_2026/ — core-scaling 20260617-123120 (+ e2e_…), group group_20260624-060144, geodesic-refined/presentation 20260614-105856; fs6_compare/true_comparison_fiducial.csv. Each run dir holds timestamped CSVs, run_meta.json, figures (PDF + 300-dpi PNG), per-hemi flatmaps, and truegeo references. See benchmark/OUTPUTS.md for the full layout, schemas, and a verified reload-and-score recipe.

TODO before un-drafting

• Continuity-only timing + group distortion runs complete; figures regenerated
• Metric corrected (fiducial) and all numbers regenerated; FS6 re-verified
• Co-author review of figures / metric choices
• Decide dev vs benchmark-autoresearch as the home for benchmark/
• Decide which large artifacts (truegeo .npz, k-ring caches) to publish vs regenerate for release

[gemini-code-assist]

Code Review

This pull request introduces significant performance and metric-fidelity improvements to the flattening pipeline, including a FreeSurfer-free orientation alignment, a distance-optimal output scale calculation, and parallelized Numba kernels for k-ring neighbor and distance computations. A performance issue was identified in _angular_kring_kernel where calling _limited_dijkstra_numba inside a parallel prange loop introduces thread-allocator contention due to internal O(n_vertices) array allocations, which should be refactored to use pre-allocated per-chunk scratch buffers.

Important

The consumer version of Gemini Code Assist on GitHub is being sunset. Starting June 18, 2026, new organization installations will be blocked, and all code review activity will officially cease on July 17, 2026.
For more details on the timeline and next steps, please review the Help Documentation.

[mvdoc]

/gemini review

[gemini-code-assist]

Code Review

This pull request parallelizes the angular k-ring kernel in autoflatten/flatten/distance.py over chunks rather than individual vertices, significantly reducing memory allocation overhead across threads. It also introduces a comprehensive benchmark and autoresearch harness under the benchmark/ directory. Feedback on the changes points out that the inlined Dijkstra algorithm in the parallel loop currently uses an inefficient linear scan to find the minimum-distance vertex, and recommends implementing a proper binary heap to achieve O(log N) complexity.

Important

The consumer version of Gemini Code Assist on GitHub is being sunset. Starting June 18, 2026, new organization installations will be blocked, and all code review activity will officially cease on July 17, 2026.
For more details on the timeline and next steps, please review the Help Documentation.