feat: rolling-origin event-aware evaluation harness for model variants

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

@@ -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())

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

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)