Refactor classification.f90: eliminate code duplication with CutDetector

[krystophny]

Refactor classification.f90 to eliminate code duplication with CutDetector

Problem Statement

The trace_orbit_with_classifiers subroutine in src/classification.f90 is a 416-line monolithic routine that reimplements cut detection logic that already exists in the CutDetector type (src/cut_detector.f90). This creates:

• Code duplication: ~200 lines of identical tip/period detection logic
• Maintenance burden: Bug fixes must be applied in two places
• Technical debt: Prevents clean Python API for classification
• Testing difficulty: Classification logic intertwined with integration and I/O

Current Architecture

trace_orbit_with_classifiers (classification.f90:48-416)

Does everything in one routine:

1. Orbit integration (like trace_orbit)
2. Banana tip detection (v_parallel sign change)
3. Periodic boundary crossing detection (toroidal cuts)
4. Stencil interpolation (Lagrange polynomials)
5. J_parallel classification (check_orbit_type)
6. Topological ideal/non-ideal classification
7. Minkowski fractal dimension computation (fract_dimension)
8. Dynamic array management (reallocating tip/period buffers)
9. File I/O (writing 10+ fort.* files)
10. Early exit logic (fast_class, class_plot, regularity detection)

Structure:

subroutine trace_orbit_with_classifiers(anorb, ipart)
  ! Lines 48-151: Setup, coordinate conversion, early exits
  ! Lines 154-207: Allocate stencil & cut storage arrays
  ! Lines 210-401: MAIN LOOP
  !   do it=2,ntimstep                    ! Outer loop
  !     do ktau=1,ntau                    ! Inner loop
  !       call orbit_timestep_sympl(...)  ! Integration
  !       
  !       ! === DUPLICATED FROM trace_to_cut ===
  !       ! Update stencil (lines 241-249)
  !       ! Detect tip (lines 252-254)
  !       ! Interpolate tip (lines 256-263)
  !       ! Detect period (lines 305-314)
  !       ! Interpolate period (lines 316-323)
  !       
  !       ! === CLASSIFICATION (unique) ===
  !       ! classify at tip (lines 288-299)
  !       ! Minkowski at ntcut (lines 348-382)
  !       
  !       ! === FILE I/O (unique) ===
  !       ! write fort.* files (lines 379-387)
  !     enddo
  !   enddo
  ! Lines 403-416: Cleanup
end subroutine

CutDetector::trace_to_cut (cut_detector.f90:70-152)

Clean, focused function:

• Integrates orbit until next tip or periodic crossing
• Returns var_cut(6) = [s, θ, φ, |v|, v_∥, J_∥]
• Returns cut_type (0=tip, 1=periodic)
• Used by examples/orbits_and_cuts.py for visualization

Key difference: Returns one cut at a time, caller manages loop.

Code Duplication Analysis

Identical Logic in Both Routines

Feature	trace_orbit_with_classifiers	CutDetector::trace_to_cut
Parameters	nplagr=6, nder=0, npl_half=3	nplagr=6, nder=0, nplagr/2=3
Stencil update	Lines 241-249	Lines 91-100
Tip detection	Lines 252-254	Lines 106-108
Tip interpolation	Lines 256-263	Lines 111-117
Period detection	Lines 305-314	Lines 126-135
Period interpolation	Lines 316-323	Lines 138-144
Parallel invariant	Lines 241, 286	Lines 91, 117
Stencil rotation	ipoi=cshift(ipoi,1)	self%ipoi=cshift(self%ipoi,1)

Total duplicated lines: ~200 (lines 241-263 + 305-343 in classification.f90)

What trace_orbit_with_classifiers Adds

Beyond cut detection, it provides:

1. Cut storage (lines 267-285, 325-343):

   ! Dynamic arrays that grow as cuts are collected
   real(dp), allocatable :: zpoipl_tip(:,:)   ! (2, nfp_tip)
   real(dp), allocatable :: zpoipl_per(:,:)   ! (2, nfp_per)
   ! Reallocate when full (13 OpenMP critical sections!)

2. J_parallel & topological classification (lines 288-299):

   fpr_in = [var_tip(1), var_tip(iangvar), var_tip(6)]
   call check_orbit_type(nturns, nfp_cot, fpr_in, ideal, ijpar, ierr_cot)
   iclass(1,ipart) = ijpar    ! 0=unclassified, 1=regular, 2=stochastic
   iclass(2,ipart) = ideal    ! 0=unclassified, 1=ideal, 2=non-ideal

3. Minkowski/fractal dimension (lines 348-382):

   if(kt == ntcut) then
     call fract_dimension(ifp_tip, zpoipl_tip(:,1:ifp_tip), fraction)
     if(fraction > 0.2d0) then
       regular = .False.
       iclass(3,ipart) = 2  ! Stochastic
     else
       iclass(3,ipart) = 1  ! Regular
     endif
   endif

4. File output (scattered):

   • fort.20000 - trapped-passing boundary (line 101)
   • fort.10000 - forced regular passing (line 146)
   • fort.10011/10012 - regular passing/trapped (lines 437, 439)
   • fort.10021/10022 - stochastic passing/trapped (lines 443, 445)
   • fort.40012/40022/40032 - J_parallel classes (lines 456-460)
   • fort.50012/50022/50032 - topological classes (lines 470-474)

5. Early exit conditions:

   • fast_class=.true. → exit after first tip (line 297)
   • class_plot=.true. + Minkowski done → exit (line 380)
   • regular=.true. → skip integration, just count (lines 211-220)

Why This Matters

Current Issues

1. Python API blocked: Tests (test_classification_api.py, test_simple_api.py, test_batch_api.py) expect high-level classify_fast() API, but implementation requires:

   • Managing global params (class_plot, fast_class, tcut, ntcut)
   • Parsing fort.* files (no data structure return)
   • Working directory changes (fort units are global)

2. Maintenance cost: Bug fix in cut detection must be applied twice:

   • Fixed in trace_to_cut for visualization
   • Separately fixed in trace_orbit_with_classifiers for classification

3. Testing: Can't test classification logic independently from:

   • Integration
   • File I/O
   • OpenMP parallelization

4. Extensibility: Adding new classification methods requires:

   • Modifying 416-line routine
   • Understanding nested loops, dynamic allocation, OpenMP critical sections

Impact on Users

• Fortran users: Works fine, just hard to maintain
• Python users: Cannot access classification features without:
  • File I/O workarounds
  • Managing working directory
  • Parsing text files

Proposed Refactoring Plan

Phase 1: Extract Cut Storage (Low Risk)

Goal: Simplify memory management, prepare for CutDetector integration.

New type:

! In classification.f90
type :: ClassificationCutStorage
  integer :: max_cuts
  integer :: n_tips, n_periods
  
  ! Fixed-size arrays (no dynamic reallocation)
  real(dp), allocatable :: tips(:,:)        ! (2, max_cuts)  - [s, θ]
  real(dp), allocatable :: periods(:,:)     ! (2, max_cuts)  - [s, θ]
  real(dp), allocatable :: tip_data(:,:)    ! (6, max_cuts)  - full var_tip
  
contains
  procedure :: init_storage
  procedure :: add_tip
  procedure :: add_period
  procedure :: is_full
end type ClassificationCutStorage

Benefits:

• Eliminates 13 OpenMP critical sections for reallocation
• Pre-allocate based on expected number of cuts
• OpenMP-friendly (no mid-loop allocation)

Changes to trace_orbit_with_classifiers:

! Replace lines 162-167, 267-285, 325-343
type(ClassificationCutStorage) :: storage
call storage%init_storage(max_cuts=1000)

! In loop, replace reallocation with:
if(.not. storage%is_full()) then
  call storage%add_tip(var_tip)
endif

Estimated effort: 2-3 hours
Risk: Low (isolated change, preserves existing logic)

---

Phase 2: Use CutDetector for Detection (Medium Risk)

Goal: Eliminate duplicated cut detection code.

Key insight: trace_to_cut integrates until next cut, but trace_orbit_with_classifiers integrates fixed timesteps. Need to bridge this:

subroutine trace_orbit_with_classifiers_v2(anorb, ipart)
  type(CutDetector) :: cutty
  type(ClassificationCutStorage) :: storage
  real(dp) :: var_cut(6)
  integer :: cut_type, ierr, n_cuts
  integer(8) :: kt, kt_start, kt_cut
  
  ! === SETUP (keep lines 48-151 mostly as-is) ===
  ! Initialize cutty instead of local variables
  call cutty%init(fper, z)
  call storage%init_storage(max_cuts=1000)
  
  ! === MAIN LOOP: Replace lines 210-401 ===
  n_cuts = 0
  kt = 0
  
  do while(n_cuts < storage%max_cuts .and. kt < ntimstep*ntau)
    kt_start = kt
    
    ! Use CutDetector instead of inline logic
    call trace_to_cut_with_counter(cutty, anorb%si, anorb%f, z, &
                                     var_cut, cut_type, kt_cut, ierr)
    kt = kt + kt_cut
    
    if(ierr /= 0) exit  ! Particle lost
    
    ! Update confined counters (need to backfill timesteps)
    do it = (kt_start/ntau)+1, (kt/ntau)
      if(passing) then
        !$omp atomic
        confpart_pass(it) = confpart_pass(it) + 1.d0
      else
        !$omp atomic
        confpart_trap(it) = confpart_trap(it) + 1.d0
      endif
    enddo
    
    n_cuts = n_cuts + 1
    
    ! === CLASSIFICATION (keep) ===
    if(cut_type == 0) then  ! Tip
      call storage%add_tip(var_cut)
      
      ! J_parallel & topological classification
      fpr_in = [var_cut(1), var_cut(iangvar), var_cut(6)]
      call check_orbit_type(nturns, nfp_cot, fpr_in, ideal, ijpar, ierr_cot)
      iclass(1,ipart) = ijpar
      iclass(2,ipart) = ideal
      
      if(fast_class) exit
      
    elseif(cut_type == 1) then  ! Periodic crossing
      call storage%add_period(var_cut)
    endif
    
    ! === MINKOWSKI AT TIME CUTOFF ===
    if(kt >= ntcut) then
      if(storage%n_tips > 0) then
        call fract_dimension(storage%n_tips, storage%tips(:,1:storage%n_tips), fraction)
        if(fraction > 0.2d0) then
          regular = .False.
          iclass(3,ipart) = 2
        else
          iclass(3,ipart) = 1
        endif
      endif
      
      if(class_plot) then
        call output_minkowsky_class(ipart, regular, passing)
        exit
      endif
      
      ! Continue but skip integration (keep regularity shortcut)
      regular = .True.
    endif
  enddo
  
  ! === CLEANUP (keep lines 403-416) ===
  ! Write files if class_plot
  if(class_plot .and. .not. passing) then
    call output_jpar_class(ipart, ijpar)
    call output_topological_class(ipart, ideal)
  endif
  
  zend(:,ipart) = z
  times_lost(ipart) = kt*dtaumin/v0
end subroutine

New helper needed:

! In cut_detector.f90
subroutine trace_to_cut_with_counter(self, si, f, z, var_cut, cut_type, n_steps, ierr)
  ! Same as trace_to_cut, but also returns n_steps taken
  integer(8), intent(out) :: n_steps
  ! ...
  n_steps = 0
  do i=1, nstep_max
    n_steps = n_steps + 1
    call tstep(...)
    ! ... rest as before
  enddo
end subroutine

Benefits:

• ✅ Eliminates ~200 lines of duplicate code
• ✅ Reuses tested cut detection from CutDetector
• ✅ Single place to fix bugs
• ✅ Easier to understand (separation of concerns)

Challenges:

• ⚠️ Need to track timesteps correctly for confpart_pass/trap arrays
• ⚠️ CutDetector has allocatable components → not thread-safe
  • Solution: Make thread-private or use !$omp critical around trace_to_cut
• ⚠️ Behavior change: currently updates confpart every ntau steps, new version would update at each cut
  • Solution: Backfill timesteps between cuts (see loop above)

Testing strategy:

1. Run existing simple.x test cases with fast_class=.true.
2. Compare iclass(:,:) arrays bit-for-bit
3. Compare fort.* file outputs
4. Check confpart_pass/trap arrays match

Estimated effort: 1-2 days
Risk: Medium (needs careful testing, OpenMP considerations)

---

Phase 3: Separate Classification Logic (High Value)

Goal: Make classification logic reusable, testable, and accessible from Python.

New module: src/orbit_classification.f90

module orbit_classification
  use params, only: dp => real64
  implicit none
  
  private
  public :: classify_from_tips, ClassificationResult
  
  type :: ClassificationResult
    integer :: ijpar      ! 0=unclassified, 1=regular, 2=stochastic
    integer :: ideal      ! 0=unclassified, 1=ideal, 2=non-ideal
    integer :: minkowski  ! 0=unclassified, 1=regular, 2=stochastic
  end type
  
contains
  
  ! Pure function: no side effects, no file I/O, thread-safe
  function classify_from_tips(tip_cuts, n_tips, nturns, tcut, dtaumin, v0) &
       result(classes)
    real(dp), intent(in) :: tip_cuts(:,:)  ! (6, n_tips) from CutDetector
    integer, intent(in) :: n_tips
    integer, intent(in) :: nturns
    real(dp), intent(in) :: tcut, dtaumin, v0
    type(ClassificationResult) :: classes
    
    integer :: nfp_cot, ierr_cot
    real(dp) :: fpr_in(3), fraction
    integer :: i, iangvar
    
    ! Initialize
    classes%ijpar = 0
    classes%ideal = 0
    classes%minkowski = 0
    
    if(n_tips == 0) return
    
    iangvar = 2  ! Use theta angle
    
    ! === J_PARALLEL & TOPOLOGICAL CLASSIFICATION ===
    ! Process each tip sequentially (like trace_orbit_with_classifiers)
    nfp_cot = 0
    do i = 1, n_tips
      fpr_in(1) = tip_cuts(1, i)  ! s
      fpr_in(2) = tip_cuts(iangvar, i)  ! theta
      fpr_in(3) = tip_cuts(6, i)  ! parallel invariant
      
      call check_orbit_type(nturns, nfp_cot, fpr_in, &
                            classes%ideal, classes%ijpar, ierr_cot)
      
      ! Early exit logic could be added here if needed
    enddo
    
    ! === MINKOWSKI FRACTAL DIMENSION ===
    if(n_tips > 0) then
      ! Extract (s, theta) for fractal dimension
      call fract_dimension(n_tips, tip_cuts(1:2, 1:n_tips), fraction)
      
      if(fraction > 0.2d0) then
        classes%minkowski = 2  ! Stochastic
      else
        classes%minkowski = 1  ! Regular
      endif
    endif
    
  end function classify_from_tips
  
  ! Helper for batch processing
  subroutine classify_batch(tip_cuts_batch, n_tips_batch, nturns, &
                            n_particles, classes_out)
    real(dp), intent(in) :: tip_cuts_batch(:,:,:)  ! (6, max_tips, n_particles)
    integer, intent(in) :: n_tips_batch(:)  ! (n_particles)
    integer, intent(in) :: nturns, n_particles
    type(ClassificationResult), intent(out) :: classes_out(:)  ! (n_particles)
    
    integer :: ipart
    
    !$omp parallel do
    do ipart = 1, n_particles
      classes_out(ipart) = classify_from_tips( &
        tip_cuts_batch(:, 1:n_tips_batch(ipart), ipart), &
        n_tips_batch(ipart), nturns, 0.0d0, 0.0d0, 0.0d0)
    enddo
    !$omp end parallel do
  end subroutine classify_batch
  
end module orbit_classification

Update trace_orbit_with_classifiers:

subroutine trace_orbit_with_classifiers_v3(anorb, ipart)
  use orbit_classification, only: classify_from_tips, ClassificationResult
  type(ClassificationResult) :: classes
  
  ! === COLLECT CUTS (using Phase 2 code) ===
  ! ... (Phase 2 loop, but without inline classification)
  
  ! === CLASSIFY FROM COLLECTED CUTS ===
  classes = classify_from_tips(storage%tip_data, storage%n_tips, &
                                nturns, tcut, dtaumin, v0)
  iclass(1,ipart) = classes%ijpar
  iclass(2,ipart) = classes%ideal
  iclass(3,ipart) = classes%minkowski
  
  ! === FILE OUTPUT (optional) ===
  if(class_plot) then
    call write_classification_files(ipart, classes, storage, passing)
  endif
  
  ! === CLEANUP ===
  ! ...
end subroutine

Benefits:

• ✅ Python-accessible: Can wrap classify_from_tips directly
• ✅ Testable: Pure function, no side effects
• ✅ Reusable: Works with any cut data source
• ✅ Composable: Can chain with other analysis
• ✅ Thread-safe: No global state, no I/O in classification

Python API enabled:

import pysimple

pysimple.init(vmec_file, trace_time=0.015)
particles = pysimple.sample_surface(n=16, surface=0.4)

# Collect cuts (Phase 2)
cuts_data = pysimple.collect_orbit_cuts(particles, max_cuts=100)

# Classify (Phase 3 - NEW!)
results = pysimple.classify_orbits(cuts_data)

# Results dict:
# {
#   'j_parallel': np.array([1, 2, 0, ...]),  # per-particle
#   'topology': np.array([1, 0, 2, ...]),    # per-particle
#   'minkowski': np.array([1, 2, 1, ...]),   # per-particle
#   'n_tips': np.array([45, 67, 32, ...])    # per-particle
# }

Estimated effort: 2-3 days
Risk: Medium (need to preserve exact classification behavior)

---

Phase 4: File I/O Separation (Optional Polish)

Goal: Decouple file writing from classification logic.

New module: src/classification_io.f90

module classification_io
  use orbit_classification, only: ClassificationResult
  implicit none
  
  private
  public :: write_classification_files, ClassificationFileWriter
  
  type :: ClassificationFileWriter
    logical :: class_plot
    integer :: units(10)  ! Fort unit numbers
    logical :: units_open(10)
  contains
    procedure :: open_files
    procedure :: close_files
    procedure :: write_particle
  end type
  
contains
  
  subroutine write_classification_files(ipart, classes, storage, passing)
    integer, intent(in) :: ipart
    type(ClassificationResult), intent(in) :: classes
    type(ClassificationCutStorage), intent(in) :: storage
    logical, intent(in) :: passing
    
    ! Write to appropriate fort.* files based on classification
    ! (Move all write() statements here from trace_orbit_with_classifiers)
    
    select case(classes%ijpar)
      case(0)
        write(iaaa_jer, *) zstart(2,ipart), zstart(5,ipart), trap_par(ipart)
      case(1)
        write(iaaa_jre, *) zstart(2,ipart), zstart(5,ipart), trap_par(ipart)
      case(2)
        write(iaaa_jst, *) zstart(2,ipart), zstart(5,ipart), trap_par(ipart)
    end select
    
    ! ... similar for topology, minkowski
  end subroutine
  
end module classification_io

Benefits:

• ✅ Clean separation: classification logic vs. output formatting
• ✅ Can disable I/O for Python API (just return data structures)
• ✅ Easier to support multiple output formats (HDF5, NetCDF, etc.)

Estimated effort: 4-6 hours
Risk: Low (pure refactoring, preserves file formats)

---

Migration Strategy

Backward Compatibility

During refactoring, maintain both versions:

! classification.f90
subroutine trace_orbit_with_classifiers(anorb, ipart)
  ! Current implementation (keep as legacy)
  ! Mark as deprecated in comments
end subroutine

subroutine trace_orbit_with_classifiers_v2(anorb, ipart)
  ! New implementation (Phase 2+)
end subroutine

Switch controlled by compile flag:

option(USE_LEGACY_CLASSIFICATION "Use old classification routine" OFF)

Testing:

! New test: test/tests/test_classification_parity.f90
! Compare outputs of both versions for same inputs

Gradual Rollout

1. Phase 1 (Week 1): Implement ClassificationCutStorage

   • Test with current trace_orbit_with_classifiers
   • Verify identical results

2. Phase 2 (Week 2-3): Implement trace_orbit_with_classifiers_v2 using CutDetector

   • Run in parallel with legacy version
   • Compare outputs, fix discrepancies

3. Phase 3 (Week 4-5): Extract orbit_classification module

   • Implement Python wrappers
   • Update tests to use new API

4. Phase 4 (Week 6): Polish

   • Extract I/O to classification_io
   • Update documentation
   • Deprecate legacy version

5. Phase 5 (Week 7+): Remove legacy after soak period

   • Delete old trace_orbit_with_classifiers
   • Clean up compile flags

Testing Requirements

Unit Tests (New)

1. test_classification_cut_storage.f90:

   • Test ClassificationCutStorage type
   • Test overflow handling
   • Test thread-safety

2. test_classify_from_tips.f90:

   • Test classify_from_tips with known inputs
   • Compare against check_orbit_type direct calls
   • Test edge cases (0 tips, 1 tip, many tips)

3. test_cut_detector_consistency.f90:

   • Verify CutDetector::trace_to_cut gives same cuts as inline version
   • Use fixed seed, identical initial conditions

Integration Tests

4. test_classification_parity.f90:
   • Run both legacy and new trace_orbit_with_classifiers
   • Compare iclass(:,:) arrays
   • Compare fort.* file contents
   • Test with multiple VMEC files, surface values

Regression Tests

5. Golden record comparison:
   • Save outputs from current version (commit hash)
   • Re-run after each phase
   • Verify bit-for-bit identical results

Python Tests

6. test_python_classification_api.py:
   • Test pysimple.collect_orbit_cuts()
   • Test pysimple.classify_orbits()
   • Compare with Fortran simple.x outputs
   • Re-enable currently skipped tests

Success Criteria

Must Have

• Phase 1-3 complete
• All existing tests pass
• New unit tests added and passing
• Regression tests show identical results (< 1e-14 difference)
• Python API can call classification
• Documentation updated

Nice to Have

• Phase 4 complete (I/O separation)
• Performance benchmarks (ensure no slowdown)
• Python tests re-enabled and passing
• Example notebooks using new API

Metrics

• Code reduction: ~200 lines eliminated (duplication)
• Module count: +2 (orbit_classification, classification_io)
• Test coverage: +6 test files
• Python API: +2 functions (collect_orbit_cuts, classify_orbits)

Open Questions

1. OpenMP threading: Should CutDetector be thread-safe (thread-private) or use critical sections?

   • Proposal: Make thread-private copy per particle in trace_parallel

2. Timestep tracking: How to maintain exact confpart_pass/trap behavior when using trace_to_cut?

   • Proposal: Add n_steps output to trace_to_cut_with_counter

3. File I/O from Python: Should Python API support writing fort.* files or only return data?

   • Proposal: Return data by default, add legacy_files=True option

4. Backward compatibility period: How long to keep legacy version?

   • Proposal: 2 release cycles (~6 months)

5. Performance impact: Will function call overhead of trace_to_cut slow down classification?

   • Need: Benchmark Phase 2 against current version
   • Mitigation: Inlining hints, profiling

Related Issues

• Restructure Python package: clean 'simple' namespace with bundled extensions #183: Clean Python API (blocked by this)
• Tests failing: test_classification_api.py, test_simple_api.py, test_batch_api.py

References

• src/classification.f90: Current implementation
• src/cut_detector.f90: Reusable cut detection
• examples/orbits_and_cuts.py: Example usage of trace_to_cut
• src/check_orbit_type_sub.f90: J_parallel classification logic

Implementation Checklist

Phase 1: Cut Storage

• Define ClassificationCutStorage type
• Implement init_storage, add_tip, add_period methods
• Replace dynamic allocation in trace_orbit_with_classifiers
• Test thread-safety
• Verify identical outputs

Phase 2: Use CutDetector

• Add trace_to_cut_with_counter to cut_detector.f90
• Implement trace_orbit_with_classifiers_v2
• Handle timestep tracking for confpart arrays
• Add OpenMP thread-safety
• Parity testing with legacy version
• Performance benchmarking

Phase 3: Extract Classification

• Create orbit_classification.f90 module
• Implement classify_from_tips function
• Implement classify_batch for parallel
• Update trace_orbit_with_classifiers_v2 to use new module
• Add Python wrappers (f90wrap)
• Python API tests
• Documentation

Phase 4: I/O Separation (Optional)

• Create classification_io.f90 module
• Extract all file writes
• Add Python legacy_files option
• Update examples

Phase 5: Cleanup

• Deprecation warnings on legacy version
• Update all examples to new API
• Remove compile flag after soak period
• Delete legacy code
• Final documentation pass

[krystophny]

Automated triage pass (Claude), grounded in the repository at d8b0e01. Verify before acting.

Still relevant, but the issue text has drifted from the code and must be re-specified before anyone implements it. Three of the four original phases changed status:

• Phase 4 (I/O separation) is already done. classification_result_t and write_classification_results exist (src/classification.f90:34-41, :604-634), added in d607ec1 (2025-10-30, "feat: refresh python API classification (feat: refresh python API classification #184)") — one day after this issue was filed.
• The "Python API blocked" motivation is void. pysimple.classify_parallel / classify_fast / classify_fractal (python/pysimple/__init__.py:873, :987, :1035) return jpar / topology / fractal as numpy arrays via params_wrapper.get_iclass_bulk (:971-972, returned at :981-983). There is no fort.* parsing and no chdir anywhere in python/pysimple/__init__.py.
• Phases 1-3 are not done, and the premise that the two code paths are numerically identical is now false — see below. That numerical divergence is the real content of this issue today.

Also stale in the body: the file is src/check_orbit_type.f90 (module check_orbit_type_sub), not src/check_orbit_type_sub.f90.

---

Purpose

Banana-tip and period-crossing detection is the Poincaré-cut machinery behind the orbit classifiers in SIMPLE: the J∥ / topology classifier (check_orbit_type) and the Minkowski box-counting classifier (fract_dimension). It is implemented twice with different numerical parameters, so the cuts simple.x classifies on and the cuts pysimple.trace_to_cut returns are not the same object. Anyone reading tip footprints out of the Python path and comparing them to class_parts.dat gets silently inconsistent data. Consolidating also removes a maintenance hazard: any fix to tip detection currently has to be made twice.

Current state

Copy A — inline, in the classifier. trace_orbit_with_classifiers (src/classification.f90:67-500, 433 lines). Stencil fill and rotation at :299-306, tip detection/interpolation at :308-368, period detection/interpolation at :370-412, Minkowski at :414-455.

Copy B — the reusable type. cut_detector_t (src/cut_detector.f90:21-30) + trace_to_cut (:71-212). Its only Fortran caller is trace_to_cut_flat (src/simple_main.f90:1010-1063, call at :1043), the f90wrap entry behind pysimple.trace_to_cut (python/pysimple/__init__.py:713-760).

Copy C (dead). test/tests/test_poiplot_classification.f90:46-64 carries a third hand-rolled stencil. It is add_executable only, never add_test (test/tests/CMakeLists.txt:390-392).

The divergences that block a naive swap:

1. Different interpolation order. params_init sets nplagr = 4, hence npl_half = 2 (src/params.f90:360-362); cut_detector hard-codes integer, parameter :: nplagr = 6 (src/cut_detector.f90:16). Cubic vs quintic Lagrange, and a different half-stencil offset for the "cut is complete" trigger (classification.f90:313 uses npl_half, cut_detector.f90:107/:170 use nplagr/2).

2. Different parallel-invariant accumulation. Classifier:

   par_inv = par_inv+z(5)**2*dtaumin ! parallel adiabatic invariant

   (src/classification.f90:298) vs cut_detector:

   self%par_inv = self%par_inv + z(5)**2 ! parallel adiabatic invariant

   (src/cut_detector.f90:92, :155). One is a Riemann sum with the time increment, the other is a bare count-weighted sum; they differ by the factor dtaumin.

   Scope of the impact today. check_orbit_type has exactly one caller, src/classification.f90:350, and it is fed from copy A. So the copy-B normalisation does not currently affect any classification — it only makes var_cut(6) from pysimple.trace_to_cut incomparable to var_tip(6) in the classifier. The hazard is on the merge: check_orbit_type tests component 3 against an absolute threshold, tol_perpinv=15.d0 (src/check_orbit_type.f90:43) in if(abs(fprs(3,k)-perpinv_beg).gt.tol_perpinv) (:140-146). dtaumin is a normalised time step, dtaumin = dtau/ntau seeded from 2.d0*pi*rbig/npoiper2 (src/params.f90:303-305 for the VMEC branch with rbig = rt0 at :299; :281-283 is the separate test-field branch with rbig = 1.0d0). It is not unity, so adopting copy B's accumulation into the classifier would rescale the J∥ criterion against a fixed tolerance. The exact magnitude and how many orbits change class must be measured, not assumed — that measurement is what the parity test in step 5 is for.

3. Different integrator entry. The classifier calls advance_symplectic_with_retry(anorb%si, anorb%f, orbit_timestep_sympl, ierr) (src/classification.f90:272-273) and refreshes z(4)/z(5) from f%mu, f%Bmod, f%vpar (:282-284). trace_to_cut calls tstep → timestep_sympl_z (src/simple.f90:293-309), which has no Newton-retry/warning-mode recovery and never updates z(4) (it sets only z(1:3) and z(5)). Switching the classifier to trace_to_cut as-is would discard the recovery work of 6eb1542 (Continue after bounded Newton maxiter without false losses #487), f60c57c (Continue numerical recoveries and report resolved confinement #494) and d8b0e01 (Restore warning-mode orbit recovery without Newton stalls #499).

4. trace_to_cut is not thread-safe. real(dp) :: phiper = 0.0d0 (src/cut_detector.f90:86) is a procedure local with an initialiser, i.e. implicitly SAVEd and therefore shared across threads and across particles. It is only benign today because trace_to_cut_flat is called serially from Python. (The module also carries a blanket save at src/cut_detector.f90:14.)

5. Dead branch. In the first loop of trace_to_cut (do i = 1, nplagr, src/cut_detector.f90:88) the guards if(i.gt.nplagr .and. …) at :107 and :135 can never fire, so no cut can ever be reported from the stencil-priming loop.

6. Latent uninitialised read in the classifier. ideal/ijpar are declared without initialisation (src/classification.f90:99) and assigned only as intent(out) arguments of check_orbit_type inside the tip branch (:350), but are read at :459-460 (output_jpar_class(ipart, ijpar), output_topological_class(ipart, ideal)). The only way to reach :459 is ierr being set to 1 at :453, i.e. the class_plot + kt == ntcut path — every genuine integration error exits earlier at :289-294. An orbit that reaches ntcut without a completed tip stencil therefore writes garbage. class_result%jpar / %topology are initialised to 0 (:112-113) — use those.

Allocation churn: eight !$omp critical regions exist purely to grow zpoipl_tip/zpoipl_per inside the hot loop (src/classification.f90:326-340, :393-407). Both buffers are subroutine-local allocatables (:88), hence already thread-private, so the criticals buy nothing but serialisation.

Reproduction / evidence

Fetch the fixture non-interactively (test/test_data/get_test_data.sh prompts for input and is not usable in a script; test/conftest.py:39-56 is the real mechanism and pins the SHA256):

cd /path/to/SIMPLE
python -c "import sys; sys.path.insert(0,'test'); from conftest import download_wout_file; download_wout_file()"
# -> test/test_data/wout.nc
#    (wout_LandremanPaul2021_QA_reactorScale_lowres_reference.nc from simsopt,
#     sha256 333ab8ce...e112e771, pinned in test/conftest.py:24)
make CONFIG=Fast                      # never fpm; run from repo root

Observed today, statically, from the sources cited above: for the same orbit, copy A interpolates the tip with a 4-point stencil and reports var_tip(6) = Σ λ² dtaumin, copy B uses a 6-point stencil and reports Σ λ². Expected after this work: one implementation, one number.

Dynamic check to add as the parity test (plan step 5): trace one trapped particle from test/golden_record/classifier_fast/simple.in (ntestpart = 32, npoiper2 = 128, isw_field_type = 2, trace_time = 1.5d-2, fast_class = .True., class_plot = .True., deterministic = .True.) through both paths and compare the first 20 (s, θ, φ, |v|, λ, J∥) cut records. Record the measured ratio in component 6 and the measured deviation in components 1-2 in the PR description; the acceptance target after consolidation is agreement to 1e-12 relative.

Golden-record baseline that must not move: make test-golden TEST=classifier_fast and TEST=classifier_combined, compared on class_parts.dat confined_fraction.dat start.dat times_lost.dat (test/golden_record/compare_golden_results.sh:152-170; fort.* are deliberately excluded because their ordering is non-deterministic under OpenMP, see the comment at :153-155). Golden tests are registered as golden_record_<case> with GOLDEN_RECORD_RTOL=1.5e-7 (test/tests/CMakeLists.txt:1083-1091), which is what TEST= regex-matches.

Behaviour (BDD)

Scenario: One cut detector, one answer
  Given test/test_data/wout.nc and the namelist of test/golden_record/classifier_fast/simple.in
  When the new test test/tests/test_cut_detector_parity.f90 traces particle 1
       through both cut_detector::trace_to_cut and the classifier's cut loop
  Then the first 20 var_cut(1:6) records agree componentwise to rtol 1e-12
   And both used nplagr taken from params (currently 4), not a module parameter

Scenario: Parallel invariant keeps one definition
  Given cut_detector_t initialised with dtau = dtaumin from params
  When trace_to_cut returns a tip record
  Then var_cut(6) equals sum(lambda**2 * dtaumin) since the previous tip
   And that is the same quantity classification.f90:298 accumulates
   And the scale fed to check_orbit_type's tol_perpinv=15.d0 test is unchanged

Scenario: Classifier golden records are stable
  Given a deterministic-FP build (make build-deterministic-nopy)
  When make test-golden TEST=classifier_fast and TEST=classifier_combined run
  Then class_parts.dat, confined_fraction.dat, start.dat and times_lost.dat
       match the reference within GOLDEN_RECORD_RTOL=1.5e-7 / GOLDEN_RECORD_ATOL=1e-12

Scenario: Cut detection is thread-safe
  Given ntestpart = 32 and OMP_NUM_THREADS = 8
  When classify_parallel traces all particles through the shared detector code
  Then iclass(:,:) is identical to the OMP_NUM_THREADS=1 run
   And no state (phiper, kper, iper, itip) leaks between particles

Scenario: Missing tip does not produce garbage classification output
  Given class_plot = .True. and an orbit that reaches ntcut with zero completed tips
  When trace_orbit_with_classifiers writes classification output
  Then the J_parallel and topology files receive class 0 (unclassified), not an
       uninitialised integer

Implementation plan

Ordered so that every step is independently mergeable and golden-record-checkable. No sibling-repo (libneo) change is required.

1. src/cut_detector.f90 — make the detector configurable and thread-safe. Replace integer, parameter :: nplagr = 6 (:16) with components self%nplagr, self%nder, self%npl_half, self%dtau on cut_detector_t (:21-30), and make orb_sten/coef/ipoi allocatable to nplagr. Extend init(self, fper, z, nplagr, nder, dtau) (:34-69), defaulting to the params values. Change self%par_inv = self%par_inv + z(5)**2 (:92, :155) to … + z(5)**2*self%dtau. Move phiper from a saved local (:86) to a self%phiper component. Delete the unreachable guards at :107 and :135 and collapse the first loop into the second by tracking a self%nfilled counter. Remove the module-level save (:14).

2. src/cut_detector.f90 — add trace_to_cut_counted(self, si, f, z, var_cut, cut_type, n_steps, ierr) returning the microstep count, and re-implement trace_to_cut as a wrapper that discards n_steps. The classifier needs n_steps to keep confpart_pass/confpart_trap bucketed per ntau exactly as at src/classification.f90:250-261, :467-473.

3. New file src/orbit_cut_storage.f90 (add to src/CMakeLists.txt next to cut_detector.f90 at line 49; also to the list in python/CMakeLists.txt at line 78). Define cut_storage_t with n_tips, n_periods, tips(6,:), periods(6,:) and init_storage / add_tip / add_period / is_full. Pre-size from nfp (src/params.f90:357) times a growth factor; grow with move_alloc, never inside !$omp critical. Store the full six-component record, not just (s, θ) as zpoipl_tip does today (src/classification.f90:342), so the J∥ series survives.

4. src/classification.f90 — replace :299-412 (stencil, tip, period, buffer growth) with cut_detector_t + cut_storage_t calls. Keep the integrator call at :262-295 as the classifier's own (it carries advance_symplectic_with_retry); feed the post-step z into a new cut_detector entry push_state(self, z, var_cut, cut_type, have_cut) so the recovery path is preserved and only the detection logic is shared. Keep the fract_dimension calls at :420/:432 fed from storage%periods(1:2,…) / storage%tips(1:2,…). Replace output_jpar_class(ipart, ijpar) / output_topological_class(ipart, ideal) at :459-460 with class_result%jpar / class_result%topology, fixing item 6 above. Delete the eight allocation !$omp critical regions at :326-340 and :393-407.

5. New file test/tests/test_cut_detector_parity.f90, registered in test/tests/CMakeLists.txt with add_test and label fast. Asserts the parity scenario above against test/test_data/wout.nc.

6. New file test/tests/test_cut_storage.f90, registered likewise: grow past the initial capacity, verify no record is lost or reordered, verify is_full semantics.

7. New file src/orbit_classification.f90 — pure classify_from_tips(tips, n_tips, nturns, iangvar) result(classification_result_t) wrapping the check_orbit_type loop (src/classification.f90:345-363) and the fract_dimension block (:414-449), with no I/O and no globals. Point trace_orbit_with_classifiers at it. Register in src/CMakeLists.txt and export through python/pysimple alongside the existing classify_parallel.

8. test/tests/CMakeLists.txt — delete the test_poiplot_classification.x target (:390-392) or port it onto cut_detector_t; it is a third unmaintained copy that no test runs.

Steps 1-2 and 5 form one PR (tier/T2, size/M). Steps 3-4 and 6 form a second. Step 7 a third. Step 8 is trivial and can ride along.

Acceptance criteria

• grep -c 'itip\|iper\|kper\|alam_prev\|phiper\|orb_sten' src/classification.f90 returns 0 — cut detection lives only in src/cut_detector.f90.
• grep -c 'parameter :: nplagr' src/cut_detector.f90 returns 0 — the stencil order comes from params.
• make CONFIG=Fast succeeds with no new warnings.
• make test passes, including the new test_cut_detector_parity and test_cut_storage.
• make test-golden TEST=classifier_fast and make test-golden TEST=classifier_combined pass unchanged on class_parts.dat, confined_fraction.dat, start.dat, times_lost.dat.
• make test-golden TEST=test_tokamak_classifier passes unchanged.
• PYTHONPATH=build/python:build:python python -m pytest test/python/test_classification_api.py test/python/test_classify_parity.py -q passes (this is the PYTHONPATH the CMake tests use, test/python/CMakeLists.txt:88). Note test_classify_parity.py is currently not registered with add_test and self-skips unless build/simple.x exists (test/python/test_classify_parity.py:12, :61-62) — register it as part of this work.
• iclass(:,:) for classifier_fast is identical at OMP_NUM_THREADS=1 and 8.
• No !$omp critical remains inside the microstep loop of trace_orbit_with_classifiers.
• var_cut(6) from pysimple.trace_to_cut and var_tip(6) in the classifier agree to rtol 1e-12 for the parity fixture.

References

• Albert, Kasilov & Kernbichler, Symplectic integration with non-canonical quadrature for guiding-center orbits in magnetic confinement devices, J. Comput. Phys. 403, 109065 (2020), doi:10.1016/j.jcp.2019.109065 — the integrator the cut detector sits on top of.
• Albert, Buchholz, Kasilov, Kernbichler & Rath, Alpha particle confinement metrics based on orbit classification in stellarators, J. Plasma Phys. 89(3), 955890301 (2023), doi:10.1017/S0022377823000351 — defines the J∥ / topology / Minkowski classifiers implemented here; e404496 (fix(classification): fast_class regression + JPP 2023 paper reproduction #323) reproduces its figures and is the behavioural contract this refactor must not break.
• Box-counting (Minkowski-Bouligand) dimension as a regular/chaotic proxy for Poincaré sections is standard; see e.g. Meiss, Symplectic maps, variational principles, and transport, Rev. Mod. Phys. 64, 795 (1992), doi:10.1103/RevModPhys.64.795.
• ASCOT5 (Varje et al., High-performance orbit-following code ASCOT5 for Monte Carlo simulations in fusion plasmas, arXiv:1908.02482) keeps orbit-following diagnostics and end conditions in dedicated modules rather than inline in the pusher — the same separation this issue asks for.
• Beyond that, the design here is an implementation detail internal to SIMPLE; no specific literature reference applies.

Scope and risk

Three PRs, roughly size/M each; do not attempt it as one. The dangerous step is 1: changing nplagr from 6 to the params value and adding the dtaumin factor will change pysimple.trace_to_cut output (copy B is the odd one out), so that behaviour change must be called out explicitly in the PR description even though the Fortran golden records should be untouched. The main correctness hazard is the interaction between the parallel-invariant scale and the absolute tol_perpinv=15.d0 test — measure it with the step-5 parity test before merging rather than reasoning about it.

Related: #442 (SPECTRE validation suite, referenced in DOC/spectre-interface-crossing.md:150, :161, :229) exercises classification and will inherit whatever this refactor lands.

The original body's Phase 4 and the "Python API blocked" section should be struck from the issue; the Phase 1-3 checklist should be replaced by the plan above.