Expose the continuous margins behind the orbit classes - #513
Conversation
The classifier forms two continuous quantities and then discards them. The J_parallel class thresholds the spread of the parallel invariant across banana tips against tol_perpinv; the ideal-orbit class is the sign of the minimum monotonicity margin over the ordered tip sequence. Only the integer codes survive into class_parts.dat. Anything that wants a smooth confinement score has to reconstruct those quantities by re-tracing, and any attempt to differentiate a metric built on the integer codes gets exactly zero, silently, because the codes are piecewise constant. Keep both quantities and write them to class_scores.dat: J_parallel spread, the reference magnitude it should be compared against, the monotonicity margin, the score status, and trap_par. class_parts.dat keeps its column contract, so existing readers and both golden records are untouched. The classifier loops no longer exit early, since a maximum and a minimum need every term. That costs a few comparisons per resolved orbit and leaves the integer classes bit-identical. score_status distinguishes the cases a consumer must not conflate: 1 means the topology margin is a real margin, 2 an early stochastic exit that forms neither quantity, 3 a resolved orbit whose ideal class came from the recurrence test, which produces no margin. The J_parallel spread is valid for 1 and 3.
The discrete barrier overlap counts integer class codes, so its gradient is zero almost everywhere and a delta on the switching set. Pointing any differentiation tool at it returns 0.0 silently. Three stacked discontinuities cause that: the trapped mask, mu bin membership, and the class code. The first two are mollified here with a logistic in trap_par and a smooth bin kernel. The third cannot be mollified after the fact, because only the code survives; it needs the margins the classifier forms and discards, which SIMPLE now writes to class_scores.dat (itpplasma/SIMPLE#513). Two continuous scores follow. The J_parallel score is a logistic in the tip spread relative to tol_perpinv, so the classifier's threshold is its zero-width limit. The topology score is a logistic in the negated monotonicity margin, falling back to the J_parallel score where the recurrence test decided and no margin exists, since a constant there would flatten the metric over a whole class of orbits. Both are finite-time over nturns bounce periods, so Lyapunov exposure stays a few bounce times rather than the 1e3-1e4 of a slowing-down trace. That is what makes this the differentiable route and not merely the smooth one. The optimizer objective is now selectable between the discrete metric and either smooth variant, and every response records the discrete value whatever is optimized, so the campaigns stay comparable. Runs without class_scores.dat report the smooth block as unavailable rather than failing, so the discrete campaign is unaffected by the SIMPLE build. Tests pin the property that makes the surrogate valid: as the three widths go to zero the smooth overlap converges to the discrete one, for both classifiers, and it moves under a sub-threshold drift change that the discrete metric cannot see at all.
The class margins are thin on real runs. On a reactor-scale candidate at nturns = 8 only 2 of 1024 orbits carried a J_parallel spread and none carried a monotonicity margin, because both need the classifier to complete nturns return periods of the tip map, which a 20 ms trace rarely reaches. A score with no data cannot steer an optimizer. The tip sequence carries a quantity that needs none of that: the radial excursion max(s) - min(s) over the banana tips collected so far. Two tips suffice, no threshold, no recurrence, no nturns. It is the width of the radial band the orbit explores, which is exactly what a phase-space barrier suppresses, so it is the natural continuous confinement measure rather than a mollified class code. The driver test now exercises it on 22 orbits where the class margins reach 4 and 0, and checks the pairing both ways: two tips imply a positive excursion, fewer than two imply none.
There was a problem hiding this comment.
Review verdict: Approve
Summary: This small PR exposes the two continuous quantities the orbit classifier computes and then discards — the J_parallel spread across banana tips and the ideal-orbit monotonicity margin — plus a reference value and a score status, writing them to a new class_scores.dat while leaving class_parts.dat (and its column contract / golden records) untouched. The classifier loops were converted from early-exit to full-scan max/min accumulation; I verified the resulting integer classes stay bit-identical (the ijpar=2/ideal=2 outcomes are still set identically), and the new module variables are threadprivate with per-particle columns written race-free (each thread owns a distinct ipart). The accompanying driver test validates that thresholding the scores reproduces the discrete classes. I could not build/run here (no NetCDF/numpy toolchain available), so this is static analysis only.
Findings:
-
[minor] test/tests/test_class_scores_driver.py:117-121 — The topology-margin path (status 1) is only "opportunistically exercised," and the test never asserts
margins > 0; the commit itself notes this field/surface "do not reliably produce" status-1 orbits. The ideal-orbit half of the headline claim ("thresholding the continuous score must reproduce the discrete class") may therefore be entirely vacuous in CI, with only the J_parallel relation load-bearing. Consider a config known to produce status-1 orbits, or assertmargins > 0with parameters/equilibrium chosen to fire it so the margin relation is actually covered. -
[minor] src/check_orbit_type.f90:159-170 — The status-3 documentation ("resolved but the ideal-orbit class came from the recurrence test, which forms no margin") is inaccurate for the
ideal == 1empty-loop case. When the monotonicity loop yields no candidate (e.g.iret(nturns)-iret(nturnm1) == 0, ornturns <= 2),idealis 1 by a vacuous monotonicity check yetscore_statusbecomes 3, which the doc ties to the recurrence test. The code is safe (status 3 just means "no real margin"), but the documented semantics of status 3 should acknowledge thatideal==1without a margin also lands there. -
[minor] src/check_orbit_type.f90:169-170 / src/classification.f90:365-366 —
jpar_ref = abs(perpinv_beg)is documented as "the reference magnitude it should be compared against," but the classifier actually thresholdsjpar_spreadagainst the signedperpinv_beg(fprs(3,2)), not its magnitude. A downstream consumer reproducing the J_parallel class from the file alone cannot reconstruct the signed reference, so the comparison is not self-contained. Either document that the basis is the signed first-tip value or emit the signed value.
Verdict: Approve — the change is small, correct, race-free, preserves bit-identical existing classes and file contracts, and ships a self-validating test; the findings are non-blocking documentation/coverage nits.
There was a problem hiding this comment.
Review verdict: Comment
Summary: This PR introduces continuous fast-classifier scores (J_parallel variation rate, rotation-number drift, achieved precession turns, resolution metadata, and legacy margins) computed per banana tip in check_orbit_type.f90, stored per particle in params::class_scores, written to a new class_scores.dat output, and adds an explicit unresolved guard when a tip interpolation yields non-finite values. It ships an analytic unit test (test_continuous_classifier_scores.f90) and an integration driver (test_class_scores_driver.py). I verified the score formulas analytically against all five unit-test expectations and cross-checked the driver thresholds (TOL_PERPINV=15.0 matches tol_perpinv=15.d0; rotation-sample counts, jpar-sample bounds, and status-1/2/3 semantics agree with the code) and found no functional defect. The OpenMP threadprivate variables and the pure subroutine are correctly scoped; the NaN guard maps to ORBIT_EXIT_NUMERICAL_EVENT (=104 >= ORBIT_EXIT_NUMERICAL_DOMAIN), so unresolved tips are correctly excluded from confinement statistics.
Findings:
- minor src/simple_main.f90:2831 — the new class_scores.dat output is not documented anywhere (CLAUDE.md output list, examples, or docs); the 12-column contract lives only in scattered code comments (params.f90:100-108, simple_main.f90, tests). Recommend documenting the format for external consumers.
- minor src/classification.f90:337-348 — the non-finite-tip guard introduces a distinct branch (iclass=0, exit_code=NUMERICAL_EVENT, ierr=1, orbit abandoned) with no dedicated integration coverage; only the pure-subroutine edge is tested. A targeted test would lock in the zero-sample-count exclusion contract.
- minor test/tests/test_class_scores_driver.py:120-144 — legacy-threshold checks (spread>15.0, margin<0.0) duplicate Fortran literals (tol_perpinv=15.d0); future tolerance changes could silently break/false-pass without a shared source of truth.
Verdict: Comment — the change is well-structured and analytically correct, but I could not execute the added tests (no build/CI evidence in this read-only review), and the minor documentation/threshold-coupling issues are worth addressing before relying on class_scores.dat as a stable output contract.
Stacked on #512 (
fix/fast-class-standalone-dispatch); review that first.Risk tier
New diagnostic output. Integer classes are unchanged.
Correctness contract
Intended behavior change
A classifying run additionally writes
class_scores.dat: the J_parallel spreadacross banana tips, the reference magnitude, the ideal-orbit monotonicity
margin, a status, and
trap_par. These are quantitiescheck_orbit_typealready forms and discards.
Behavior that must not change
The integer classes in
class_parts.dat, its column contract, and both goldenrecords. The classifier's two decision loops no longer exit early — a maximum
and a minimum need every term — but the decisions they make are identical: the
J_parallel test still sets
ijpar = 2on any term above tolerance, and themonotonicity test still sets
ideal = 2on any inversion.Coordinate / unit conventions
jpar_spreadandjpar_refare in the same units as the tip-interpolatedparallel invariant
var_tip(6), so the ratio is dimensionless and the existingtol_perpinv = 15applies to the spread unchanged.topology_marginis anangle difference in radians.
Numerical invariants
jpar_spread >= 0,jpar_ref >= 0. Thresholdingjpar_spreadattol_perpinvreproducesijpar; the sign oftopology_marginreproducesidealwheneverscore_status == 1.Tests added
test_class_scores_driver.py. The oracle is the integerclassification SIMPLE already writes — thresholding each continuous score
must reproduce the discrete class it was derived from, a property of the pair
of files rather than of either one's values. It fails rather than passing
vacuously when no orbit resolves.
Coverage is honest about its limit: the run exercises the J_parallel
comparison on 4 orbits and the topology margin on 0, because
score_status == 1needs an orbit the recurrence test calls one-line with a non-empty monotonicity
interval, which this field and surface do not reliably produce. The test prints
both counts.
Golden-record impact
class_parts.datis byte-identical;class_scores.datis a new file.Failure modes considered
a margin will read 0 for orbits decided by the recurrence test. Status 3
exists to make that explicit, and status 1 is defined as "a margin was
actually formed" rather than "the orbit resolved".
threadprivatelist indetect_oneline_mod; the classifier is called from an OpenMP loop.orbit and nothing on unresolved ones.
Manual validation
make CONFIG=Fast;make test TEST="class_scores|fast_class|bminmax"— 5/5 pass.make test-fast: the same 4 failures as cleanmainat04b31af(
test_chartmap_startmode1,test_spectre_sympl_volume,test_spectre_sympl_crossing,test_spectre_validation), pre-existing andunrelated.