fix(classification): fast_class regression + JPP 2023 paper reproduction#323
Merged
Conversation
ⓘ You are approaching your monthly quota for Qodo. Upgrade your plan Review Summary by QodoRestore volume sampling and bminmax output for classification
WalkthroughsDescription• Restore volume sampling for num_surf=0 to enable classification diagnostics • Fix fast_class and tcut compatibility by removing incorrect mutual exclusivity check • Generate bminmax.dat output file with 101 rows for classification plotting • Migrate classification binning from Fortran postprocessor to Python script • Fix hardcoded array index and VMEC file reference in classification example Diagramflowchart LR
A["startmode=1<br/>num_surf=0"] -->|"volume sampling<br/>uniform s in [0,1]"| B["sample_particles"]
B --> C["compute_pitch_angle_params"]
C -->|"get_bminmax<br/>for all num_surf"| D["bmin/bmax arrays"]
D --> E["write_output"]
E -->|"generates"| F["bminmax.dat<br/>101 rows"]
G["class_parts.dat<br/>6 columns"] --> H["plot_classification.py"]
F --> H
H -->|"bin in Python"| I["class_jpar.pdf<br/>class_ideal.pdf"]
File Changes1. examples/classification/postprocess_class.f90
|
Code Review by Qodo
1. Plot binning normalization bug
|
![qodo-code-review[bot]](https://avatars.githubusercontent.com/in/484649?s=60&v=4){kind=small}
Comment on lines
+10
to
+22
| hs = 1.0 / ns | ||
| pmin = 0.0 | ||
| pmax = np.max(perp_inv) | ||
| hp = (pmax - pmin) / nperp | ||
|
|
||
| prompt = np.zeros((nperp, ns)) | ||
| regular = np.zeros((nperp, ns)) | ||
| stochastic = np.zeros((nperp, ns)) | ||
|
|
||
| for ipart in range(len(s)): | ||
| i = min(ns, max(1, int(np.ceil(s[ipart] / hs)))) - 1 | ||
| k = min(nperp, max(1, int(np.ceil(perp_inv[ipart] / hp)))) - 1 | ||
|
|
There was a problem hiding this comment.
1. Plot binning normalization bug 🐞 Bug ✓ Correctness
bin_classification() normalizes perp_inv by its sample max, but doplot_inner() overlays bmin/bmax boundary curves in a different normalization, so the heatmap bins and boundary curves are on inconsistent y-scales and the plot can be physically misleading. If all perp_inv values are identical (e.g., all 0), hp becomes 0 and the script divides by zero when computing bin indices.
Agent Prompt
### Issue description
`examples/classification/plot_classification.py` bins `perp_inv` using `pmax = np.max(perp_inv)`, but overlays curves computed from `bminmax.dat` that assume normalization by the global minimum B (`bmin_global`). This makes the plotted heatmap and boundary curves inconsistent, and can also crash when `pmax == pmin` (hp=0).
### Issue Context
The boundary curves use `bmin_global / bminmax[:, 1|2]`, which is consistent with plotting `perp_inv * bmin_global` (dimensionless, ideally in [0,1]). The binning should use the same quantity, and must guard against `hp == 0`.
### Fix Focus Areas
- examples/classification/plot_classification.py[5-31]
- examples/classification/plot_classification.py[33-56]
### Suggested approach
- Load `bminmax` first in `main()` and compute `bmin_global`.
- Define `perp_norm = perp_inv * bmin_global` and clip to `[0, 1]` (or set explicit bounds).
- Update `bin_classification()` to bin `perp_norm` in a fixed range `[0,1]` (so `hp = 1.0/nperp`), avoiding `pmax=np.max(...)`.
- Add a defensive guard so `hp` cannot be 0 (even if you keep data-driven pmax).
ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools
Comment on lines
771
to
775
| !$omp critical | ||
| bmod = compute_bmod(z(1:3)) | ||
| if (num_surf > 1) then | ||
| if (num_surf /= 1) then | ||
| call get_bminmax(z(1), bmin, bmax) | ||
| end if |
There was a problem hiding this comment.
2. Classifier skips bminmax for num_surf=0 🐞 Bug ✓ Correctness
compute_pitch_angle_params() now calls get_bminmax when num_surf/=1, but trace_orbit_with_classifiers() still only calls get_bminmax when num_surf>1. If num_surf remains 0 (e.g., startmode=2 reading a start.dat with varying radii), classification can use stale bmin/bmax and compute incorrect passing/trap parameters.
Agent Prompt
### Issue description
`src/simple_main.f90` updated the bmin/bmax lookup guard to `num_surf /= 1`, but `src/classification.f90` still uses `num_surf > 1` before calling `get_bminmax`. This can skip per-radius bmin/bmax updates when `num_surf==0`.
### Issue Context
Classifier computations (`passing`, `trap_par`, `perp_inv`) depend on accurate `bmin/bmax`. For `num_surf=0` (volume / mixed-radius starts), bmin/bmax should be obtained via `get_bminmax`.
### Fix Focus Areas
- src/simple_main.f90[762-777]
- src/classification.f90[148-155]
### Suggested approach
- Change the guard in `trace_orbit_with_classifiers` from `if(num_surf > 1)` to `if (num_surf /= 1)` (matching simple_main), or explicitly handle `num_surf==0` as a case that must call `get_bminmax`.
ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools
krystophny
changed the title
krystophny
force-pushed
the
feature/venv-setup
branch
from
31cfad8 to
50a0fce
Compare
krystophny
force-pushed
the
fix/classification-regression
branch
from
a3210dd to
cec6270
Compare
krystophny
force-pushed
the
fix/classification-regression
branch
from
4105177 to
23db82f
Compare
krystophny
added a commit
that referenced
this pull request
May 21, 2026
The PR #323 volume configs diverged from the published paper recipe in every knob except the equilibrium file. Bring them back into line with Albert+ 2023 JPP doi:10.1017/S0022377823000351, page 11: - ntestpart 100000 -> 200000 (paper used 500000; compromise for runtime) - trace_time 1.5d-2 -> 1.5d-1 (paper-era poster setting) - tcut 1d-2 -> 1d-1 (classification cut at 100 ms per paper §3) - class_plot .True. -> .False. so particles keep their random (s,theta,phi) from num_surf=0 volume sampling instead of being pinned to phi=0. class_plot=.True. in classification.f90 forces z(3)=cut_in_per*fper for every particle, which collapses the sampling onto a single phi plane; for QI in particular the per-surface bmin lives at phi != 0, so the volume J_perp data tops out well below the deeply-trapped boundary. - add multharm = 7 (paper value; current default differs) - remove cut_in_per (no longer used with class_plot=.False.) plot_paper_results.py: nperp 100 -> 300 to match the paper's 50 x 300 (s, J_perp) grid.
krystophny
force-pushed
the
fix/classification-regression
branch
from
0641818 to
7e7e952
Compare
krystophny
changed the title
fast_class previously set ierr = ierr_cot, which propagated the classifier completion status to the integration error flag. Every classified orbit (including regular ones) was then marked lost, so confined fractions were wildly wrong (e.g. fc=0.57 instead of 0.94 for QI s=0.3 at 100k particles). Fix: when fast_class is on and class_plot is off, set regular=.True. for regular J_parallel orbits instead of stomping ierr. That stops integration early, counts the orbit in confpart for the remaining time steps, and records times_lost = trace_time. When class_plot is on, the old ntcut termination is preserved so classification output is still written. Also remove the read_config guard rejecting fast_class + tcut > 0: the JPP 2023 recipe combines fast classification with a tcut, and the guard blocks that valid combination. Add a classifier_combined golden record exercising both flags together.
Add per-configuration input files for the three stellarators studied
in Albert et al., J. Plasma Phys. 89(3), 955890301 (2023):
- QI (Drevlak 2014, wout_23_1900_fix_bdry.nc)
- QH (Drevlak 2018, wout_qh_8_7.nc)
- QA (Henneberg 2019, wout_henneberg_qa.nc)
For each, four cases reproduce the published figures:
- q?_volume: Figs 9-11, volume (s, J_perp) classification maps
(startmode=5, num_surf=0, ntestpart=200000,
multharm=7, trace_time=0.15 s, tcut=0.1 s,
fast_class=.True.)
- q?_s03: Figs 5-7, loss curves at s=0.3
- q?_s06: Figs 5-7, loss curves at s=0.6
- q?_fig8: Fig 8, classification dashboard at s=0.6
plot_paper_results.py generates fig5-7 losses, fig8 dashboard, and
volume_classification + volume_topology (s, J_perp) maps (J_par and
ideal-topology classifier columns) from class_parts.dat /
times_lost.dat / bminmax.dat / confined_fraction.dat. run_all.sh
downloads the wout files from the proxima-simple-classification repo
on first call and runs all 12 cases under /tmp/simple_classification.
Equilibrium files, plots, and committed plot-essential run outputs
live in the companion repository:
https://github.com/itpplasma/proxima-simple-classification
Remove the legacy in-place runner (plot_classification.py,
postprocess_class.f90, simple.in, README.md) superseded by the new
script set. Makefile, plot_paper_results.py and run_all.sh headers
consistently cite the 2023 paper as the reproduction target.
…vailable Wrap the numpy/matplotlib/netCDF4 imports in a try/except so the plot test exits 0 with a 'Skipping plot test: <reason>' message when the optional Python deps are not installed, instead of ImportError crashing the test job.
krystophny
force-pushed
the
fix/classification-regression
branch
from
7e7e952 to
52639bd
Compare
…r_fast The golden record compares the current branch to a freshly built `main` binary. When a bug fix on this branch legitimately changes a case's output, the comparison will always fail until the fix lands on main. Add GOLDEN_RECORD_SKIP_CASES (space-separated list of case names) to compare_golden_results.sh and consume it from CMake test registration. Use it to skip classifier_fast while the fast_class + ierr fix in classification.f90 is in flight: ref/cur divergence is the intended outcome (regular orbits now lasting trace_time instead of being marked as early losses). classifier_combined (new) and classifier_fast both exercise fast_class with class_plot=.True., so the new code path is still covered by classifier_combined and the new path agrees between ref and cur because main's read_config used to reject fast_class+tcut and is skipped on the reference side. Once this fix lands on main and main is rebuilt as the reference, remove classifier_fast from GOLDEN_RECORD_SKIP_CASES.
3 tasks
krystophny
added a commit
that referenced
this pull request
May 22, 2026
## Summary Follow-up to #323. The squash-merged classification fix means main now matches the output that this branch's `simple.x` produces for `golden_record_classifier_fast`, so the temporary `GOLDEN_RECORD_SKIP_CASES=classifier_fast` escape hatch added in #323 can be removed. Two small changes: 1. `test/tests/CMakeLists.txt`: drop `GOLDEN_RECORD_SKIP_CASES=classifier_fast` from the per-test ENVIRONMENT. The `SKIP_CASES` plumbing in `compare_golden_results.sh` stays in place for future intentional-divergence bug fixes (just don't set the env var unless a fix needs it). 2. `test/golden_record/compare_golden_results.sh`: drop `avg_inverse_t_lost.dat` from the classifier file list. That file is only written when at least one particle is actually lost; with the fast_class fix the small `classifier_fast` test loses zero particles, so neither ref nor cur writes it. The comparator's "reference file missing" check then spuriously fails even though both runs agree byte-for-byte on the files that do exist. ## Verification Full regression suite, against freshly rebuilt main reference (cached `runs/run_main` and `simple_main` rebuilt to pick up #323): ``` $ make test-regression ... 13/13 Test #50: golden_record_sanity .................... Passed 0.57 sec 100% tests passed, 0 tests failed out of 13 Total Test time (real) = 411.09 sec ``` Both classifier cases: ``` 1/1 Test #45: golden_record_classifier_fast ... Passed 16.73 sec 1/1 Test #44: golden_record_classifier_combined ... Passed 16.18 sec ``` ## Test plan - [x] `make test-regression` passes 13/13 locally - [x] `golden_record_classifier_fast` (the test #323 had to skip) passes - [x] `golden_record_classifier_combined` (new in #323) still passes
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
1 participant
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
Three logical commits that together reproduce the orbit-classification
figures of the 2023 JPP paper (Albert, Buchholz, Kasilov,
Kernbichler, Rath, J. Plasma Phys. 89(3), 955890301,
doi:10.1017/S0022377823000351)
on current
mainand fix two classification regressions blocking it.fix(classification): allow fast_class + tcut togetherfast_classno longer marks regular orbits as lost. Previouslyierr = ierr_cotpropagated the classifier completion status tothe integration error flag and every classified orbit was reported
as lost (confined fraction ~ 0.57 instead of 0.94 for QI s=0.3).
The fix sets
regular = .True.for regular J_par orbits underfast_class .and. .not. class_plot;class_plotmode keeps theold ntcut-termination so classification output still gets written.
read_configno longer rejectsfast_class .and. tcut > 0. Thepaper recipe combines both flags.
classifier_combinedgolden record exercises both together.feat(examples/classification): JPP 2023 paper reproduction recipeq{a,h,i}_{volume,s03,s06,fig8}.in) coveringFigs 5-11 of the paper.
plot_paper_results.pyproduces fig5-7 losses, fig8 dashboard,and volume_classification + volume_topology (s, J_perp) maps.
run_all.shdownloads the wout files from the companion repo andruns all 12 cases under
/tmp/simple_classification.plot_classification.py,postprocess_class.f90,simple.in, and exampleREADME.mdsuperseded by the new script set.
test(magfie): skip orbit chartmap comparison plot when matplotlib unavailableare missing; it now exits 0 with a
Skipping plot test:message.Companion repository with equilibrium files, plot-essential run
outputs, and generated figures:
https://github.com/itpplasma/proxima-simple-classification
Verification
Full
make test-fastlocally (100% pass)Failing test before the fix (
test_bminmax_driver, casebminmax_classifier_num_surf_zero) now passes; was caused by anearlier branch revision that bumped
num_surfforstartmode=1 + num_surf=0. The volume configs now usestartmode = 5(the standardvolume-sampling entry) instead, which keeps the bminmax cache
semantics unchanged from #348.
Confined fractions match paper recipe (volume runs, 200k particles)
QA, QH, QI volume
(s, J_perp)classification maps regenerated withthis branch are at:
https://github.com/itpplasma/proxima-simple-classification/tree/main/plots
Single-surface loss curves (Figs 5-7) and dashboard (Fig 8) are being
regenerated this week and will land in the same
plots/directory.Test plan
make test-fastpasses 100% (44/44) locallytest_bminmax_driver(which was failing on prior branch state)passes after switching volume configs to
startmode=5classifier_combinedgolden record coversfast_class + tcut(s, J_perp)plots reproduce the paper's QI, QH, QAmaps with both J_par and ideal-topology classifier columns
multharm=7(in flight, ~4-6 days serial)