Volume classification plot: trapped-passing and deeply-trapped boundaries do not enclose data pixels

[krystophny]

Symptom

In the volume classification plots generated by
examples/classification/plot_paper_results.py:plot_volume_classification
(figures volume_classification_q[ahi].png), the two analytic boundary
curves do not enclose the per-pixel J_perp data:

• Dashed deeply-trapped boundary sits visibly below the lower edge
  of the populated trapped band. Worst case is the QI config: dashed
  line at J_perp ≈ 0.70, data lower edge at J_perp ≈ 0.85
  (~0.15 offset).
• Solid trapped-passing boundary cuts through the upper edge of
  the data rather than enclosing it from above.
• All three configurations under
  examples/classification/configs/{qi,qh,qa}_volume.in show the same
  pattern, just at different magnitudes.

The boundaries should enclose the trapped-classified pixels precisely;
they currently do not.

Likely cause: inconsistent normalization

In plot_paper_results.py, the boundaries are drawn as:

bmin_global = np.min(bminmax[:, 1])
ax.plot(bminmax[:, 0], bmin_global / bminmax[:, 1], "k", lw=1.5)     # solid
ax.plot(bminmax[:, 0], bmin_global / bminmax[:, 2], "k--", lw=1.0)   # dashed

i.e. both curves are normalized to the global minimum bmin over
all surfaces.

The per-particle J_perp plotted as image pixels comes from
class_parts.dat[:, 2], written by classification.f90:

perp_inv(ipart) = z(4)**2 * (1.d0 - z(5)**2) / bmod

i.e. v_perp^2 / B_local evaluated at each particle's starting point.

The two are not on the same scale. The boundary curves assume J_perp
is normalized so that a deeply-trapped particle at the surface bmin has
J_perp = bmin_global / bmin(s), but the per-particle quantity is
absolute v_perp^2 / B_local (or some other Fortran-side convention).
A global rescale by bmin_global is missing on one side, or the two
quantities are being computed against different references entirely.

The asymmetry across configs is consistent with this hypothesis: the
mismatch is largest where bmin(s) varies most across surfaces (QI),
smallest where the variation is mild (QH).

What to fix

Pick one canonical normalization and use it on both sides:

• Either change plot_volume_classification to draw the boundaries
  using the same reference field as perp_inv (e.g. drop the
  bmin_global factor and use a normalization derived from each
  particle's own bmod), or
• Change classification.f90 to write a normalized perp_inv that
  matches the bmin_global / bmin(s) convention used by the plot.

Either way, document the convention in a one-line comment next to
perp_inv = ... in classification.f90 and next to the boundary
plot lines so the two stay in sync.

Test

Regenerate the three volume plots from
examples/classification/configs/{qi,qh,qa}_volume.in after the fix
and verify that the dashed (deeply-trapped) and solid (trapped-passing)
curves enclose the non-NaN pixel region in all three cases.

Context

• Sampling and bminmax.dat generation for num_surf=0 were fixed
  in PR fix(classification): fast_class regression + JPP 2023 paper reproduction #323 (bugs Matlab interface working #2-4 in the PR description).
• find_bminmax cache initialization was fixed in PR find_bminmax: initialize cache in active tracing coordinates #348
  (find_bminmax: initialize cache in active tracing coordinates),
  which only affects the num_surf > 1 initialization path and
  explicitly does not touch the num_surf=0 classifier semantics.
• The boundary/data offset is a separate plotting/normalization
  inconsistency that survives both fixes.

References

• examples/classification/plot_paper_results.py — plot_volume_classification
• src/classification.f90 — line ~154, perp_inv definition
• examples/classification/configs/{qi,qh,qa}_volume.in

[krystophny]

Root cause located

Tracing back through the canonical 2022 post-processing scripts (Fortran dens_cl_glob.f90 and the matching plot_classification.py used to make the JPP 2020 / ISHW 2022 figures), the boundary curves are drawn from bminmax.dat at the physically correct locations:

• solid: bmin_global / B_max(s) (trapped-passing)
• dashed: bmin_global / B_min(s) (deeply trapped)

The bin range for the per-particle J_perp (= z(4)^2 (1 - z(5)^2) / bmod as written by classification.f90) is computed from the sample:

pmax = np.max(perp_inv)
hp   = pmax / nperp
# bin into [0, pmax], display with extent=[0,1]

For consistency with the boundary curves, this needs pmax = 1 / bmin_global (the theoretical maximum of v_perp^2 / B_local over the device, for an |p|=1 particle). The old Fortran post-processor had the same construct (pmax = maxval(perp_inv)), but with ntestpart = 500000 and uniform volume sampling the sample max sat within ~1% of the theoretical max, so the inconsistency was invisible.

With smaller particle counts (e.g. ntestpart = 100000) and/or non-uniform sampling, max(perp_inv) < 1/bmin_global. The data pixels get displayed at y = perp_inv / max(perp_inv), while the boundary curves stay at y = bmin_global / B_{min|max}(s) — so the data float above the dashed line by a factor (1/bmin_global) / max(perp_inv) - 1. Worst for configurations with large B_max/B_min variation across surfaces (QI > QA > QH in the test set).

Fix (small)

Replace

pmax = np.max(perp_inv)

with the theoretical maximum

pmax = 1.0 / np.min(bminmax[:, 1])

in plot_paper_results.py:plot_volume_classification.

Fix (preferred)

Make the Fortran side write J_perp already normalized so that J_perp in [0, 1] by construction and the python plotter just bins on [0, 1]. Add a one-liner in classification.f90:

perp_inv(ipart) = z(4)**2 * (1.d0 - z(5)**2) / bmod * bminmax_global_min

with bminmax_global_min taken from the cached bminmax(:,1) array. Then no plotter ever has to know about the convention again. Document the convention next to the assignment.

Verification

Regenerate volume_classification_{qi,qh,qa}.png from examples/classification/configs/{qi,qh,qa}_volume.in. The dashed and solid curves must enclose the non-NaN pixel region in all three plots, with the QI dashed line moving from J_perp ~ 0.70 up to ~ 0.85 (matching the data lower edge).

Provenance

• Canonical 2022 Fortran post-processor: original of examples/classification/postprocess_class.f90 (deleted in 01aedf0); see dens_cl_glob.f90 in the ISHW 2022 poster sources for the same construct.
• Canonical 2022 plotter: plot_classification.py in the same poster source tree; uses identical boundary formulas to the current plot_paper_results.py:plot_volume_classification but with the larger-sample assumption built in implicitly.