Classification improvement and robustness.

[GernotMaier]

VERITAS Gamma/Hadron Classification Review and Implementation Plan

Implementation status

The code-only hardening described below is implemented in the standard
eventdisplay-ml-train-xgb-classify path. It uses only existing ROOT products
and configured input lists: no simulation production, background generation,
data reprocessing, or TMVA-style behavior was added. Runtime changes include
the robust/extended profile selector (including coarse zenith conditioning), provenance-aware capped sampling,
telescope-schema validation, inactive-telescope masking, grouped splitting when
enough source files exist, reserved gamma threshold calibration, capped zenith
weights, explicit feature-schema failures, and conservative background upper
limit diagnostics. Sparse file groups are recorded as an event-level fallback
instead of being presented as independent validation.

Operational validation remains: run the repository tests in the project
environment and inspect serialized support/calibration metadata for each
existing VERITAS model bin before deployment.

Regression isolation

Classification safeguards in the shared loader/flattening module are explicitly gated on
analysis_type == "classification". Regression retains its historical telescope sorting,
inactive-slot values, event-cap quota, random stream, telescope-configuration update behavior,
and exception wrapping. The optional sorting argument preserves the original four-argument
regression helper call. Regression application and stereo-training tests were run after this
isolation pass.

Scope and fixed constraints

This review covers only the standard VERITAS classification path invoked through
eventdisplay-ml-train-xgb-classify.

The implementation must use the files that already exist:

• signal consists of optimal-condition gamma-ray simulations;
• background consists of observed events and is sparse in parts of energy, zenith, and the
  wider observing-condition space;
• no new simulations, simulation reruns, low-level reprocessing, externally prepared samples,
  or new observational data are available;
• all recommendations must therefore be implementable in this repository while reading the
  existing ROOT files and file lists.

The code cannot prove gamma efficiency under observational conditions absent from the signal
simulation. It can, however, substantially reduce the ability to learn observational shortcuts,
use scarce background efficiently, detect weak support, and avoid representing extrapolation as
measured performance. Those are the goals of this plan.

Nothing in this plan requires:

• new or rerun simulations;
• adding NSB/noise or throughput effects to existing simulations;
• low-level waveform, image-cleaning, or reconstruction work;
• new background observations or newly prepared OFF regions;
• Crab or other real-gamma validation data;
• manual conversion of the existing inputs into a new format.

An optional flattened cache is produced transparently by the training code from the same input
files and is only a computational optimization.

Recommended architecture

The production classifier should be a composed system rather than an unconstrained XGBoost
score:

existing ROOT files
    -> validated deterministic flattening
    -> provenance-aware grouped folds
    -> robust morphology classifier
    -> final-model gamma-efficiency calibration
    -> out-of-fold background-efficiency estimate
    -> support assessment
    -> supported / borrowed / fallback decision

The important separation of responsibilities is:

• morphology features are allowed to determine gamma/hadron likeness;
• nuisance and provenance variables control splitting, weighting, diagnostics, support,
  and thresholds but are excluded from the baseline classifier;
• deterministic preselection handles invalid and trivially background-like events;
• fallback logic still classifies events in sparse regions without claiming unsupported
  rejection performance.

Summary of current risks

Critical: random event splitting measures domain recognition as valid performance

Signal simulation and observed background are concatenated and randomly split by event. Events
from the same files, runs, simulation production, and conditions therefore occur in training and
test samples. The nominal test sample reproduces the same simulation-versus-observation
differences as training.

This is especially unsafe with optimal-condition signal MC. A stable random-split AUC can result
from NSB, throughput, trigger, cleaning, telescope activity, reconstruction failure, missingness,
or production differences rather than particle morphology.

Critical: the nominal test sample selects the model and calibrates it

The nominal test set is the last XGBoost eval_set and controls early stopping. The same rows
are then used for reported metrics, SHAP, efficiency curves, and thresholds. It is a validation
sample, not an independent test sample.

Critical: class-dependent schemas and missingness can encode the label

Signal and background are loaded independently. Each class drops its own all-NaN columns
before concatenation. A feature present only in one class is restored by concatenation and is
missing for the other class. XGBoost can use that missingness as a nearly perfect label.

No common-schema check currently compares branch availability, flattened columns, dtypes,
missing fractions, telescope IDs, mirror areas, or processing configuration.

Critical: the current feature set exposes observation-sensitive handles

The full per-telescope feature set contains physically useful variables, but also strong
condition sensitivity:

• tel_active_* and missing-value patterns depend on trigger, cleaning, disabled hardware,
  and threshold behavior;
• loss, timing gradient, width, and length are sensitive to NSB and image cleaning;
• size-based ordering changes when images fluctuate across cleaning thresholds;
• ze_bin can identify class composition when gamma MC and background zenith populations
  differ;
• configuration-dependent constants and telescope geometry can identify an epoch or file set.

Removing only Erec does not remove energy information: morphology, multiplicity, loss, and
size ordering are all energy-dependent.

High: zenith weighting cannot repair missing support

The current optional weighting equalizes class zenith fractions only where both classes exist.
If one class is absent from a bin, the code logs the condition and continues while ze_bin
remains a classifier input. In such a cell, zenith can become a direct class label.

Weights are uncapped, no minimum effective count is required, and validation metrics do not use
the corresponding validation weights.

High: sparse background makes current efficiencies and thresholds unstable

Efficiency is evaluated on one random split and a grid of only 101 score thresholds. Sparse
high-energy or zenith subsets can have few or zero surviving background events. The code does
not store a confidence interval, effective count, or support state. Zero observed survivors can
therefore look like zero background efficiency.

High: preselection and ML domains are not composed explicitly

Training removes events outside the MSCW, MSCL, emission-height, and selected energy window.
This is efficient and physically reasonable for trivially background-like or invalid events,
but it means the ML score is not validated outside that domain. A complete classification path
must define what happens there instead of implicitly extrapolating the model.

High: background provenance and source contamination are not checked

The trainer accepts background file lists but does not know whether events are from source-free
regions, whether exclusion regions were applied, or whether the same events appear in several
files/lists. No new preparation is proposed, but the code should inventory what it receives and
make duplication and file-level dependence visible.

If the existing files do not contain direction/region or run identifiers, the trainer must record
that provenance and source-exclusion status are unknown. It should not require new files or
pretend that source contamination was tested.

High: max_events can unexpectedly load all events

The cap is divided by the number of files using integer division. If max_events is smaller than
the file count, the per-file cap is zero; zero is treated as false and all events are loaded.
Equal per-file quotas also discard the remainder and distort the natural mixture of file sizes.

Medium: inactive VERITAS telescope entries are not explicitly masked

An activity mask is constructed from DispTelList_T, but fixed-index telescope variables are
not masked with it. Correct behavior depends on inactive ROOT entries already being NaN.
Zeros or sentinels can become real features and simulation/data shortcuts.

Medium: telescope configuration changes are incompletely validated

Configuration is replaced only when a later file has a larger maximum telescope ID. Changes in
IDs, mirror areas, or positions with the same maximum ID are ignored. Signal and background
loading also overwrite the shared saved configuration independently.

Medium: zenith definitions silently clip overflow

Angles below or above the configured range are clipped into the first or last bin. Empty,
unsorted, overlapping, or gapped definitions are not rejected.

Medium: independent energy-model scores are not guaranteed to share calibration

Each energy-window model has its own class ratio, weights, sample composition, and early-stopped
iteration. Its logistic output is a discriminant score, not a physical gamma probability.
Linear interpolation assumes comparable score calibration without establishing it.

Medium: invalid values can be converted into boundary features

Some sentinels and reconstruction failures are clipped to finite boundaries rather than changed
to NaN. Boundary piles can differ by class and are easy for trees to recognize.

Medium: preprocessing is repeated for every overlapping energy bin

Every bin invocation rereads and flattens the same files. This is expensive and can give
neighboring models different random samples of the same scarce background.

Detailed implementation plan

Phase 1: make input handling deterministic and auditable

1.1 Validate classification arguments

Update configure_training("classification") to:

• require input_signal_file_list, input_background_file_list, and model_parameters;
• require 0 < train_test_fraction < 1 until that option is replaced by grouped-fold options;
• reject negative max_events, non-positive core counts, and invalid read sizes;
• use a deterministic default random_state in production configuration;
• validate every energy bin has finite E_min < E_max;
• validate zenith bins are finite, ordered, non-overlapping, and contiguous;
• distinguish configuration, schema, ROOT-read, and data-validation exceptions instead of
  wrapping all failures as FileNotFoundError.

Add classification-specific options:

• --feature_profile robust|extended, default robust;
• --group_folds, default 5;
• --weight_cap, default selected conservatively, for example 10;
• --minimum_background_effective_count;
• --minimum_signal_effective_count;
• --support_policy conservative|warn, default conservative;
• --classification_cache, optional internal cache path;
• --run_crossfit_diagnostics enabled by default for production training.

1.2 Preserve internal provenance

During loading, add reserved columns that are never passed to XGBoost:

• __class_label;
• __input_file_id, a stable index or hash rather than an absolute path;
• __entry_index, the original ROOT entry number;
• __sample_id = (__input_file_id, __entry_index);
• __energy_bin and __ze_bin for routing and diagnostics;
• __preselection_state;
• available condition values already present in the ROOT file.

If run or shower identifiers already exist in the input tree, use them. Do not require new
branches. If they do not exist, file ID is the strongest available grouping key. When only one
file exists for a class, use deterministic contiguous entry blocks as a documented weaker
fallback and mark the grouping quality in metadata.

Before training, reject duplicate __sample_id values and report repeated physical identifiers
when available.

1.3 Establish one schema before class concatenation

Do not drop all-NaN columns independently. Instead:

1. inspect all signal and background files before flattening;
2. create the expected common branch and telescope schema;
3. compare complete telescope configurations for every file;
4. flatten both classes to the same ordered columns;
5. calculate missing and clipped fractions by class and file;
6. drop a feature only through one common decision;
7. fail if a required robust feature is absent or all-NaN in either class;
8. quarantine extended features with extreme class-dependent missingness.

Persist the schema, feature order, transformation version, and a configuration hash in the
model.

1.4 Correct VERITAS activity and sentinel handling

Before size sorting or derived-feature calculation:

• use DispTelList_T to set every inactive telescope measurement to NaN;
• validate telescope IDs are non-negative and present in telconfig;
• validate fixed-array widths against configured IDs;
• translate known Eventdisplay sentinels to NaN before clipping;
• keep physical clipping separate from invalid-value conversion;
• report boundary-clipped and missing fractions by class and file.

The tel_active_* columns may remain available for diagnostics but should not be part of the
default robust feature profile.

1.5 Fix exact and efficient sampling

Replace the per-file integer quota with an exact streaming reservoir over each complete class
file list. The reservoir must:

• return exactly min(max_events, available_events) after preselection;
• use one deterministic random stream per class derived from the configured seed;
• preserve the natural contribution of file sizes unless explicit stratified sampling is
  requested;
• retain provenance for every selected row;
• never interpret a zero quota as unlimited.

For sparse background, allow separate caps for signal and background. Never discard background
merely to make raw class counts equal; control class influence with weights.

Phase 2: introduce a robust feature contract

2.1 Define a robust baseline profile

The default profile should contain physically motivated shower morphology variables already
available after current preprocessing:

• MSCW;
• MSCL;
• log-transformed energy-consistency quantity;
• emission height;
• log-transformed emission-height consistency;
• core distance, after sentinel validation.

It also includes the per-telescope image-morphology parameters needed for gamma/hadron
separation: image size, width, length, asymmetry, loss, distance, image orientation
(cosphi, sinphi), and time-gradient (tgrad_x) columns. These are retained because they
carry the primary shower-shape information. They are not the same as detector-condition
proxies.

Core distance should be ablated from the minimal profile as well as tested in it, because its
distribution depends on collection geometry and trigger/selection acceptance in addition to
shower reconstruction.

The coarse ze_bin is included as a conditioning feature because atmospheric depth and
projection change many image parameters. It is not treated as a physical gamma/hadron
discriminator: class/zenith balancing, a ze_bin-only negative-control diagnostic, and
per-zenith performance gates are required. --ignore_ze_bin remains available for ablation.

The baseline must exclude:

• input file, run, or production identity;
• telescope configuration constants;
• tel_active_*;
• mirror-area and telescope-position geometry (mirror_area_*, tel_rel_*, R_core_*);
• explicit missingness masks;
• raw condition or pointing variables;
• reconstructed energy as a direct discriminator.

This does not make the baseline immune to domain shift, because morphology itself can shift,
but it removes the easiest observation-condition shortcuts.

2.2 Keep the current expanded profile as an ablation

The current per-telescope feature set becomes extended, not the unexamined default. Train both
profiles with identical folds and weights. The extended profile is eligible for production only
if it improves all of the following:

• grouped out-of-fold background efficiency at fixed gamma efficiency;
• worst-file and worst-zenith performance;
• stability of background score distributions across files after energy/zenith conditioning;
• missingness-only and nuisance-only negative-control results;
• ensemble/fold stability in sparse cells.

A gain in random-split AUC is not sufficient.

2.3 Add negative-control classifiers

Using the same grouped folds, automatically train small diagnostic classifiers from:

1. ze_bin alone;
2. telescope activity and feature-missingness masks only;
3. provenance grouping proxies that can be represented without using unique file identifiers as
   categorical labels;
4. energy and zenith routing variables only.

These models are never serialized for application. Store their grouped AUC and feature
importance. Strong separation is evidence that the main classifier has a domain shortcut
available. It does not prove that every morphology feature is invalid, but it blocks claims
based solely on inclusive performance.

Phase 3: replace current zenith weighting with regularized balancing

3.1 Build training-only balance cells

Within each selected model energy window, construct coarse cells from:

• a finer reconstructed-energy index;
• broad zenith bins;
• optionally telescope multiplicity for diagnostics, not necessarily for weighting.

Compute the table from the training fold only. Merge adjacent cells deterministically until
both classes reach a configurable minimum or no further merge is possible.

3.2 Use shrunk and capped weights

For cells containing both classes:

• define an explicit reference distribution, preferably the average of normalized signal and
  background cell fractions;
• calculate class/cell weights toward that reference;
• shrink weights toward 1 as cell counts decrease;
• cap weights with weight_cap;
• normalize total signal and background weight to equal totals unless another prior is
  explicitly configured.

For a cell lacking one class:

• do not fabricate a compensating weight;
• retain only the coarse ze_bin conditioning variable when support is adequate;
• mark the cell as weak support and require the ze_bin-only diagnostic to be reported;
• allow its events to contribute to the pooled class morphology distribution with capped base
  weight;
• mark the cell as weak support for calibration and reporting.

Report, by class and cell:

• raw count;
• sum of weights;
• maximum weight;
• effective sample size (sum(w) ** 2) / sum(w ** 2);
• whether the cell was merged or weakly supported.

3.3 Make early-stopping evaluation consistent

Derive validation weights using mappings fitted on the training fold only. Pass corresponding
weights through sample_weight_eval_set. Store both:

• reference-weighted metrics, which test the intended morphology comparison;
• natural-distribution metrics, which describe the supplied background sample.

Phase 4: use scarce background with grouped cross-fitting

4.1 Construct grouped folds

Use StratifiedGroupKFold where the available groups permit it. Groups must be based on file,
run, or simulation identity—not random rows. Reduce the number of folds automatically if a
class has too few groups, and record the reduction.

Each fold must report counts by class, energy sub-bin, zenith bin, and file. Empty evaluation
cells are marked undefined; they are never reported as zero efficiency.

4.2 Separate three distinct products

Training should produce:

1. cross-fitted assessment: every background event and every non-reserved gamma event gets
   one prediction from a model that did not train on its group;
2. final production model: trained after hyperparameters and tree count are fixed;
3. final-model gamma calibration: a deterministic gamma subset never used to fit the final
   model, used to derive its score thresholds.

Gamma simulation is less scarce, so reserving gamma calibration costs much less than reserving
background. All available background can train the final model after its out-of-fold assessment
prediction has been produced.

4.3 Select tree count without reusing assessment scores

For each grouped fold, use an internal validation partition for early stopping. Aggregate the
best iterations, for example with a robust median. Train the final model with that fixed tree
count and no production early-stopping sample.

Hyperparameter tuning, if enabled, must be nested or restricted to designated development
folds. The aggregated out-of-fold assessment must not select hyperparameters.

4.4 Evaluate operating points without transferring raw thresholds

Raw thresholds from fold models do not necessarily transfer to the final model. Instead:

• in each assessment fold, derive the threshold for a target gamma efficiency from that fold's
  held-out gamma scores;
• count held-out background survivors at that fold-specific threshold;
• aggregate survivor and trial counts across folds;
• train the final model;
• derive its production threshold at the same target efficiency from the reserved gamma
  calibration subset.

This uses all scarce background for honest performance estimation while calibrating the actual
deployed model on abundant held-out gamma MC.

Phase 5: make sparse-region calibration conservative

5.1 Replace the 101-point threshold scan

For each target gamma efficiency, use exact empirical gamma-score quantiles with a documented
tie convention. Store the threshold, gamma count, and quantile uncertainty.

Do not describe predict_proba as a physical probability unless a separate calibration test
supports that interpretation. Treat it as a monotonic classifier score.

5.2 Estimate background efficiency with group-aware uncertainty

For every energy/zenith region and operating point, store:

• background trials and survivors from out-of-fold predictions;
• raw and effective weighted counts;
• point estimate;
• conservative confidence interval;
• number of contributing files/runs;
• between-group dispersion.

Use grouped bootstrap over files/runs as the primary uncertainty estimate. For unweighted or
approximately unweighted counts, also store a binomial interval. Zero survivors must yield a
finite upper limit, never a zero-efficiency claim.

5.3 Borrow strength without pretending to have data

Smooth threshold and background-efficiency summaries across neighboring energy and zenith
regions using a low-complexity, regularized surface. For sparse cells:

• shrink toward neighboring and inclusive estimates;
• enlarge the uncertainty;
• record which cells contributed;
• forbid extrapolation beyond the configured energy/zenith envelope without fallback status.

The raw cell counts and unsmoothed estimates must remain in the model metadata so smoothing
cannot hide the lack of data.

5.4 Define support states

Assign each model region one of:

• supported: minimum effective signal/background counts and group coverage are met;
• borrowed: the pooled classifier is usable, but calibration relies on regularized neighboring
  information;
• fallback: calibration support is insufficient even after pooling.

Threshold choice should use the upper confidence bound on background efficiency in borrowed
regions. Significance optimization should include a configurable background-systematic floor
and minimum event requirements.

Phase 6: guarantee a classification path everywhere required

6.1 Compose deterministic and ML decisions

For every input row:

1. invalid reconstruction or failed quality requirements produce invalid_quality;
2. events outside loose gamma-like MSCW/MSCL preselection produce preselection_background;
3. events inside the per-energy model domain use that robust classifier;
4. events in a sparse cell use the same pooled score with borrowed calibration;
5. events without defensible local calibration use a global robust fallback model or the
   established deterministic MSCW/MSCL selection;
6. positive energies outside the trained model envelope use fallback, not the nearest model
   with an unqualified nominal threshold.

Preserve output row counts. Add output or accompanying metadata for:

• classifier score;
• gamma decision;
• decision origin;
• support state;
• model/bin identifier;
• optional fold/ensemble spread if an ensemble is retained.

6.2 Train a global fallback from the same data

Train one compact robust classifier over the complete configured energy range using the same
files, grouped folds, and capped energy/zenith weights. It may include transformed reconstructed
energy only as a conditioning variable after balancing; compare both with and without it.

This global model is not expected to outperform local models. Its purpose is stable behavior
where a local energy model has insufficient background. If even the global support criteria
fail, use the deterministic classical selection and mark the result fallback.

Phase 7: add robustness diagnostics and release gates

7.1 Required diagnostic tables

Every training run should save:

• class counts and weighted effective counts by energy, zenith, file, and fold;
• common schema and missingness comparison;
• clipped/sentinel fractions;
• training and validation curves;
• grouped out-of-fold ROC and operating-point tables;
• per-file and per-zenith background survival;
• negative-control AUCs;
• core-versus-extended feature comparison;
• feature and SHAP importance from held-out rows only;
• support-state map;
• final-model gamma thresholds and their calibration counts.

7.2 Stability tests using only existing background

At fixed reconstructed-energy and zenith strata, compare background score and survival across:

• input files;
• available pointing/azimuth values;
• multiplicity;
• telescope activity state;
• missingness burden;
• any condition proxy already available in the ROOT files.

Use grouped bootstrap intervals and report the worst group. A smooth inclusive curve is not a
substitute for this test.

7.3 Production acceptance gates

A model should not be marked production-ready unless:

1. no training/assessment group overlap exists;
2. required schemas and telescope configurations match;
3. weight caps and minimum effective counts are satisfied or the cell is marked borrowed or
   fallback;
4. zero-survivor cells have finite upper limits;
5. the robust profile passes missingness/nuisance negative controls;
6. extended features, if used, improve worst-group results without increasing shortcut
   diagnostics materially;
7. all thresholds are calibrated on predictions from models that did not train on those gamma
   events;
8. the final model stores preprocessing, fold, weighting, feature, and support metadata;
9. repeated training with the same seed reproduces folds, samples, thresholds, and predictions;
10. every event follows a tested supported, borrowed, or fallback path.

Phase 8: improve computational efficiency

8.1 Flatten once, train all bins

Add an internal, versioned columnar cache produced automatically by the training command. This
does not require external data preparation. The command should:

1. stream each existing ROOT file once;
2. validate and flatten it;
3. store features, routing quantities, and hashed provenance;
4. partition by class and coarse reconstructed energy;
5. reuse the cache for all overlapping energy models and cross-fitting folds.

Reject a cache when input file fingerprints, preprocessing version, feature profile, telescope
configuration, or cut configuration changes.

8.2 Bound diagnostic cost

• implement classification support for eval_max_events using deterministic group-aware
  subsampling;
• cap SHAP rows per fold and sample them reproducibly;
• train negative controls with shallow trees and early stopping;
• parallelize energy bins or folds, but avoid multiplying XGBoost threads by process count;
• keep the final production model compact by fixing its tree count from fold results;
• log peak memory and elapsed time for loading, cache creation, folds, and final fitting.

Module-level change map

config.py

• add and validate classification robustness options;
• validate bin definitions and required inputs;
• set deterministic seeds;
• distinguish feature, routing, calibration, and provenance configuration.

data_processing.py

• attach reserved provenance and original entry indices;
• build one common schema across both classes;
• validate every telescope configuration;
• mask inactive fixed-index VERITAS values;
• normalize sentinels before clipping;
• implement exact class-level reservoir sampling;
• optionally write/read the internal versioned cache;
• return invalid zenith/energy states rather than clipping silently.

features.py

• define named robust and extended classification profiles;
• keep nuisance/routing variables separate from XGBoost inputs;
• document every feature's physical role and observation sensitivity.

models.py

• replace event-wise split with grouped folds;
• implement training-fold-only regularized weights;
• pass validation weights to XGBoost;
• generate out-of-fold predictions and aggregate best iterations;
• train the final local and global fallback models;
• derive final thresholds from reserved gamma calibration;
• store counts, intervals, support state, preprocessing hash, and diagnostics;
• prevent raw threshold transfer between fold and final models.

evaluate.py

• replace the fixed threshold grid with exact operating points;
• calculate grouped bootstrap and binomial intervals;
• represent undefined metrics as NaN;
• calculate worst-group and between-group stability;
• evaluate negative controls and robust/extended ablations.

Model loading and classification application

• validate the stored preprocessing/feature contract;
• load local and global fallback models;
• route supported, borrowed, and fallback events explicitly;
• preserve row counts and expose decision origin/support status;
• reject missing required features instead of silently creating an all-NaN matrix.

Required tests

1. Exact max_events behavior, including fewer requested events than files.
2. Reproducible class sampling and independent random streams.
3. Common-schema failure when one class lacks a required feature.
4. Detection of class-dependent all-NaN and extreme missingness columns.
5. Telescope-configuration mismatch detection.
6. Inactive fixed-index VERITAS entries converted to NaN.
7. Sentinel conversion preceding clipping.
8. Invalid, empty, overlapping, gapped, and out-of-range zenith bins.
9. No group overlap across training, validation, gamma calibration, and assessment.
10. One honest out-of-fold assessment score per background event.
11. Fold-count reduction when too few files/groups exist.
12. Training-only weight estimation with no fold leakage.
13. Weight cap, equal total class weight, and effective-count calculations.
14. One-class sparse cells marked borrowed/fallback without crashing training.
15. Validation weights passed consistently to XGBoost.
16. Exact gamma-score quantile thresholds with ties.
17. Zero background survivors producing a non-zero upper limit.
18. Grouped bootstrap reproducibility.
19. Raw fold thresholds never copied to the final model.
20. Final thresholds derived only from reserved gamma events.
21. Robust feature profile excludes nuisance and activity columns.
22. Nuisance-only and missingness-only diagnostic models.
23. Deterministic preselection, supported, borrowed, and fallback routing.
24. Global fallback behavior outside reliable local bins.
25. Output row-count preservation and decision-origin metadata.
26. Cache invalidation after input/configuration/preprocessing changes.
27. End-to-end synthetic VERITAS ROOT -> train -> save -> load -> classify contract.

Irreducible limitation

With only optimal-condition gamma simulation, no code-only method can measure the true gamma
efficiency under all real observing conditions. The proposed implementation must therefore not
claim that it has removed simulation-to-data systematics. It provides the strongest defensible
use of the available samples by:

• removing direct nuisance handles from the baseline classifier;
• detecting missingness and condition shortcuts;
• testing generalization across existing background groups;
• using every scarce background event out of fold;
• assigning conservative uncertainty and fallback states;
• never equating absence of observed background survivors with proven rejection.

This limitation should be serialized in model metadata and documented with every performance
result.

Reference basis

The VERITAS BDT study by Krause, Pueschel & Maier
(2017) demonstrates why condition matching matters for IACT
classification: NSB changes discriminating variables, major hardware epochs affect training
distributions, and condition subdivision becomes limited by available background. Because the
present constraints do not allow matched gamma simulation, this plan substitutes feature
restriction, grouped assessment, uncertainty, and fallback behavior; it does not claim those
steps reproduce condition-matched MC.

Verification limitations

This is a source-level and methodological review. Representative ignored VERITAS artifacts
were not present in this checkout. The complete test suite could not be executed because the
named project environment was unavailable and the other available environments lacked required
dependencies. No package source code was modified as part of this review.

[Copilot]

🟡 Not ready to approve

Two issues need fixing: the “extended” feature profile can still leak tel_active_* into training/apply, and _efficiency_dataframe can miscompute the 95% background-efficiency upper bound due to float-derived survivor rounding.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

This review doesn't count toward merge requirements. Sign up for the private preview to control whether Copilot approvals count.

Pull request overview

This PR hardens the gamma/hadron classification pipeline by preventing provenance/routing leakage into training features, improving split robustness (optionally grouped by source file), adding held-out signal threshold calibration, and tightening data/metadata validation during training and apply.

Changes:

• Add feature-profile selection for classification (“robust” vs “extended”) and corresponding unit tests for the hardening contract.
• Update classification training to use a train/validation/test split (optionally grouped by source file), record split/support metadata, and compute held-out signal threshold calibrations.
• Improve robustness in preprocessing and evaluation: mask inactive telescope slots for classification, validate telescope-config consistency across files, and add background-efficiency upper limits to efficiency diagnostics.

File summaries

File	Description
tests/test_classification_robustness.py	Adds focused regression tests for feature-profile filtering, grouped splitting, and threshold-calibration behavior.
src/eventdisplay_ml/models.py	Implements grouped/stratified splitting, held-out signal threshold calibration, stricter feature-schema checks on apply, and adds classifier nuisance diagnostics/support metadata.
src/eventdisplay_ml/features.py	Introduces classification feature-profile selection to exclude provenance/routing columns from model inputs.
src/eventdisplay_ml/evaluate.py	Adds quantile-based threshold calibration helper and computes an upper confidence bound for background efficiency.
src/eventdisplay_ml/data_processing.py	Masks inactive telescope slots for classification, adds stricter input/config validation for classification training, and strengthens zenith-bin input validation.
src/eventdisplay_ml/config.py	Adds CLI flags to select classification feature profile and enable/disable grouped splits, and logs the chosen settings.

Review details

• Files reviewed: 6/6 changed files
• Comments generated: 2
• Review effort level: Lite

We're testing this review assessment. Please use 👍 or 👎 to tell us if it's correct.

[Copilot]

🟡 Not ready to approve

The new split helper can fail with cryptic sklearn errors for some valid train_test_fraction/small-sample scenarios and needs explicit guards to make failures deterministic and actionable.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

This review doesn't count toward merge requirements. Sign up for the private preview to control whether Copilot approvals count.

Review details

Suppressed comments (3)

src/eventdisplay_ml/models.py:1297

• In the stratified event-level fallback, splitting label_hold into validation/test will fail when a class has too few events (e.g., 2 events total with train_fraction=0.5 leaves 1 holdout event). Add a clear guard so small synthetic inputs fail with a helpful error instead of a sklearn exception.

            label_validation, label_test = train_test_split(
                label_hold,
                train_size=0.5,
                random_state=rng,
            )

src/eventdisplay_ml/models.py:1280

• In grouped mode, train_test_split(g_hold, ...) can raise a cryptic ValueError when train_test_fraction leaves fewer than 2 holdout groups (e.g., 6 groups with train_fraction=0.9). Add an explicit guard so the failure is deterministic and actionable (or fall back to event-level splitting).

This issue also appears on line 1293 of the same file.

            g_validation, g_test = train_test_split(
                g_hold,
                train_size=0.5,
                random_state=rng,
            )

src/eventdisplay_ml/models.py:1327

• _classification_nuisance_diagnostics only considers columns ending in _<0..63> when computing feature_missing_fraction, so any telescope-indexed columns with suffixes >=64 are silently ignored (possible for CTAO / larger arrays). Use a regex for a trailing numeric suffix instead of a fixed range.

    telescope_columns = [
        column
        for column in df.columns
        if column.endswith(tuple(f"_{index}" for index in range(64)))
    ]

• Files reviewed: 7/7 changed files
• Comments generated: 0 new
• Review effort level: Lite

We're testing this review assessment. Please use 👍 or 👎 to tell us if it's correct.

[Copilot]

🟡 Not ready to approve

Confirmed issues in the updated classification weighting and missing-bin application logic can produce incorrect behavior (cap not enforced; borrowed-bin interpolation inconsistencies).

Once you've addressed the issues Copilot identified, you can request another Copilot review.

This review doesn't count toward merge requirements. Sign up for the private preview to control whether Copilot approvals count.

Review details

Suppressed comments (5)

src/eventdisplay_ml/models.py:1449

• _class_zenith_balance_weights clips weights to weight_cap and then normalizes by the mean, which can push some weights back above weight_cap (so the cap is not actually enforced on the returned weights). If the cap is meant as a hard robustness guard, apply the final clip after any normalization step (or remove mean-normalization).

    weight_values = np.clip(weights.to_numpy(dtype=np.float32), 0.0, float(weight_cap))
    mean_weight = weight_values.mean()
    if mean_weight > 0:
        weight_values /= mean_weight

src/eventdisplay_ml/models.py:590

• Related to missing-bin borrowing: threshold interpolation currently uses e_bin_lo == e_bin_hi (requested bins) to decide whether to blend thresholds, even if one/both bins were resolved to a different (borrowed) model. This can yield thresholds that are inconsistent with how the score was computed when bins are missing.

        thresholds_lo = models[resolved_lo].get("thresholds", {})
        thresholds_hi = models[resolved_hi].get("thresholds", {})
        for eff in threshold_keys:
            if eff in is_gamma:
                thr_lo = thresholds_lo.get(eff)
                if thr_lo is None:
                    continue
                if e_bin_lo == e_bin_hi:
                    threshold = thr_lo
                else:
                    thr_hi = thresholds_hi.get(eff)
                    if thr_hi is None:
                        continue
                    alpha = group_df["e_alpha"].to_numpy(dtype=np.float32)
                    threshold = (1.0 - alpha) * thr_lo + alpha * thr_hi
                is_gamma[eff][group_df.index] = (class_probs >= threshold).astype(np.uint8)

src/eventdisplay_ml/models.py:1172

• weights_validation is recomputed from the validation split’s own label/zenith distribution. For a balancing scheme intended to define a fixed target distribution, validation weights should be derived from the training split only and then applied to validation; otherwise early stopping/eval metrics depend on validation composition and can vary run-to-run in a way unrelated to model quality.

    if model_configs.get("balance_class_zenith_weights", False):
        weights_train = _class_zenith_balance_weights(full_df.iloc[train_idx], y_train)
        weights_validation = _class_zenith_balance_weights(
            full_df.iloc[validation_idx], y_validation
        )

src/eventdisplay_ml/data_processing.py:1572

• zenith_in_bins still clips out-of-range zenith angles into the first/last bin via np.clip(...). The PR description calls out returning invalid zenith states rather than clipping silently; if that’s a requirement for the hardening contract, callers likely need an explicit invalid marker (or an exception) instead of silent clipping.

    if bins.ndim != 1 or len(bins) < 2 or not np.all(np.isfinite(bins)):
        raise ValueError("Zenith-bin edges must be a finite one-dimensional sequence.")
    if np.any(np.diff(bins) <= 0):
        raise ValueError("Zenith-bin edges must be strictly increasing.")
    idx = np.clip(np.digitize(zenith_angles, bins) - 1, 0, len(bins) - 2)
    return idx.astype(np.int32)

src/eventdisplay_ml/models.py:564

• When an energy-bin model is missing, _resolve_classification_bin borrows the nearest available model, but the score interpolation still uses e_bin_lo/e_bin_hi and e_alpha computed from the requested bin centers. This can blend models whose calibration/energy-centers don’t match the interpolation weights, producing inconsistent scores in sparse/missing-bin cases.

This issue also appears on line 575 of the same file.

        resolved_lo = _resolve_classification_bin(models, e_bin_lo)
        resolved_hi = _resolve_classification_bin(models, e_bin_hi)
        model_lo = models[resolved_lo]["model"]
        model_hi = models[resolved_hi]["model"]
        missing_lo = sorted(set(models[resolved_lo]["features"]) - set(flatten_data.columns))
        missing_hi = sorted(set(models[resolved_hi]["features"]) - set(flatten_data.columns))
        if missing_lo or missing_hi:
            raise ValueError(
                "Classification model/input feature schema mismatch: "
                f"low-bin missing={missing_lo}, high-bin missing={missing_hi}."
            )
        flatten_lo = flatten_data.loc[:, models[resolved_lo]["features"]]
        flatten_hi = flatten_data.loc[:, models[resolved_hi]["features"]]

• Files reviewed: 7/7 changed files
• Comments generated: 0 new
• Review effort level: Lite

We're testing this review assessment. Please use 👍 or 👎 to tell us if it's correct.

[Copilot]

🟡 Not ready to approve

Two confirmed edge-case bugs in new zenith dict-bin handling and grouped splitting can silently mis-bin zenith values or break group-disjoint splits under overlapping group IDs.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

This review doesn't count toward merge requirements. Sign up for the private preview to control whether Copilot approvals count.

Review details

Suppressed comments (2)

src/eventdisplay_ml/data_processing.py:1571

• zenith_in_bins() accepts dict-style bins but currently converts them to numeric edges by taking all Ze_min values plus the last Ze_max. This implicitly assumes bins are contiguous, ordered, and non-overlapping; if a config provides gapped/overlapping/out-of-order dict bins, the conversion can silently produce incorrect bin edges and wrong ze_bin assignments. Consider validating the dict bins (finite, Ze_min < Ze_max, contiguous) and constructing edges from the dicts’ boundaries explicitly so invalid definitions fail fast.

    if isinstance(bins[0], dict):
        if any("Ze_min" not in b or "Ze_max" not in b for b in bins):
            raise ValueError("Zenith-bin dictionaries require Ze_min and Ze_max.")
        bins = [b["Ze_min"] for b in bins] + [bins[-1]["Ze_max"]]
    bins = np.asarray(bins, dtype=float)

src/eventdisplay_ml/models.py:1374

• _classification_split_indices() splits groups independently per class using groups[label_mask].unique(). If the same group IDs appear in both classes but in a different order, the per-class train_test_split() calls can assign the same group to different splits across classes, breaking the intended group-disjoint property. Sorting label_groups (or performing one global split over the union of groups) makes the grouping deterministic and prevents cross-split leakage when group IDs overlap.

        for label in sorted(y_data.unique()):
            label_mask = y_data.to_numpy() == label
            label_groups = np.asarray(groups[label_mask].unique())
            n_train_groups = int(np.ceil(len(label_groups) * train_fraction))
            n_hold_groups = len(label_groups) - n_train_groups

• Files reviewed: 10/10 changed files
• Comments generated: 0 new
• Review effort level: Lite

We're testing this review assessment. Please use 👍 or 👎 to tell us if it's correct.

[Copilot]

🟡 Not ready to approve

Provenance/grouping is still keyed off file-path strings and the new per-row provenance index does not match the “original ROOT entry” contract described in the PR, which can undermine reproducibility/auditability.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

This review doesn't count toward merge requirements. Sign up for the private preview to control whether Copilot approvals count.

Review details

Suppressed comments (2)

src/eventdisplay_ml/models.py:1210

• Grouped classification splitting is currently keyed off __source_file (a path string). Since load_training_data() already creates a stable numeric __source_file_id, using that for groups avoids leaking absolute paths into split logic and makes splits reproducible even when file paths change (e.g., different mount points).

        full_df.get("__source_file"),

src/eventdisplay_ml/data_processing.py:1192

• __source_row is assigned as np.arange(len(file_df)), which is not the original ROOT entry index (it changes after cuts/sampling). The PR description calls for preserving the original entry number for provenance/duplicate detection; the current implementation can’t support that use case and may be misleading.

                    file_df["__source_row"] = np.arange(len(file_df), dtype=np.int64)

• Files reviewed: 10/10 changed files
• Comments generated: 0 new
• Review effort level: Lite

We're testing this review assessment. Please use 👍 or 👎 to tell us if it's correct.

[Copilot]

🟡 Not ready to approve

Telescope-configuration compatibility checks use exact float equality, which can spuriously fail across files due to minor floating-point representation differences and abort classification training.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

This review doesn't count toward merge requirements. Sign up for the private preview to control whether Copilot approvals count.

Review details

Suppressed comments (1)

src/eventdisplay_ml/data_processing.py:87

• _telescope_config_signature compares float arrays (mirror_area/tel_x/tel_y) via exact tuple equality. Small per-file floating-point representation differences (or NaNs) can make otherwise-identical telescope configurations appear different and abort classification training with a ValueError.

    """Return fields that determine the flattened classification schema."""
    return tuple(
        tuple(np.asarray(config[key]).tolist())
        for key in ("tel_ids", "mirror_area", "tel_x", "tel_y")
    )

• Files reviewed: 9/9 changed files
• Comments generated: 0 new
• Review effort level: Lite

We're testing this review assessment. Please use 👍 or 👎 to tell us if it's correct.

[Copilot]

🟡 Not ready to approve

Classification training applies zenith-balance sample weights but does not pass corresponding evaluation-set weights for early stopping/metrics, making weighted training inconsistent with unweighted validation.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

This review doesn't count toward merge requirements. Sign up for the private preview to control whether Copilot approvals count.

Review details

Suppressed comments (1)

src/eventdisplay_ml/models.py:1136

• When balance_class_zenith_weights is enabled, training passes sample_weight but does not pass corresponding weights for the eval_set. This means early-stopping and reported eval metrics are computed on unweighted validation data even though the training objective is weighted, which can select a different best iteration than intended and contradicts the PR description’s “validation weights” behavior.

        model = xgb.XGBClassifier(**cfg.get("hyper_parameters", {}))
        fit_kwargs = {"eval_set": eval_set, "verbose": True}
        if weights_train is not None:
            fit_kwargs["sample_weight"] = weights_train
        model.fit(x_train, y_train, **fit_kwargs)

• Files reviewed: 9/9 changed files
• Comments generated: 0 new
• Review effort level: Lite

We're testing this review assessment. Please use 👍 or 👎 to tell us if it's correct.