Explore smooth Poincare-map proxy for alpha-loss optimization

[krystophny]

Stage: orbit classification / optimizer proxy
Source language: Fortran with optional Python plotting or validation helpers
Manuals to read first: CLAUDE.md, DOC/coordinates-and-fields.md, src/classification.f90, src/orbit_symplectic_base.f90, src/orbit_symplectic.f90, src/field/field_base.f90
Related: #39

Goal

Explore a smooth finite-time Poincare or bounce-map proxy for alpha-loss optimization. The proxy should give optimizers a differentiable or cheaply directionally differentiable target before raw long-time SIMPLE loss runs are used for validation.

Motivation

The hard SIMPLE outputs are not good local derivative targets:

• confined versus lost,
• trapped versus passing,
• regular versus chaotic,
• prompt versus non-prompt.

These labels change at class boundaries. A fixed-branch adjoint of a hard label misses the boundary motion, which is exactly what optimization needs to see. The useful first object is a smooth class-risk margin based on a finite-time map.

Exploratory design

Add a value-only prototype mode first. It should trace a fixed number of bounce tips or section crossings and emit a smooth proxy score, without changing existing hard-class outputs.

Candidate map:

z_j* -> y_j = (s_j, alpha_j)
y_{j+1} = P(y_j, field_coefficients)

Candidate score:

nu       = (alpha_N - alpha_0) / (2 pi N)
D_nu     = drift between first-half and second-half rotation numbers
Delta_s  = smooth radial width over the event sequence
d_pq     = distance to low-order p/q resonance
S_pq     = Delta_s / (d_pq + d0)
C        = sigmoid(smoothmax_pq(S_pq - 1) / eps)

For trapped-particle studies, combine C with a smooth trapped weight and use hard loss only as the validation metric.

Files to inspect

• src/classification.f90: existing tip, crossing, class, and diagnostic logic.
• src/orbit_symplectic*.f90: event tracing and implicit-step residual/Jacobian structure.
• src/field/: available field values and derivatives.
• python/ or test/python/: optional plotting or regression helpers.

First deliverable

File a follow-up implementation issue or PR that answers these questions from the current code:

• Which existing routine already records enough event data for (s, alpha)?
• Which event should define the first prototype: bounce tip or Poincare section crossing?
• Which angle definition is least ambiguous for VMEC-backed SIMPLE runs?
• Which fixed small case should be used for a golden/proxy regression?
• What diagnostics should be emitted for nu, D_nu, Delta_s, nearest resonance, smooth overlap score, and hard class?

Success criteria

make test-fast

and a short design note or follow-up issue identifies the smallest value-only implementation path.

Non-goals

• Do not autodifferentiate the raw long-time loss fraction.
• Do not replace existing hard classification outputs.
• Do not implement a full adjoint in the first step.
• Do not require SIMSOPT integration in this issue.

Later path

1. Value-only Poincare resonance scan.
2. SPSA or random-direction finite differences with common particles.
3. Tangent-linear directional derivative of the smooth margin.
4. Full finite-time Poincare-map adjoint with event-time correction.

[krystophny]

Hybrid FD/AD split note

This issue fits a mixed differentiability strategy.

A single solver can expose AD or analytical products for branch-stable subkernels while finite differences or zeroth-order samples remain around hard outputs and driver logic.

Pattern:

• use AD or analytical products for cheap local pieces: profiles, splines, interpolation, field kernels, residual kernels, or smooth proxy maps;
• keep finite differences or zeroth-order sampling for full executable calls, root and topology changes, orbit-class events, convergence logic, restart logic, or long-time loss labels;
• use cheap derivative products to choose reduced finite-difference directions or anisotropic random directions for the expensive target;
• report event changes as diagnostics, not as failed gradients.

Cost model: cheap proxy Jacobians shape where expensive evaluations are spent. This gives useful optimization signal before full end-to-end differentiability exists.

[krystophny]

Branch-margin refinement for the SensOpt connection:

The smooth Poincare/bounce-map proxy should expose explicit branch-margin residuals, not only a scalar proxy score. These margins can be used as protected residuals in SensOpt so that the SVD/search geometry is aware of nearby orbit-class or loss-channel boundaries.

Suggested margin residual vector:

m_branch(x, particle) = [
  m_loss,          # smooth distance/margin to lost-vs-confined boundary
  m_trapped,       # distance to trapped/passing boundary or trapped weight margin
  m_prompt,        # prompt-vs-nonprompt loss margin
  m_regular,       # regular-vs-chaotic / map-regularity margin
  m_resonance_pq,  # distance to nearest low-order p/q resonance
  m_overlap,       # smooth island/overlap risk margin
  m_radial_width   # finite-time radial excursion margin
]

Possible concrete definitions consistent with this issue:

nu       = (alpha_N - alpha_0) / (2 pi N)
D_nu     = |nu_first_half - nu_second_half|
Delta_s  = smooth radial width over the event sequence
d_pq     = min_{p/q} |nu - p/q|
S_pq     = Delta_s / (d_pq + d0)
C        = sigmoid((smoothmax_pq(S_pq) - 1) / eps)

Then expose something like:

smooth_proxy_score = weighted sum or smoothmax of margins
branch_margin_residual = [D_nu, Delta_s, d_pq, S_pq, C, trapped_weight, prompt_margin, loss_margin]
hard_validation_labels = [lost/confined, trapped/passing, prompt/nonprompt, regular/chaotic]

SensOpt usage:

r_p = [QS / geometry proxies from upstream, SIMPLE branch margins]
J_p = cheap FD / tangent-linear / future adjoint derivative of these smooth margins
hard SIMPLE loss labels = validation target or black-box objective only

Important distinction:

• The hard labels should not be silently differentiated.
• The smooth margins can be differentiated, finite-differenced, or used as value-only protected residuals.
• If the SVD should avoid a branch boundary, the corresponding smooth margin must be included in r_p; otherwise plain QS/proxy SVD has no reason to know about that boundary.

First implementation path:

1. value-only branch-margin residual output for a fixed particle ensemble;
2. regression plot comparing margins against hard labels;
3. directional finite differences of the margins with common particles;
4. later tangent-linear or adjoint of the finite-time map if the value-only proxy proves useful.

[krystophny]

Phase-aware SensOpt refinement:

The Poincare/bounce-map proxy can become the first concrete instance of the outer "region graph" layer described in itpplasma/sensopt#30.

Suggested SIMPLE-specific interpretation:

region node = fixed orbit/loss-class signature for a fixed particle ensemble
edge        = crossing of one smooth branch margin

Possible branch signature for a particle or ensemble:

(confined/lost status, prompt/nonprompt flag, trapped/passing weight/bin,
 regular/chaotic proxy bin, nearest resonance family p/q)

Do not attempt to trace a global boundary in phase space. Instead, when one margin becomes small or interesting, perform local one-dimensional boundary probes:

m_i(x + alpha d) = 0

where m_i is one of the smooth SIMPLE margins, e.g.

m_loss, m_prompt, m_trapped, m_regular, m_resonance_pq, m_overlap, m_radial_width

Useful probe directions:

d = -grad m_i / ||grad m_i||

or a SensOpt-preserving direction that crosses m_i while keeping QS/geometry and other margins approximately fixed.

Then evaluate candidates on both sides:

x_minus = x_boundary - eps d
x_plus  = x_boundary + eps d

using the same fixed/common particle ensemble. If the branch signature changes and the hard validation loss improves, create a new local region and restart local SensOpt there.

Policy for this issue:

• Smooth margins are optimizer/proxy quantities.
• Hard SIMPLE labels are validation and branch signatures, not AD targets.
• Phase changes should be explored by boundary probes and local restarts, not by pretending the hard label has a derivative.

Minimal follow-up after value-only proxy:

1. output branch signature and margin vector per particle/ensemble;
2. implement one margin line-search probe for a selected resonance or radial-width margin;
3. compare local-only SensOpt versus one boundary-probe restart on a tiny fixed case.

[krystophny]

Detailed implementation refinement after inspecting src/classification.f90, src/check_orbit_type.f90, and src/orbit_symplectic.f90:

SIMPLE already has the raw ingredients for a smooth proxy layer, but they are currently mostly used to produce hard labels.

Current hard classifiers / event data

classification.f90 exposes:

classification_result_t:
  passing        # trapped/passing initial class
  lost           # integration failure / loss
  fractal        # regular/chaotic footprint classifier
  jpar           # regular/stochastic by J_parallel-like invariant
  topology       # ideal/non-ideal by recurrence ordering

The code determines the initial trapped/passing status using

passing = z(5)**2 > 1 - bmod / bmax

and also writes the trapped-passing boundary in classification plots as

sqrt(1 - bmod / bmax)

This gives an explicit trapped/passing margin:

m_tp = z(5)**2 - (1 - bmod / bmax)

The code records event sequences:

zpoipl_tip   # bounce-tip / lambda-zero event footprints
zpoipl_per   # periodic boundary / Poincare-like footprints
par_inv      # accumulated parallel invariant-like quantity

and later classifies regular/chaotic by fractal dimension of these footprints. It also calls check_orbit_type, which detects recurrences, recurrence ordering, and J_parallel variation.

Smooth margin outputs to add before AD

Add a value-only diagnostic mode that returns continuous or nearly continuous margins, before attempting tangent-linear or adjoint derivatives:

m_tp                  = z5^2 - (1 - B/Bmax)
loss_margin           = min_t distance_to_loss_condition, or signed distance to radial/domain boundary
prompt_margin         = time_to_loss - prompt_time_cut, if loss occurs; positive censoring margin if not
radial_width          = smooth width of s(t) or r(t)
footprint_width_tip   = smooth diameter/variance of zpoipl_tip
footprint_width_per   = smooth diameter/variance of zpoipl_per
freq_tip              = finite-time rotation/frequency from tip sequence
freq_per              = finite-time rotation/frequency from Poincare/period sequence
freq_drift            = |freq_first_half - freq_second_half|
resonance_distance    = min_{p/q low order} |freq - p/q|
overlap_score         = smooth radial_width / (resonance_distance + d0)
jpar_variation        = max/min/variance of par_inv over recurrence groups
recurrence_margin     = signed distance in the recurrence checks used by check_orbit_type

The hard labels remain validation outputs. The smooth margins are optimizer/UQ quantities.

FD / AD policy inside SIMPLE

Use the derivative-source taxonomy from sensopt#34:

1. Analytical / existing local derivatives

   • use existing field derivatives and implicit-step Jacobians where available;
   • orbit_symplectic.f90 already contains Jacobian routines for implicit symplectic steps, e.g. jac_sympl_euler1, jac_sympl_euler2.

2. Forward/tangent-linear AD later

   • good for directional derivatives of finite-time margins along a small number of SensOpt directions;
   • propagate event sequence only while event identities remain fixed.

3. Reverse AD later

   • only for scalar smooth proxy sums over many particles, e.g.

   J_smooth = mean_i sigmoid(-m_loss_i/tau) + weights * overlap_score_i
   

   • do not reverse-AD the hard labels lost, passing, fractal, jpar, or topology.

4. Branch-preserving FD

   • valid if the same hard labels and event sequence are unchanged.

5. Cross-boundary FD / event probe

   • if labels change, classify the boundary as smooth/kink/jump/ensemble-flux instead of using the FD quotient as a derivative.

Per-particle separatrix handling

For fixed particle i, treat each margin as

m_i(p) = 0

and for a direction d estimate crossing distance

alpha_i = -m_i(p) / (grad m_i dot d)

For a finite ensemble, the hard loss/classification fraction is a step function. The useful finite-step estimator is a crossing count:

Delta L(d, eps) ≈ (1/Np) sum_{0 < alpha_i < eps} Delta label_i

This should be compared against the actual fixed-ensemble hard loss fraction.

Frequency-map / resonance-overlap bridge to NEO-RT

SIMPLE can expose a finite-time frequency-map layer analogous to NEO-RT's canonical-frequency layer:

NEO-RT:  Omega_theta, Omega_phi, R = m_theta Omega_theta + m_phi Omega_phi
SIMPLE:  finite-time frequencies from zpoipl_tip/zpoipl_per sequences

For near-QS/QH/QI fields, define a good field-line/orbit angle, compute finite-time rotation numbers, and build margins:

m_pq = min_{p/q low order} |nu - p/q|
D_nu = |nu_first_half - nu_second_half|
Delta_s = robust radial width
overlap_score = smoothmax_pq(Delta_s / (m_pq + d0))

Validate these against:

fractal classifier
J_parallel stochastic classifier
hard loss outcome
trapped/passing initial class

This provides the common ground with sensopt#35: NEO-RT uses canonical resonances; SIMPLE uses finite-time frequency/resonance-overlap margins. Both feed SensOpt as branch/proxy margins.

[krystophny]

Code-grounded design pass for the smooth proxy, based on the current main (48d1742). Line references below are from a full read of the event-detection, integrator, field, and classification paths. Transferable patterns come from the rabe bootstrap-offset differentiability campaign (Albert et al., arXiv:2407.21599) and from DESC's handling of roots and variable-count features; everything is restated self-contained.

1. Where the proposed map P(y_j) already lives

Two event families exist, both detected inside trace_orbit_with_classifiers:

• Bounce tips: sign change of z(5) (= v_par/v) from minus to plus, src/classification.f90:278. Localization is not a root solve in time: a countdown itip centers a nplagr-point stencil on the event (classification.f90:282), plag_coeff builds Lagrange weights in the variable z(5) evaluated at zerolam = 0 (classification.f90:285, src/params.f90:245), and the full state plus par_inv is interpolated by var_tip = matmul(orb_sten, coef) (classification.f90:287). Footprint (s, theta) is stored in zpoipl_tip (classification.f90:311).
• Field-period crossings (the Poincare-like section): z(3) crossing integer multiples of fper (classification.f90:340-348), same Lagrange localization in the variable phi - phiper (classification.f90:351-358), footprint into zpoipl_per (classification.f90:378).

A standalone, reusable map application already exists: trace_to_cut in src/cut_detector.f90:71-212 steps the symplectic integrator until the next event and returns var_cut = (z(1:5), par_inv) plus cut_type (0 = tip, 1 = period crossing). It is the answer to the issue's first deliverable question "which existing routine already records enough event data". Today it is only wrapped into pysimple (python/CMakeLists.txt:77); the Fortran production path re-implements the same detection inline in classification.f90. The prototype should build on trace_to_cut and retire the duplication.

The map parameter ("field_coefficients") enters through spline coefficient arrays built at init (src/get_canonical_coordinates.F90:302-380 for the canonical-flux batch splines; analogous Boozer path via splint_boozer_coord, src/field/field_can_boozer.f90:112). Field evaluation is linear in these coefficients, so they are the natural differentiable parameter interface; vmec_B_scale (src/params.f90:113) is a cheap scalar parameter for fixtures.

2. Discrepancies with the issue text

1. y_j = (s_j, alpha_j): no field-line label alpha is computed anywhere. The recorded second coordinate is canonical theta (Boozer theta in BOOZER mode), wrapped mod 2pi at classification.f90:288-289 and cut_detector.f90:112-113,141-142. The modulo destroys the winding number, so the rotation number nu = (alpha_N - alpha_0)/(2 pi N) cannot be computed from the stored footprints as-is. The prototype must record unwrapped angles alongside (trivial: keep a copy before the modulo).
2. "Poincare section" in the code means phi crossing any multiple of one field period fper, not a single fixed surface. For the proxy that is fine (it is a stroboscopic map per field period) but the write-up should say "period-crossing map", and cut_in_per (classification.f90:137) only relocates start points for plotting.
3. Two stencil conventions coexist: nplagr = 4 in the classification path (params.f90:246) vs nplagr = 6, nder = 0 in cut_detector.f90:16-17. Unify before adding tangents; otherwise fixtures validate one localization and production runs another.
4. Which event for the prototype: tips for trapped particles, period crossings for passing. check_orbit_type already consumes tip data only (classification.f90:315-319), and trapped alphas dominate losses, so tips first; the same machinery covers both.

3. Differentiability hazards along the proxy path, with specs

3a. Event detection and localization. The sign-change scan (classification.f90:278, cut_detector.f90:100-104,164-168) is the combinatorial part: which steps bracket event j is an inactive integer frame, frozen per stage. The localization itself is already smooth in the orbit states: var_cut is a rational function of the stencil (Lagrange weights from plag_coeff depend smoothly on the node values), so the tangent flows through by chain rule once the stencil membership is frozen. The hazard is degeneracy: Lagrange weight denominators contain products of node differences in the event variable, which vanish when the crossing is non-transversal (tip with d(v_par)/dt ~ 0, i.e. near the trapped-passing transition; period crossing with phidot ~ 0). Spec: either regularize those denominators, or replace localization by a Newton solve for the event time with the implicit-function tangent dt* = -g/(g_dot - eps) and dz* = z_dot dt* + dz, with g = v_par for tips and g = phi - phi_sec for crossings, regularized against g_dot -> 0. Both reduce to the standard event-time sensitivity of hybrid systems; the regularized-root-tangent pattern is the same one used as a custom derivative rule in DESC's polynomial-root code and in the bootstrap-offset campaign. Emit the transversality margin |g_dot| as a diagnostic.

3b. Implicit symplectic step. Every integrator mode solves stage equations F(x*; z_n, dt, p) = 0 by Newton with data-dependent iteration counts and tolerance exits: newton1 (src/orbit_symplectic.f90:357-385), newton2 (387-443), newton_midpoint (445-496, dgesv at 481), newton_rk_gauss (614-659, stage residual f_rk_gauss at 498-531, dgesv at 643), newton_rk_lobatto (890-939). Max-iteration exhaustion is counted, not failed (count_event at 384, 442). Do not AD through these loops. Spec: hand-written tangent of the converged step via the implicit function theorem, dx* = -F_x^{-1} (F_{z_n} dz_n + F_p dp). The F_x Jacobians already exist and are exactly the Newton Jacobians: sympl_euler1_jacobian (src/orbit_symplectic_euler1.f90:29-48), jac_sympl_euler2 (orbit_symplectic.f90:209-223), jac_midpoint_part1/2 (272-353), jac_rk_gauss (537-612). Only F_{z_n} (dependence on si%pthold and si%z) and F_p (field sensitivity through dH, dpth, vpar in get_val, src/field_can.f90:217-229) are new code. This is the same "tangent of the fixed point, not of the iteration" rule as for any implicit scheme; SIMPLE's fixed-step symplectic integrators are otherwise exactly the right substrate, since the step count and quadrature are data-independent.

3c. Post-Newton field extrapolation. extrap_field = .True. (src/orbit_symplectic_base.f90:10) makes the realized map depend on the last Newton iterate: field values are Taylor-extrapolated from the final evaluation point (sympl_euler1_extrapolate_field, orbit_symplectic_euler1.f90:80-94; inline at orbit_symplectic.f90:1271-1286 and 1334-1342) before the explicit angle advance (orbit_symplectic_euler1.f90:96-103). The primal is therefore the implicit map only up to O(tol), and its dependence on iteration count is real. Spec for tangent runs: tighten rtol and either disable extrap_field or accept and document the O(tol) primal/tangent inconsistency; the IFT tangent is of the exact implicit map.

3d. Spline field evaluation. Chain: evaluate pointer (src/field_can.f90:23-66) -> evaluate_flux/boozer/meiss/albert -> splint_can_coord (src/get_canonical_coordinates.F90:501-777) or splint_boozer_coord -> libneo batch splines. Cell selection is interval_index(j) = max(0, min(n-2, int(x_norm(j)))) (libneo src/interpolate/batch_interpolate_3d.f90:1118); within a cell the evaluation is polynomial, order 3..5 enforced at spline build (get_canonical_coordinates.F90:302,335,380). Spec: cell indices are inactive integers; the tangent is the in-cell polynomial derivative plus linear sensitivity to the coefficient arrays. No smoothing needed; C^(order-1) continuity across cells is enough. Two genuine kinks sit on the axis path: r_eval = abs(r) (get_canonical_coordinates.F90:560-563) and rho = sqrt(r) with drhods = 0.5/rho (581,591), singular at r = 0. The proxy ensemble must stay off-axis; emit a floor flag if any event lands below a small r_min.

3e. Domain guards. Loss is declared by x(1) > 1 inside each timestep routine (orbit_symplectic.f90:1207-1210, 1259-1262, 1322-1325, 1386-1389, and the Lobatto block) and newton1 even aborts mid-iteration on x(1) > 1 (375), so the lost state can be an unconverged iterate. The inner radial guard resets x(1) = 0.01 discontinuously (376, 1212-1215, 1264-1267, 1327-1330, 1391-1394). Spec: keep all of this in the primal as-is; for the proxy, loss must never be read off these branches. Use a smooth radial margin from the recorded states (3f below), mask the particle when the guard fires, and flag the tangent invalid for that particle from that event on (NaN-plus-flag, not a hard stop).

3f. Hard classifiers. All four labels are step functions and stay validation-only, as the issue demands. What the code computes underneath, and the continuous margin to expose instead:

• Trapped/passing: passing = z(5)**2 > 1 - bmod/bmax (classification.f90:150). Margin m_tp = z(5)^2 - (1 - bmod/bmax) is free. Caveat: bmax = maxval(bstart) over a discrete field-line scan (src/samplers.f90:50; get_bminmax for num_surf > 1, classification.f90:147-149) is a hard max whose argmax can jump under field perturbation. Unlike the reference choices in the bootstrap-offset campaign, where an audit showed the hard max cancels exactly in all outputs, bmax does not cancel here, so it needs a smooth max (log-sum-exp) or a frozen argmax per stage. should_skip (params.f90:427-434) thresholds the continuous trap_par (classification.f90:151) against contr_pp: frame label, freeze it.
• Regular/chaotic (Minkowski): box counting in fract_dimension (cut_detector.f90:214-274, integer bins at 245-248, integer nboxes), threshold fraction > 0.2 at the single instant kt == ntcut (classification.f90:384,391,403). Integer-valued; no useful margin can be extracted from box counts. Replace as optimizer signal by the finite-time frequency drift D_nu (section 4), which measures the same thing continuously.
• J_par regular/stochastic: check_orbit_type jumps on |par_inv(k) - par_inv(2)| > tol_perpinv = 15 (src/check_orbit_type.f90:43,139-147). The continuous residual is right there: max_k |fprs(3,k) - fprs(3,2)| is a smooth function of the tip states. Expose it as m_jpar instead of the binary ijpar.
• Topology (ideal/non-ideal): recurrence detection by sign product dt1*dt2 < 0 (check_orbit_type.f90:88), early stochastic exit for return period <= ipermin = 10 (check_orbit_type.f90:4,92-99), ordering checks with signed violations theta_2 - theta_1 (127-133), integer recurrence counts in check_recurs (154-178). Hard through and through; per the issue's later comments these belong in the branch-signature/validation layer, not in the differentiable margin vector. The signed ordering distances are available if a topology margin is ever wanted.

4. Concrete smooth margin observables, each tied to existing ingredients

All computable from a fixed number N of trace_to_cut applications, per particle:

margin	definition	existing ingredient
m_jpar	drift of per-bounce J_par: slope or (max-min)/mean of var_cut(6) over tips	par_inv accumulation classification.f90:267, per-bounce reset :312, cut_detector.f90:92,114; continuous analog of ijpar (check_orbit_type.f90:143)
nu	(Theta_N - Theta_0)/(2 pi N), unwrapped angle at period crossings (passing) or precession angle per tip (trapped)	event states; needs the unwrapped-angle fix from section 2
D_nu	abs(nu_firsthalf - nu_secondhalf)	same; continuous analog of the fractal label
Delta_s	detrended dispersion of footprints: variance of zpoipl_tip/zpoipl_per radii about a linear fit in event index	classification.f90:311,378
d_pq, C	softmin distance of nu to low-order p/q; overlap score sigmoid(smoothmax_pq(Delta_s/(d_pq+d0)) - 1)	as proposed in the issue; pure postprocessing of nu and Delta_s
m_loss	1 - smoothmax_t s over the N-crossing window (log-sum-exp)	event and step states; hard r > 1 at orbit_symplectic.f90:1207 stays the validation label
m_tp	z(5)^2 - (1 - bmod/bmax), smooth-max bmax	classification.f90:150, samplers.f90:50
m_transversal	min over events of abs(g_dot)	new diagnostic from 3a

5. Frozen combinatorial frame

Per optimization stage: fixed particle ensemble (common random numbers, zstart from src/samplers.f90), fixed N map applications, fixed event family per particle (tips vs crossings, chosen by sign of m_tp at stage start), fixed stencil memberships and spline cell indices as inactive integers. Variable event counts go into fixed-size (N)-arrays with a validity mask, the same pattern DESC uses for a variable number of magnetic wells. A particle lost before N crossings is masked from event j on and contributes (mask, m_loss), never a branch. Per-step survivor counters confpart_pass/confpart_trap (classification.f90:203-209,436-442) and times_lost (classification.f90:452, src/simple_main.f90:848) are untouched validation outputs.

6. What does not transfer (limits of the playbook)

• No IFT trick fixes exp(lambda T). Tangent and adjoint norms of the long-time flow grow like exp(lambda t); for reactor alphas over slowing-down time lambda*T ~ 1e3-1e4. Regularized event tangents make each map application differentiable; they do nothing about the product of N map Jacobians. The finite-N proxy itself is the fix, which is exactly the framing of this issue.
• Indicator outputs have zero pathwise derivative a.e. confined_fraction.dat (simple_main.f90:1109) and times_lost cannot be rescued by smoothing the integrator; the ensemble crossing-count estimator from the SensOpt comment above is the right object for the loss fraction.
• Usable N is bounded by chaos, and the bound should be measured, not assumed. Rough estimate: tangent noise ~ eps_mach * exp(Lambda N) with per-map Lyapunov exponent Lambda; demanding amplification below 1e10 (keeping ~6 digits) gives N <= ~23/Lambda, i.e. N ~ 25-230 for Lambda ~ 0.1-1. This is an order-of-magnitude statement; measure Lambda per configuration from tangent-norm growth and set N per stage. Regular orbits grow tangents only polynomially, so the proxy stays clean precisely where orbits are good, which is the regime an optimizer converges into.
• Outputs near the noise floor have noise gradients. Margins assembled from near-cancelling quantities (m_jpar from long par_inv sums) must be reported as (value, gradient, validity floor).

7. Validation plan

End-to-end AD-vs-FD over many crossings is meaningless once exp(Lambda N) dominates, and cross-machine comparison is unsound at the subsystem level: in the bootstrap-offset campaign, near-cancelling intermediates differed by up to 5e3 relative between machines running identical Fortran while end outputs held ~2e-7. Therefore subsystem tangent fixtures on same-machine oracles only:

1. Single implicit step: IFT tangent vs central FD wrt vmec_B_scale and wrt one spline coefficient, per integmode (euler1 first, then midpoint, Gauss). Target ~1e-6 relative, fixed seed, extrap_field off.
2. Single event: event-state tangent across one tip and one period crossing vs FD with the stencil frame held fixed on both sides of the FD step; include a near-degenerate tip case to exercise the regularization.
3. N-crossing proxy tangent at moderate N (start N = 16) vs FD; record degradation against the measured exp(Lambda N) curve and document the usable N. Wired into make test-fast with the reactor-scale LandremanPaul2021 QA wout.nc already used as the standard small case.
4. Value-only regression (before any tangent work): scatter of m_jpar, D_nu, Delta_s, C against the hard fractal/jpar/topology/lost labels for a fixed ensemble; the margins must separate the classes or the proxy is wrong regardless of differentiability.

8. Ordered work items (PR-sized)

1. Value-only prototype (this issue's first deliverable): new proxy module driving trace_to_cut for fixed N; unify nplagr (4 vs 6, params.f90:246 vs cut_detector.f90:16-17); record unwrapped angles and per-event (s, theta, phi, par_inv, t); emit the section-4 margin vector plus existing hard labels; default outputs unchanged. Includes validation item 7.4.
2. Degeneracy-regularized event localization plus the m_transversal diagnostic.
3. Fixed-frame masked N-crossing driver with early-loss masking and m_loss.
4. F_p and F_z tangent of one implicit step (euler1), reusing the existing Newton Jacobians; fixture 7.1.
5. Event tangent (7.2), chained N-crossing proxy tangent (7.3); measure Lambda, document usable N.
6. Smooth-max bmax and m_tp; freeze should_skip as a frame label.
7. Directional derivatives of the margin vector along optimizer directions (SPSA cross-check per the hybrid FD/AD comment above); adjoint only if the tangent proves out and the parameter count demands it.

[krystophny]

Design corrections from a ground-truth study on a model banana-tip twist map (rotation number profile crossing low-order resonances, per-bounce radial kicks from residual harmonics, exact lost/confined labels by long iteration):

1. Chaos-gate the resonance-distance margin. An orbit librating inside a p/q island has its rotation number locked at the rational, so d_pq is near zero while the orbit is perfectly confined. Raw S_pq = Delta_s/(d_pq + d0) false-positives every island-locked orbit; the margin counts only when the regularity indicator says chaotic.

2. Add a diffusive time-to-edge margin. t_esc = (s_edge - s_max)^2 / D_loc, with D_loc estimated from increments of block-averaged tips (blocks of ~20 bounces remove libration), is the strongest single per-orbit feature for the delayed-loss tail in the marginal regime.

3. Revise the per-orbit acceptance gate to AUC >= 0.7. Escape times in the sticky marginal regime are heavy-tailed; from a 200-bounce window the best feature reaches AUC 0.73 (n = 768 orbits, 45 delayed losses). The 0.9 target moves to the ensemble level (loss fraction, rank over configurations), which is what an optimizer consumes.

4. Let the observation window do the bulk of the work. Most delayed losses in the model occur within a few hundred bounces; a window of ~600 bounces resolves them by observation, leaving only the heavy tail to the proxy. This matches the existing early-classification design: classify at a fixed bounce count of order 1e2-1e3, proxy ranks the remainder.

Rotation numbers should use the weighted Birkhoff average (Das-Yorke exponential bump window; superpolynomial convergence on regular orbits) or the RRE variant (Ruth, Bindel, arXiv:2403.19003) instead of the plain mean.