Feature Count Sweep

Martanto edited this page Jul 17, 2026 · 1 revision

Feature Count Sweep

⚠️ Experimental — not production-ready. The feature-count sweep is a post-hoc diagnostic added on the ft/feature-count-sweep branch. Its public surface (FeatureCountSweep, sweep_feature_count), defaults (mode="forecast", estimator_mode="default"), on-disk layout (sweep/{mode}/{classifier}/), and CLI wrappers may change without a deprecation cycle. Do not rely on the results as a substitute for a retrained model — treat the recommended N* as a hypothesis to verify with a fresh TrainingModel.fit(top_n_features=N*) run before adopting it.

Post-hoc "which top_n_features should I have picked?" utility. Reuses per-seed FDR rankings, resampled ids, and tuned hyperparameters that TrainingModel.fit() already wrote to disk, then re-scores each seed at several candidate N values.

Two entry points, both re-exported from eruption_forecast.features:

from eruption_forecast.features import FeatureCountSweep, sweep_feature_count

sweep_feature_count(source, ...) — cold-flow wrapper. Discovers the per-seed inputs from a completed training/ directory (or a live TrainingModel instance) and returns a fitted FeatureCountSweep (or a {classifier: sweep} dict).
FeatureCountSweep(...) — the underlying engine, callable directly when you already have (X, y, per-seed rankings) in memory.

Quick Start

Distilled from example_feature_count_sweep.py. Prerequisite: the scenario has already run ForecastModel.predict(...) (or PredictionModel.forecast(...) with save_model=True), so a {hash}.PredictionModel.pkl cache pickle exists.

import glob, os
from eruption_forecast import EvaluationModel
from eruption_forecast.features import sweep_feature_count
SCENARIO_DIR = "output/VG.OJN.00.EHZ/scenarios/scenario-1"
TRAINING_DIR = os.path.join(SCENARIO_DIR, "training")
ERUPTION_DATES = ["2025-04-10", "2025-04-22"]
# 1) Prepare an EvaluationModel in prediction-reuse mode. build_label()
# populates em.y_true cheaply (no per-classifier predict_proba pass).
prediction_pkl = glob.glob(os.path.join(SCENARIO_DIR, "prediction", "*.PredictionModel.pkl"))[0]
em = EvaluationModel.from_file(prediction_pkl, eruption_dates=ERUPTION_DATES)
em.build_label(window_step=em.window_step, window_step_unit=em.window_step_unit)
# 2) Sweep every classifier in the training directory.
results = sweep_feature_count(
 source=TRAINING_DIR,
 mode="forecast", # score against the prediction window (default)
 evaluation_source=em,
 n_candidates=list(range(1, 26)),
 scoring="roc_auc",
 save=True, # writes cv_scores.csv, curve.png, ...
)
# 3) Recommended N* per classifier.
for name, sweep in results.items():
 print(f"{name}: N* = {sweep.n_features_}")
# Reload later without re-running the sweep:
# from eruption_forecast.features import FeatureCountSweep
# sweep = FeatureCountSweep.load(
# f"{TRAINING_DIR}/features/{cv_slug}/sweep/forecast/RandomForestClassifier/FeatureCountSweep.pkl"
# )

Pass classifier_name="RandomForestClassifier" to focus on a single classifier (returns a FeatureCountSweep instead of a dict).

Fast within-training proxy — no evaluation_source needed:

results = sweep_feature_count(
 source=TRAINING_DIR,
 mode="cv", # RFECV-analog inside the resampled subset
 n_candidates=list(range(1, 26)),
 cv=5,
 scoring="roc_auc",
 n_jobs=4, # parallel over CV folds (cv mode only)
 save=True,
)

Two Scoring Modes

Mode	Scoring surface	Cost	When to use
`"forecast"` (default)	Held-out prediction-window features + ground-truth `y_true` from an `EvaluationModel` in prediction-reuse mode	One `fit()` per `(seed, N)`	Match the project's real forecasting evaluation
`"cv"`	`RFECV`-analog — cross-validation inside the trainer's resampled subset	`cv_splits` fits per `(seed, N)`	Fast proxy, but does not measure forecast-window generalisation

Motivation for forecast being the default: the CV proxy runs inside the same tiny under-sampled subset the trainer used, and on realistic eruption datasets it saturates at AP ≈ 1.0 across every N, producing no signal about temporal generalisation. Forecast mode answers "which N gives the best score on the held-out forecast window?" directly.

Fail-fast on missing features

In forecast mode, if any per-seed top-N list contains a feature that does not appear in the prediction feature matrix, the sweep raises KeyError immediately — it is never silently intersected or skipped. This surfaces the mismatch caused by PredictionModel.extract_features(select_tremor_columns=..., exclude_features=...) narrowing the extraction set. Fix by widening the prediction extraction (or by re-training with a narrower top_n_features pool).

`estimator_mode` — Default vs Tuned

Controls how each seed's best_model.pkl is prepared before every refit:

Value	Behaviour	Bias
`"default"` (recommended)	Reset the tuned estimator via `estimator.__class__()` — every candidate `N` is compared with the same untuned base learner	None (matches sklearn `RFECV` / yellowbrick convention)
`"tuned"`	Keep the trainer's `GridSearchCV`-picked hyperparameters	Biases the curve toward the trained `top_n_features` because those params were selected at that specific `N`

A large gap between the two curves is a signal that the trainer's hyperparameters are strongly N-specific — a smell worth investigating before adopting the recommended N*.

Picking `N*`

The recommendation is the N with the highest mean score across seeds, with an optional tie-break toward the smallest N:

parsimony=False — strict argmax on the mean.
parsimony=True, parsimony_tolerance=None (default) — 1-SE rule (adaptive): pick the smallest N whose mean is within one standard error of the peak.
parsimony=True, parsimony_tolerance=t — fractional tolerance: pick the smallest N whose mean is ≥ peak ×ばつ (1 − t).

For single-seed sweeps the standard error is undefined; the tolerance band collapses to zero and the peak N is returned.

The per-seed argmax N is exposed on sweep.seed_argmax_ — a diagnostic that reveals whether a shared N* across seeds is a reasonable assumption in the first place.

Forecast-Mode Wiring

sweep_feature_count(mode="forecast", ...) requires an EvaluationModel in prediction-reuse mode whose y_true has already been populated. The cheapest way to prepare one:

from eruption_forecast import EvaluationModel
em = EvaluationModel.from_file(
 "output/VG.OJN.00.EHZ/prediction/{hash}.PredictionModel.pkl",
 eruption_dates=["2025-04-10", "2025-04-22", "..."],
)
em.build_label(
 window_step=em.window_step,
 window_step_unit=em.window_step_unit,
)

build_label(...) populates em.y_true and writes y_true.csv under evaluation/prediction/labels/ without running the full per-classifier predict_proba pass that em.evaluate() would trigger. The sweep needs only em.features_df and em.y_true.

The _resolve_forecast_inputs guard validates the source and raises ValueError when:

evaluation_source is None
evaluation_source.model_kind != "prediction" (training-reuse labels are the training labels themselves and do not form a valid held-out set)
evaluation_source.y_true is empty (need to call build_label / evaluate first)
evaluation_source.features_df is empty (upstream PredictionModel never ran extract_features())

Complete Example

import glob
import os
from eruption_forecast import EvaluationModel
from eruption_forecast.features import sweep_feature_count
SCENARIO_DIR = "output/VG.OJN.00.EHZ/scenarios-new/scenario-11"
TRAINING_DIR = os.path.join(SCENARIO_DIR, "training")
PREDICTION_DIR = os.path.join(SCENARIO_DIR, "prediction")
ERUPTION_DATES = ["2025-04-10", "2025-04-22", "2025-05-18"]
# Auto-discover the persisted PredictionModel cache pickle.
prediction_pickle = sorted(
 glob.glob(os.path.join(PREDICTION_DIR, "*.PredictionModel.pkl"))
)[0]
em = EvaluationModel.from_file(prediction_pickle, eruption_dates=ERUPTION_DATES)
em.build_label(window_step=em.window_step, window_step_unit=em.window_step_unit)
results = sweep_feature_count(
 source=TRAINING_DIR,
 mode="forecast", # default
 evaluation_source=em,
 classifier_name=None, # sweep every classifier
 n_candidates=list(range(1, 26, 1)),
 scoring="roc_auc",
 parsimony=True,
 parsimony_tolerance=None, # → 1-SE rule
 estimator_mode="default", # RFECV convention
 random_state=42,
 save=True,
 verbose=True,
)
for classifier_name, sweep in results.items():
 print(f"[{classifier_name}] N* = {sweep.n_features_}")
 print(sweep.cv_scores_[["mean", "std"]].round(4).to_string())

The working script lives at example_feature_count_sweep.py on the branch. It ships run_sweep_forecast() (recommended default), run_sweep_cv() (fast proxy), and single-classifier variants.

Outputs

sweep_feature_count(save=True) writes under {training_dir}/features/{cv-slug}/sweep/{mode}/{classifier}/ — mode is part of the path so forecast/ and cv/ runs never collide:

×ばつ seed) score matrix ├── seed_argmax_hist.csv # per-seed argmax N (diagnostic) ├── support.json # {"n_features": N*, "seeds": {seed: [features...]}} ├── curve.png # two-panel plot: mean±std curve + argmax histogram └── FeatureCountSweep.pkl # full instance, reloadable via FeatureCountSweep.load(...)">

training/features/{cv-slug}/sweep/{mode}/{classifier-name}/
├── cv_scores.csv # aggregated summary: N, mean, std, n_seeds
├── cv_scores_raw.csv # full (N ×ばつ seed) score matrix
├── seed_argmax_hist.csv # per-seed argmax N (diagnostic)
├── support.json # {"n_features": N*, "seeds": {seed: [features...]}}
├── curve.png # two-panel plot: mean±std curve + argmax histogram
└── FeatureCountSweep.pkl # full instance, reloadable via FeatureCountSweep.load(...)

FeatureCountSweep.load(path) restores the full instance so you can re-plot or re-inspect without re-running the sweep.

Attributes After `fit()`

Attribute	Type	Notes
`sweep.cv_scores_`	`pd.DataFrame`	Indexed by `N`; columns `mean`, `std`, `n_seeds`
`sweep.cv_scores_raw_`	`pd.DataFrame`	Full `(N ×ばつ seed)` matrix — the source data for aggregation
`sweep.n_features_`	`int`	Recommended `N*`
`sweep.seed_argmax_`	`pd.Series`	Per-seed argmax `N`
`sweep.support_`	`dict[int, list[str]]`	For each seed, the top-`N*` feature list
`sweep.mode`	`str`	`"cv"` or `"forecast"` (echoed back for downstream code)
`sweep._classifier_name`	`str \| None`	Informational, set by the wrapper

Constructor / Function Reference