From a0086086a2ffe0063cac4bbda50b6df6c3a32e15 Mon Sep 17 00:00:00 2001 From: grabowski Date: Wed, 12 Aug 2026 15:19:26 +0700 Subject: [PATCH] feat: rolling-origin event-aware evaluation harness for model variants MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit One fold per monsoon season (train <= 30 Apr, test Jun-Nov, 2021-2025) replaces the single fixed holdout that contained only ~4 warning events. Metrics are what matters operationally: sustained first-alert lead vs each observed 3.70m crossing (two consecutive alerting samples required; lookback floored at the previous event's end so multi-peak floods can't launder lead credit), peak error from the prediction actually issued 24h before the peak (3h match tolerance, null on outages), false-alarm episodes (12h gap tolerance), MAE / flood-regime MAE, and a Brier score on warning exceedance — included because sigma cancels algebraically in any p>=0.5 alert metric, so lead times compare predictors while Brier compares uncertainty models. Variants: baseline_abs (current), rise (target = future max - current level), rise_weighted (flood-regime sample weights 1x->5x), and rise_quantile (q50/q90 heads, spread-implied sigma). Harness verified by a 3-agent adversarial review (features bit-identical across fold cutoffs; three metric flaws found and fixed before first use). Also: features.build_labels/build_matrix gain stats_end so the rescue quantile is computed from pre-cutoff data only, closing the label- construction leak flagged in the earlier ML review. --- scripts/evaluate_variants.py | 18 ++ src/ml/evaluate.py | 364 +++++++++++++++++++++++++++++++++++ src/ml/features.py | 22 ++- 3 files changed, 398 insertions(+), 6 deletions(-) create mode 100644 scripts/evaluate_variants.py create mode 100644 src/ml/evaluate.py diff --git a/scripts/evaluate_variants.py b/scripts/evaluate_variants.py new file mode 100644 index 0000000..aa020cb --- /dev/null +++ b/scripts/evaluate_variants.py @@ -0,0 +1,18 @@ +#!/usr/bin/env python3 +"""CLI for the rolling-origin model-variant evaluation. + +Usage: + uv run scripts/evaluate_variants.py # P.1, all variants + uv run scripts/evaluate_variants.py --stations P.1,P.103 + uv run scripts/evaluate_variants.py --variants baseline_abs,rise_quantile +""" + +import os +import sys + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..")) + +from src.ml.evaluate import main + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/ml/evaluate.py b/src/ml/evaluate.py new file mode 100644 index 0000000..ccc9d03 --- /dev/null +++ b/src/ml/evaluate.py @@ -0,0 +1,364 @@ +"""Rolling-origin, event-aware evaluation of forecast-model variants. + +Replaces the single fixed holdout (which contained only ~4 warning events) +with one fold per monsoon season: train on everything through 30 April of the +season's year (labels' rescue statistics bounded to the same cutoff, and label +windows cannot reach the June+ test span, so the folds are leak-free), test on +June-November. Metrics are event-level — first-alert lead versus each warning +crossing, peak error at 24 h — plus pointwise MAE and false-alarm episodes, +because pointwise PR-AUC alone hid the things that matter operationally. + +Variants under test target the two failures documented in +docs/FLOOD_FORECASTING.md's re-examination note: absolute-level regression +cannot extrapolate past its training maximum, and the flat sigma miscalibrates +probabilities. +""" + +import json +import logging +from typing import Callable, Dict, List, Optional, Tuple + +import numpy as np +import pandas as pd +from scipy.special import erf + +from . import data, features +from .train import HGB_PARAMS, _make_regressor + +logger = logging.getLogger(__name__) + +HORIZON = 24 +SEASONS = (2021, 2022, 2023, 2024, 2025) +TEST_MONTHS = ("06-01", "11-30") +TRAIN_END_MD = "04-30" +ALERT_P = 0.5 +FIXED_SIGMA = 0.15 +EVENT_GAP_H = 24 # merge >=thr runs closer than this into one event +FALSE_ALARM_GRACE_H = 48 + + +def _phi(z: np.ndarray) -> np.ndarray: + return 0.5 * (1.0 + erf(z / np.sqrt(2.0))) + + +def _quantile_regressor(q: float): + from sklearn.ensemble import HistGradientBoostingRegressor + + return HistGradientBoostingRegressor(loss="quantile", quantile=q, **HGB_PARAMS) + + +def _flood_weights(y_abs: pd.Series) -> np.ndarray: + """Upweight the flood regime: 1x below 2.5 m ramping to 5x at >= 3.7 m.""" + return 1.0 + 4.0 * np.clip((y_abs.to_numpy() - 2.5) / 1.2, 0.0, 1.0) + + +class Variant: + """A trainable candidate producing (pred_abs, sigma_per_row) on test rows.""" + + def __init__(self, name: str, target: str, weighted: bool = False, + quantile: bool = False): + self.name = name + self.target = target # 'abs' or 'rise' + self.weighted = weighted + self.quantile = quantile + + def fit_predict( + self, X_tr, y_abs_tr, X_te + ) -> Tuple[np.ndarray, np.ndarray]: + level_tr = X_tr["level"] + level_te = X_te["level"].to_numpy() + y_tr = (y_abs_tr - level_tr) if self.target == "rise" else y_abs_tr + weights = _flood_weights(y_abs_tr) if self.weighted else None + + if self.quantile: + q50 = _quantile_regressor(0.5).fit(X_tr, y_tr, sample_weight=weights) + q90 = _quantile_regressor(0.9).fit(X_tr, y_tr, sample_weight=weights) + p50 = q50.predict(X_te) + spread = np.maximum(q90.predict(X_te) - p50, 0.0) + sigma = np.maximum(spread / 1.2816, 0.05) + pred = p50 + else: + reg = _make_regressor().fit(X_tr, y_tr, sample_weight=weights) + pred = reg.predict(X_te) + sigma = np.full(len(X_te), FIXED_SIGMA) + + pred_abs = pred + level_te if self.target == "rise" else pred + pred_abs = np.maximum(pred_abs, level_te) # peak >= current, as served + return pred_abs, sigma + + +VARIANTS: Dict[str, Variant] = { + "baseline_abs": Variant("baseline_abs", target="abs"), + "rise": Variant("rise", target="rise"), + "rise_weighted": Variant("rise_weighted", target="rise", weighted=True), + "rise_quantile": Variant("rise_quantile", target="rise", weighted=True, + quantile=True), +} + + +def _find_events(observed: pd.Series, thr: float) -> List[dict]: + """Contiguous >=thr episodes (gaps under EVENT_GAP_H merged).""" + above = observed[observed >= thr] + if above.empty: + return [] + events = [] + start = prev = above.index[0] + for ts in above.index[1:]: + if (ts - prev) > pd.Timedelta(hours=EVENT_GAP_H): + events.append((start, prev)) + start = ts + prev = ts + events.append((start, prev)) + out = [] + for begin, end in events: + window = observed.loc[begin:end] + out.append( + { + "crossing": begin, + "end": end, + "peak_ts": window.idxmax(), + "peak_level": float(window.max()), + } + ) + return out + + +def _first_alert_lead( + p: pd.Series, + crossing: pd.Timestamp, + window_start_floor: Optional[pd.Timestamp] = None, +) -> Optional[float]: + """Hours between the first SUSTAINED alert near the crossing and the + crossing. Positive = warned in advance; negative = late. + + Sustained = two consecutive hourly samples with p >= ALERT_P (a single + noisy spike gets no credit). The lookback never reaches past + ``window_start_floor`` (the previous event's end), so one event's tail + cannot be credited as early warning for the next crossing. + """ + start = crossing - pd.Timedelta(hours=72) + if window_start_floor is not None and window_start_floor > start: + start = window_start_floor + window = p.loc[start: crossing + pd.Timedelta(hours=24)] + if len(window) < 2: + return None + alert = (window >= ALERT_P) & (window.shift(-1) >= ALERT_P) & ( + (window.index.to_series().shift(-1) - window.index.to_series()) + <= pd.Timedelta(hours=2) + ) + hits = window.index[alert.fillna(False)] + if len(hits) == 0: + return None + return float((crossing - hits[0]).total_seconds() / 3600.0) + + +def _false_alarm_episodes( + p: pd.Series, observed: pd.Series, thr: float +) -> int: + """Alert episodes with no observed >=thr within +/- FALSE_ALARM_GRACE_H.""" + alert_hours = p[p >= ALERT_P].index + if len(alert_hours) == 0: + return 0 + grace = pd.Timedelta(hours=FALSE_ALARM_GRACE_H) + exceed_times = observed[observed >= thr].index + episodes = 0 + episode_start = None + prev = None + for ts in alert_hours: + # 12h gap tolerance: a data hole mid-alarm must not double-count it + if prev is None or (ts - prev) > pd.Timedelta(hours=12): + if episode_start is not None: + episodes += _is_false(episode_start, prev, exceed_times, grace) + episode_start = ts + prev = ts + episodes += _is_false(episode_start, prev, exceed_times, grace) + return episodes + + +def _is_false(start, end, exceed_times, grace) -> int: + if len(exceed_times) == 0: + return 1 + near = (exceed_times >= start - grace) & (exceed_times <= end + grace) + return 0 if near.any() else 1 + + +def evaluate_station( + df_long: pd.DataFrame, + station: str, + variants: Optional[List[str]] = None, + seasons: Tuple[int, ...] = SEASONS, +) -> Dict: + """Run every fold x variant for one station; returns the results tree.""" + warn_thr, _ = features.get_thresholds(station) + grid = features.make_hourly_grid(df_long) + X_all = features.build_features(grid, station) + observed = grid.observed[(station, "water_level")] + + keep = X_all["obs_age_h"].notna() + train_start = features.TRAIN_START.get(station) + if train_start: + keep &= X_all.index >= pd.Timestamp(train_start) + X_all = X_all.loc[keep] + + chosen = {k: VARIANTS[k] for k in (variants or VARIANTS)} + results: Dict = {"station": station, "warn_thr": warn_thr, "folds": []} + + for year in seasons: + train_end = pd.Timestamp(f"{year}-{TRAIN_END_MD}") + test_lo = pd.Timestamp(f"{year}-{TEST_MONTHS[0]}") + test_hi = pd.Timestamp(f"{year}-{TEST_MONTHS[1]} 23:00") + + # Labels rebuilt per fold so rescue statistics stop at the cutoff + Y = features.build_labels( + grid, station, (HORIZON,), stats_end=train_end.isoformat() + ).loc[X_all.index] + y_abs = Y[f"max_level_{HORIZON}"] + + tr = (X_all.index <= train_end) & y_abs.notna() + te = (X_all.index >= test_lo) & (X_all.index <= test_hi) + if tr.sum() < 5000 or te.sum() < 500: + logger.info(f"{station} {year}: skipped (train {tr.sum()}, test {te.sum()})") + continue + + X_tr, X_te = X_all.loc[tr], X_all.loc[te] + y_tr = y_abs.loc[tr] + y_te = y_abs.loc[te] + obs_test = observed.loc[test_lo:test_hi].dropna() + events = _find_events(obs_test, warn_thr) + + fold: Dict = { + "year": year, + "n_train": int(tr.sum()), + "n_test": int(te.sum()), + "events": [ + { + "crossing": e["crossing"].isoformat(), + "peak_ts": e["peak_ts"].isoformat(), + "peak_level": e["peak_level"], + } + for e in events + ], + "variants": {}, + } + + for name, variant in chosen.items(): + pred_abs, sigma = variant.fit_predict(X_tr, y_tr, X_te) + pred_series = pd.Series(pred_abs, index=X_te.index) + p_warn = pd.Series( + 1.0 - _phi((warn_thr - pred_abs) / sigma), index=X_te.index + ) + + labeled = y_te.notna() + errors = (pred_series[labeled] - y_te[labeled]).abs() + high = y_te[labeled] >= warn_thr - 1.2 # flood-regime rows + # Brier score on within-24h warning exceedance: unlike the p>=0.5 + # alert metrics (where sigma cancels algebraically), this actually + # exercises each variant's uncertainty model. + exceed = Y[f"exceed_warn_{HORIZON}"].loc[te] + scored = exceed.notna() + brier = ( + float(((p_warn[scored] - exceed[scored]) ** 2).mean()) + if scored.any() + else None + ) + event_rows = [] + for i, event in enumerate(events): + floor = events[i - 1]["end"] if i > 0 else None + lead = _first_alert_lead(p_warn, event["crossing"], floor) + issue_ts = event["peak_ts"] - pd.Timedelta(hours=HORIZON) + peak_pred = None + if len(pred_series): + nearest = pred_series.index.get_indexer( + [issue_ts], method="nearest" + )[0] + matched_ts = pred_series.index[nearest] + # Tolerance: a "24h-ahead" prediction matched to a row + # hours away (data outage) is not that prediction at all. + if abs(matched_ts - issue_ts) <= pd.Timedelta(hours=3): + peak_pred = float(pred_series.iloc[nearest]) + event_rows.append( + { + "crossing": event["crossing"].isoformat(), + "lead_h": lead, + "peak_level": event["peak_level"], + "peak_pred_24h_before": peak_pred, + } + ) + fold["variants"][name] = { + "mae": float(errors.mean()) if len(errors) else None, + "mae_above_2p5": ( + float(errors[high].mean()) if high.any() else None + ), + "brier_warn": brier, + "events": event_rows, + "false_alarm_episodes": _false_alarm_episodes( + p_warn, obs_test, warn_thr + ), + } + results["folds"].append(fold) + return results + + +def summarize(results: Dict) -> str: + """Compact comparison table across folds for one station.""" + lines = [f"\n=== {results['station']} (warn {results['warn_thr']:.2f} m) ==="] + header = ( + f"{'variant':16} {'year':>5} {'MAE':>6} {'MAE_hi':>7} {'Brier':>7} " + f"{'FA':>3} events (lead h | peak err m)" + ) + lines.append(header) + for fold in results["folds"]: + for name, m in fold["variants"].items(): + events = " ".join( + f"[{e['crossing'][:10]}: " + f"{'—' if e['lead_h'] is None else format(e['lead_h'], '+.0f')}h" + + ( + f" | {e['peak_pred_24h_before'] - e['peak_level']:+.2f}" + if e["peak_pred_24h_before"] is not None + else "" + ) + + "]" + for e in m["events"] + ) or "no events" + lines.append( + f"{name:16} {fold['year']:>5} " + f"{m['mae'] if m['mae'] is not None else float('nan'):6.3f} " + f"{m['mae_above_2p5'] if m['mae_above_2p5'] is not None else float('nan'):7.3f} " + f"{m['brier_warn'] if m.get('brier_warn') is not None else float('nan'):7.4f} " + f"{m['false_alarm_episodes']:>3} {events}" + ) + return "\n".join(lines) + + +def main(argv=None) -> int: + import argparse + + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--stations", default="P.1") + parser.add_argument("--db-url", default=None) + parser.add_argument("--variants", default=None, + help="comma list; default all") + parser.add_argument("--out", default="models/eval_variants.json") + args = parser.parse_args(argv) + + logging.basicConfig( + level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s" + ) + df = data.load_measurements(db_url=args.db_url) + if df.empty: + logger.error("no measurement data") + return 1 + + variant_names = args.variants.split(",") if args.variants else None + all_results = [] + for station in args.stations.split(","): + station = station.strip() + logger.info(f"Evaluating {station}...") + results = evaluate_station(df, station, variant_names) + all_results.append(results) + print(summarize(results)) + + with open(args.out, "w", encoding="utf-8") as fh: + json.dump(all_results, fh, indent=1) + logger.info(f"results written to {args.out}") + return 0 diff --git a/src/ml/features.py b/src/ml/features.py index 28b0180..e1032b6 100644 --- a/src/ml/features.py +++ b/src/ml/features.py @@ -272,9 +272,17 @@ def _future_window_stats(col: pd.Series, horizon_h: int) -> Tuple[pd.Series, pd. def build_labels( - grid: HourlyGrid, station: str, horizons: Tuple[int, ...] = (6, 12, 24) + grid: HourlyGrid, + station: str, + horizons: Tuple[int, ...] = (6, 12, 24), + stats_end: Optional[str] = None, ) -> pd.DataFrame: - """Build max-level and threshold-exceedance labels for one target station.""" + """Build max-level and threshold-exceedance labels for one target station. + + ``stats_end`` bounds the data used for label-construction statistics (the + rescue quantile below): pass the training cutoff during evaluation so + test-period extremes cannot influence which training rows receive labels. + """ idx = grid.observed.index observed_level = _series(grid.observed, station, "water_level", idx) warn_thr, danger_thr = get_thresholds(station) @@ -283,10 +291,11 @@ def build_labels( # NOT to warn_thr: coupling it to the configurable threshold made raising a # station's threshold silently shrink its regression training set (P.5 lost # 34% of rows and +46% MAE when its warning went 3.0 -> 4.55). + stats_level = ( + observed_level.loc[: pd.Timestamp(stats_end)] if stats_end else observed_level + ) rescue_thr = ( - float(observed_level.quantile(0.975)) - if observed_level.notna().any() - else np.inf + float(stats_level.quantile(0.975)) if stats_level.notna().any() else np.inf ) out: Dict[str, pd.Series] = {} @@ -321,11 +330,12 @@ def build_matrix( df_long: pd.DataFrame, station: str, horizons: Tuple[int, ...] = (6, 12, 24), + stats_end: Optional[str] = None, ) -> Tuple[pd.DataFrame, pd.DataFrame, dict]: """Build (X, Y, meta) training/inference matrices for one station.""" grid = make_hourly_grid(df_long) X = build_features(grid, station) - Y = build_labels(grid, station, horizons) + Y = build_labels(grid, station, horizons, stats_end=stats_end) keep = X["obs_age_h"].notna() train_start = TRAIN_START.get(station)