Skip to content

Vectorized, modular EAS optical Cherenkov engine (cphotang refactor) - #124

Open
Areustle wants to merge 10 commits into
mainfrom
cphotang_refactor
Open

Vectorized, modular EAS optical Cherenkov engine (cphotang refactor)#124
Areustle wants to merge 10 commits into
mainfrom
cphotang_refactor

Conversation

@Areustle

@Areustle Areustle commented Jul 30, 2026 •
edited
Loading

Copy link
Copy Markdown
Contributor

Summary

Full rewrite of the EAS optical Cherenkov computation (cphotang.py) into a vectorized, modular, pure-Python engine. Closes #122 and closes #113.

The monolithic per-event CphotAng loop is decomposed along its call-graph structure — each concern now lives in a single-purpose module, and cphotang.py itself is down from a 2400-line monolith (~1600 after vectorization) to a ~950-line orchestrator:

Concern (#113 checklist) Module
Atmospheric interactions atmospheric_models.py — ozone/aerosol profile tables, exact analytic layer-walk slant columns (ozone_column, aerosol_column), refractive_index
Shower propagation propagation.py — canonical slant_depth / length_at_depth oracles plus the GL-in-slant-depth longitudinal grid (propagation_grid, visibility_window, gl_node_grid, node_geometry) as pure functions; shower_properties.py — per-node shower state (shower_state_at_nodes, shower_age_at_depth, greisen_particles, hillas_e2)
Cherenkov light generation hillas_kernel.py / hillas_batch_kernel.py — θ-collapsed single-integral Gaisser-Hillas kernel + the pluggable photon-yield factory make_hillas_photon_yield (entry point for CHASM & other engines, #115); quadrature.py — GL primitives
Multiple shower dispatch dispatch.py — map_showers_distributed (chunking, dask cluster lifecycle, progress); utils/distributed.py — BackgroundCluster warm-start helper

Vague heritage names got honest ones (tracklen → secondary_track_fraction, e0 → hillas_scale_energy, node_fields → shower_state_at_nodes, model → hillas_photon_yield, …); cphotang re-exports the old names for backward compatibility.

Key changes

  • Vectorized batch kernel: ~0.012–0.018 ms/shower in run(); the event loop, per-event Python overhead, and the C++ zsteps extension are gone (pure Python + NumPy; simplifies wheels/packaging).
  • Accurate energy quadrature (science impact — see below): the secondary-electron energy integral ∫ F(u(E))·(−dT/dE) dE is computed with two-panel Gauss-Legendre quadrature in ln E (3 low + 8 high nodes, per-(shower,node) bounds) instead of the legacy 9-bin decade Riemann sum.
  • Detector/cloud visibility window: the longitudinal integration window is bounded by detector altitude and cloud top analytically; _apply_clouds post-hoc masking removed.
  • Config: optional [simulation.cherenkov_quadrature] node-count knobs (defaults preserve behavior; older config files load unchanged); unit metadata declared once via Annotated types (~30 boilerplate validators/serializers removed); generic FITS config reader — fixes a latent bug where initial_position.longitude was read from the latitude key on FITS round-trip.
  • Radio: EAS radio path updated to consume the refactored interfaces.
  • Tau sampling (utils/cdf.py): inverse-CDF interpolation made robust at grid boundaries (bit-identical in the interior).
  • Tests: 170 pass — golden anchors, property tests, quadrature-convergence and propagation-length suites (addresses long-standing eas_optical unit tests #60 coverage gaps for eas_optical). Every modularization commit is regression-gated bit-identical.

⚠️ Science-result change (deliberate, documented)

This PR changes optical Cherenkov results relative to main. Main's legacy decade-grid energy sum evaluates the angular-capture CDF at each decade's linear midpoint, over-counting photons by ~1.6× on the energy axis. With the accurate integral, on a POEMMA-like 1e6-event run:

  • photon density dphots: ~0.62× main
  • Optical MC integral: ~0.50× main
  • passing events: ~0.55× main
  • GEO-only integral: ~0.87× main (via the √(2 ln numPEs) effective-Cherenkov-angle feedback)

This was isolated to the energy quadrature alone: substituting only the GL quadrature into main's own harness reproduces these ratios within ~2%. The bare Cherenkov angle distribution is unchanged. We consider the accurate integral the physically-correct choice; downstream sensitivity baselines should be re-derived after merge.

Compatibility

  • Config files: fully backward compatible (new sections optional).
  • CLI: unchanged.
  • Output FITS: same columns/keys.
  • Python API: cphotang re-exports all moved/renamed symbols under their old names.
  • Requires no compiled extensions.

🤖 Generated with Claude Code

Areustle and others added 3 commits July 30, 2026 11:44
@Areustle
cphotang refactor: propagation, quadrature, Hillas kernel, radio
866dfb5 
- Initial photon_sum refactor specification with tests
- Profiling: simplify profile_run.py for single-threaded run() timing
- perf(propagation): faster slant-depth inverse and atmosphere density
- refactor(cphotang): asymmetric energy quad, pluggable yield model, cang reformulation
- Kernel singularity analysis, CDF-transform exploration, run() narrative
- perf(cphotang.__call__): collapsed path + adaptive core-aware chunk size
- HPC: dask distributed local cluster startup early in compute
- Radio refactor
@Areustle
remove C++ dependencies
06dd59d
@Areustle @claude
config: dedupe unit metadata, generic FITS reader, quadrature knobs
e1e409a 
The config system declared structure and values redundantly. Three
structural fixes, all behavior-preserving (existing config tests pass
unchanged):

1. Unit metadata declared once. Add reusable Annotated quantity types
   (Kilometers, Radians, MegaHertz, Decibels, SquareMeters) via a
   _unit_float factory bundling BeforeValidator + PlainSerializer. Each
   unit-bearing field now names its unit exactly once (e.g.
   `altitude: Kilometers = 525.0`), replacing the prior trio of an
   explicit default + per-field field_validator + per-field
   field_serializer. Deletes ~30 boilerplate methods, including 4
   copy-pasted radians validator/serializer pairs.

2. config_from_fits is now the generic inverse of the writer. Replace the
   45-line hand-transcribed key list with unflatten_dict (new, beside
   flatten_dict) over the `Config <path>` HIERARCH headers. Auto-syncs
   with the model and reconstructs the full config (the old reader was
   lossy). Fixes a latent bug where the hand-written reader loaded
   initial_position latitude into longitude.

3. Add optional simulation.cherenkov_quadrature config (n_nodes,
   n_slant_sub, n_energy_low, n_energy_high) for the CphotAng GL kernel
   knobs. Defaults match run(), so omitting the section (any existing
   config file) reproduces current behavior. Wired config -> eas ->
   CphotAng.__call__ -> run(); sample_input_file.toml gets an example
   block.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
@Areustle
Areustle force-pushed the cphotang_refactor branch from cf19bff to e1e409a Compare July 30, 2026 15:44
Areustle and others added 7 commits July 30, 2026 13:23
@Areustle @claude
refactor: move BackgroundCluster to utils/distributed
eca964f 
Call-graph analysis shows BackgroundCluster has zero edges to any physics --
it is pure cluster-lifecycle infrastructure, yet compute.py imported it from
the cphotang physics module. Home it with the other cross-cutting utilities;
cphotang re-exports it for backward compatibility. No behavior change.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
@Areustle @claude
refactor: home ozone/aerosol models in atmospheric_models
eec9a54 
The ozone profile tables (OzZeta/OzDepth/OzDsum), aerosol OD table (aOD55),
their derived exact layer-walk tables, ozone_losses/ozone_rate, and the
analytic ozone_column/aerosol_column/refractive_index functions are all
atmosphere physics with no coupling to the Cherenkov orchestration -- the
call graph clusters them apart from CphotAng.run. Move them (verbatim) to
atmospheric_models.py as module-level functions and constants; CphotAng
keeps thin ozone_losses/ozone_rate delegates and re-exports for
compatibility. Physics goldens unchanged (bit-identical).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
@Areustle @claude
refactor: home per-node shower fields in shower_properties
0f10099 
shower_age (energy-aware t,s variant -> renamed shower_age_at_depth to
avoid colliding with the legacy fixed-1e8-GeV shower_age(T)),
greisen_particles, hillas_e2, and node_fields (renamed
shower_state_at_nodes -- it computes the per-node shower state, not
generic "fields") form a self-contained call-graph community with no
edges into the Cherenkov orchestration. Home them with the other
longitudinal-profile physics. Verbatim moves; goldens bit-identical.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
@Areustle @claude
refactor: home the photon-yield kernel adapter in hillas_batch_kernel
d9ca421 
hillas_single_integral_model, its nested closure, tracklen, and e0 cluster
with the batch kernel in the call graph -- the factory's only real callees
are precompute_kernel_NE16 / delta_nphots_single_integral_NE16, defined in
hillas_batch_kernel. Moving them there turns those cross-module edges into
local calls and gives the vague names honest ones:

  hillas_single_integral_model -> make_hillas_photon_yield  (factory)
  model (nested closure)       -> hillas_photon_yield
  tracklen                     -> secondary_track_fraction  (Hillas eqn 8 T(E))
  e0                           -> hillas_scale_energy       (E0(s) of eqn 8)

cphotang keeps backward-compatible aliases. Goldens bit-identical.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
@Areustle @claude
refactor: home the longitudinal grid in propagation
b44f5ba 
_propagation_grid/_valid_window/_node_grid/_node_geometry cluster with
propagation.py's functions in the call graph -- they are path geometry and
quadrature, not Cherenkov orchestration. Move them there as pure functions
(propagation_grid, visibility_window, gl_node_grid, node_geometry). The
three shower-physics depths they consumed (e2hill age-floor, Greisen death
depth, shower-maximum depth) are now explicit parameters computed by a thin
CphotAng._propagation_grid wrapper, so the propagation layer stays free of
shower models and no import cycle forms.

Also homes the generic cached_leggauss GL cache in quadrature.py, replacing
four deferred function-local imports propagation used to dodge the old
layering. Goldens bit-identical.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
@Areustle @claude
refactor: extract distributed dispatch; rename opaque photon stages
3afc500 
__call__'s body was an isolated call-graph community: chunk sizing, cluster
lifecycle, map_blocks fan-out, and the progress bar -- pure execution
machinery with no physics edges. It now lives in dispatch.py as
map_showers_distributed(chunk_fn, ...); CphotAng.__call__ keeps only the
physics packing (its run() closure) and delegates. All dask/rich imports
leave cphotang.

The three "photon" stages whose names didn't distinguish them get honest
ones: _photon_yield -> _transmitted_yield (per-node transmitted yield
weights), _photon_sum -> _kernel_energy_integral (the energy-quadrature
kernel evaluation), _photon_density -> _density_at_detector.

cphotang.py: 1583 -> 958 lines across the six moves. Goldens bit-identical.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
@Areustle @claude
refactor: unify us_std_atm_density on the propagation implementation
1db0c6d 
atmospheric_models carried a second, legacy copy of the US Standard
Atmosphere density (ratio-scaled base pressures, searchsorted layer
lookup, boolean-mask pressure branches). Its only remaining caller on
this branch was scripts/atmosphere_rho_curve_fit.py; the vectorized
cphotang path already used propagation's version. The two agree to
<2e-7 relative on a dense 0-90 km grid (the residual is rounding of
the base-pressure table); the legacy copy additionally wrapped z < 0
to the top layer via searchsorted - 1, which the propagation version
handles correctly.

Delete the legacy def and re-export propagation's from
atmospheric_models for compatibility; pin the unification with an
identity test.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

No reviews

Assignees

No one assigned

Labels

None yet

Projects

None yet

Milestone

No milestone

Development

Successfully merging this pull request may close these issues.

EAS Refactor & Performance Improvement Modularize cphotang EAS for CHASM and other implementations

1 participant

@Areustle