diff --git a/docs/FLOOD_FORECASTING.md b/docs/FLOOD_FORECASTING.md index 979f407..f2d1478 100644 --- a/docs/FLOOD_FORECASTING.md +++ b/docs/FLOOD_FORECASTING.md @@ -451,8 +451,9 @@ acts before any gauge rises, and `rain_fc24` — a weather *forecast* — acts before the rain itself falls. **Remaining honest limits:** marginal just-over-threshold crests (2025: +2 h) are intrinsically short-notice; the rain series only exists from 2021-03, so older training rows are rain-blind; -forecast-rain quality bounds what the feature can add; and Mae Ngat/Mae Kuang -dam releases remain uningested (see `docs/DATA_SOURCES.md`). +forecast-rain quality bounds what the feature can add; and Mae Ngat reservoir +state, though now ingested daily (see `docs/DATA_SOURCES.md`), measurably +*hurts* alert lead as a model feature — see the 2026-08-13 experiment below. **Danger-level skill at P.1 is unproven.** P.1 never crossed 4.5 m in the 2025-01-01 → 2026-08-10 test span (`base_rate_danger` is 0.0, so every danger @@ -475,6 +476,46 @@ P.1 additionally reports `stages`: exceedance probability for each of the seven official inundation stages (3.70–4.60 m), computed from the regression head and its calibration sigma, so they need no retrain and no per-stage classifiers. +### 2026-08-13: Mae Ngat dam features — a documented negative result + +With `rid_reservoir_daily` backfilled to 2018 (daily Mae Ngat storage/inflow/ +outflow, `src/ml/dam.py`), the obvious v4 experiment was to feed reservoir +state to the mainstem models: during the Oct 2024 flood the dam hit 113% of +usable capacity with 19–22 MCM/day inflow spikes on the crossing days. + +**It fails the acceptance gate.** On the 2024 record-flood backtest (train +< 1 Sep 2024, belt-and-braces alerting, identical to the deployed pipeline): + +| dam features | first-alert lead | record-peak err (24 h ahead) | +|----------------------------|------------------|------------------------------| +| none (deployed v3 config) | **+13 h** (PASS) | +0.24 m | +| all four | +10 h (FAIL) | +0.22 m | +| storage % + 3-day delta | +12 h | +0.21…+0.27 m | +| inflow + outflow | +10 h (FAIL) | +0.35 m | +| outflow only | +12 h | +0.20 m | + +Every subset costs 1–3 h of warning for at most a ~3 cm peak-error gain. The +mechanism is the publication lag: RID posts the daily report on the morning of +its own date (features apply it from 07:00, `dam.py`'s leakage rule), so at the +04:00 first-alert hour of 24 Sep 2024 the freshest dam row still described +23 Sep — a benign reservoir quietly absorbing inflow (outflow 0.13 MCM/day). +The columns therefore argue *against* imminent flooding exactly when the rain +features are (correctly) raising the alarm. The rolling-origin harness agrees: +`rise_rain_dam` matches `rise_rain` on leads and false alarms, only nudging +event-peak amplitude (−0.11 → −0.03 m on the Sep 2024 event), and `rise_dam` +(dam without rain) is strictly worse with alarm-latch artifacts. + +**Disposition:** dam features are OFF by default (`train_all(use_dam=False)`; +opt-in via `--dam` on the training CLI, `scripts/backtest_render.py --dam`, +and the `rise_rain_dam` / `rise_dam` harness variants). The collector keeps +accruing daily rows; revisit post-monsoon when the 2026 season adds dam-era +flood events — an intraday scrape (the lsim.rid.go.th source, reachable only +from Thai networks) would remove the publication-lag objection entirely. + +**Shipped from the same work:** the HII gap-fill merge in the data loader +(`fill_from_hii`, +9,341 h at P.81, +682 h at P.92, +810 h at P.20) is +lead-neutral — the gate holds at 13 h with fill on — and ships enabled. + ## 6. Deployment ### API diff --git a/docs/img/backtest-2024-p1-detail.png b/docs/img/backtest-2024-p1-detail.png index f960de3..95e0128 100644 Binary files a/docs/img/backtest-2024-p1-detail.png and b/docs/img/backtest-2024-p1-detail.png differ diff --git a/docs/img/backtest-2024-p1.png b/docs/img/backtest-2024-p1.png index 0052032..f6cf1eb 100644 Binary files a/docs/img/backtest-2024-p1.png and b/docs/img/backtest-2024-p1.png differ diff --git a/docs/img/backtest-2025-p1.png b/docs/img/backtest-2025-p1.png index 21225f1..4db948c 100644 Binary files a/docs/img/backtest-2025-p1.png and b/docs/img/backtest-2025-p1.png differ diff --git a/scripts/backtest_render.py b/scripts/backtest_render.py index 987243f..3c23d35 100644 --- a/scripts/backtest_render.py +++ b/scripts/backtest_render.py @@ -45,19 +45,23 @@ AMBER = "#c07d10" RED = "#d9534f" -def fit_backtest_model(df_long: pd.DataFrame, train_end: str): +def fit_backtest_model(df_long: pd.DataFrame, train_end: str, use_dam: bool = False): """Train the 24 h regression + warning heads on rows <= train_end only. Mirrors the deployed hgb-v3 pipeline: the regression head learns the RISE over the current level, with Open-Meteo catchment-rain features (trailing sums + the forward-24h forecast sum); label statistics are bounded to the - training cutoff. + training cutoff. use_dam=True adds the Mae Ngat reservoir columns — an + ablation-only configuration (2026-08-13 result: costs 1-3 h of lead). """ + from src.ml import dam as dam_mod from src.ml import rain as rain_mod rain_series = rain_mod.catchment_mean(rain_mod.load_history()) + dam_frame = dam_mod.load_history() if use_dam else None X, Y, _meta = features.build_matrix( - df_long, STATION, (HORIZON,), stats_end=train_end, rain=rain_series + df_long, STATION, (HORIZON,), stats_end=train_end, rain=rain_series, + dam=dam_frame, ) train_mask = X.index <= pd.Timestamp(train_end) X_train, Y_train = X.loc[train_mask], Y.loc[train_mask] @@ -200,16 +204,23 @@ def main(argv=None) -> int: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--db-url", default=None) parser.add_argument("--out-dir", default=os.path.join("docs", "img")) + parser.add_argument("--dam", action="store_true", + help="ablation: include Mae Ngat reservoir features " + "(2026-08 result: costs 1-3 h of alert lead)") + parser.add_argument("--no-hii-fill", action="store_true", + help="ablation: load without the HII gap-fill merge") args = parser.parse_args(argv) - df = data.load_measurements(db_url=args.db_url) + df = data.load_measurements( + db_url=args.db_url, hii_fill=not args.no_hii_fill + ) if df.empty: print("no measurement data available", file=sys.stderr) return 1 os.makedirs(args.out_dir, exist_ok=True) # --- October 2024 record flood: trained only on data before 1 Sep 2024 --- - X, reg, clf = fit_backtest_model(df, "2024-08-31") + X, reg, clf = fit_backtest_model(df, "2024-08-31", use_dam=args.dam) obs, fc, flood_start, first_alert = event_series( df, X, reg, clf, "2024-09-10", "2024-10-14 23:00") peak = float(obs.max()) @@ -235,7 +246,7 @@ def main(argv=None) -> int: detail=True) # --- September 2025 flood: the deployed configuration (trained <= 2024) --- - X25, reg25, clf25 = fit_backtest_model(df, "2024-12-31") + X25, reg25, clf25 = fit_backtest_model(df, "2024-12-31", use_dam=args.dam) obs25, fc25, flood25, alert25 = event_series( df, X25, reg25, clf25, "2025-09-22", "2025-10-02 12:00") pred_at_alert = float(fc25.loc[alert25:, "pred_max"].iloc[:24].max()) if alert25 is not None else None diff --git a/src/ml/dam.py b/src/ml/dam.py new file mode 100644 index 0000000..2eaef8b --- /dev/null +++ b/src/ml/dam.py @@ -0,0 +1,111 @@ +"""Mae Ngat reservoir series for the flood models. + +rid_reservoir_daily (collected hourly by src/rid_reservoir.py, backfilled to +2018) holds daily storage/inflow/outflow for every RID large dam. Mae Ngat +Somboon Chon (DAM_ID 200103) is the only large dam upstream of Chiang Mai: +in Oct 2024 its inflow hit 19-22 MCM/day and storage 114% of usable capacity +days around the P.1 crossing — upstream state no river gauge carries. + +Leakage rule: RID publishes the daily report for date D on the morning of D, +so the row becomes visible to features at D 07:00 local time, never earlier. +Forward-fill is capped at FFILL_LIMIT_H so a stalled collector degrades to +NaN (HGB-native) instead of silently serving stale reservoir state. + +Known residual optimism: the collector upserts keep-last (and re-fetches +yesterday), so the stored row for date D is RID's FINAL revision, which +training then back-dates to D 07:00 — values live serving may not have had +that morning. This bias works IN FAVOR of dam features, so the 2026-08-13 +negative result (they cost 1-3 h of alert lead) holds a fortiori; but any +future POSITIVE result must first validate intraday row stability or shift +the flow columns to D+1 07:00. +""" + +import datetime +import logging +from pathlib import Path +from typing import Optional + +import pandas as pd + +from ..rid_reservoir import MAE_NGAT_DAM_ID +from .data import CACHE_DIR, resolve_db_url + +logger = logging.getLogger(__name__) + +REPORT_HOUR = 7 # daily value valid from 07:00 local on its own date +FFILL_LIMIT_H = 48 # two missed daily reports -> NaN, not stale state +DAM_COLUMNS = ("storage_pct", "inflow_mcm", "outflow_mcm") +CACHE_FILE = f"dam_{MAE_NGAT_DAM_ID}.csv.gz" + + +def load_daily( + db_url: Optional[str] = None, + dam_id: str = MAE_NGAT_DAM_ID, + start: Optional[datetime.date] = None, + cache_dir: Path = CACHE_DIR, +) -> Optional[pd.DataFrame]: + """Daily dam rows indexed by date. DB first, on-disk cache as fallback.""" + cache_path = Path(cache_dir) / CACHE_FILE + resolved = resolve_db_url(db_url) + if resolved: + try: + from sqlalchemy import create_engine, text + + query = ( + "SELECT date, storage_pct, inflow_mcm, outflow_mcm " + "FROM rid_reservoir_daily WHERE dam_id = :dam_id" + ) + params = {"dam_id": dam_id} + if start is not None: + query += " AND date >= :start" + params["start"] = start + engine = create_engine(resolved, pool_pre_ping=True) + with engine.connect() as conn: + daily = pd.read_sql( + text(query + " ORDER BY date"), conn, params=params + ) + daily["date"] = pd.to_datetime(daily["date"]) + daily = daily.set_index("date") + for col in DAM_COLUMNS: + daily[col] = pd.to_numeric(daily[col], errors="coerce") + # Only full, NON-EMPTY loads refresh the cache: a truncated or + # freshly-recreated table must not wipe a good fallback archive. + if start is None and not daily.empty: + cache_path.parent.mkdir(parents=True, exist_ok=True) + daily.to_csv(cache_path, compression="gzip") + return daily + except Exception as error: + logger.warning(f"dam series DB load failed: {error}") + if cache_path.exists(): + logger.warning("falling back to on-disk cache for the dam series") + return pd.read_csv(cache_path, index_col=0, parse_dates=True) + return None + + +def hourly_frame(daily: Optional[pd.DataFrame]) -> Optional[pd.DataFrame]: + """Step the daily rows onto an hourly grid, each valid from D 07:00.""" + if daily is None or daily.empty: + return None + frame = daily.copy() + frame.index = pd.to_datetime(frame.index) + pd.Timedelta(hours=REPORT_HOUR) + frame = frame[~frame.index.duplicated(keep="last")].sort_index() + hourly_index = pd.date_range( + frame.index.min(), + frame.index.max() + pd.Timedelta(hours=FFILL_LIMIT_H), + freq="h", + ) + return frame.reindex(hourly_index).ffill(limit=FFILL_LIMIT_H) + + +def load_history(db_url: Optional[str] = None) -> Optional[pd.DataFrame]: + """Full hourly Mae Ngat history for training; None when unavailable.""" + return hourly_frame(load_daily(db_url)) + + +def serving_frame( + db_url: Optional[str] = None, days: int = 21 +) -> Optional[pd.DataFrame]: + """Recent hourly dam state for inference (covers the 336 h feature window + plus the 72 h storage-delta lag).""" + start = datetime.date.today() - datetime.timedelta(days=days) + return hourly_frame(load_daily(db_url, start=start)) diff --git a/src/ml/evaluate.py b/src/ml/evaluate.py index f13140a..57e9e2a 100644 --- a/src/ml/evaluate.py +++ b/src/ml/evaluate.py @@ -56,12 +56,14 @@ 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, use_rain: bool = False): + quantile: bool = False, use_rain: bool = False, + use_dam: bool = False): self.name = name self.target = target # 'abs' or 'rise' self.weighted = weighted self.quantile = quantile self.use_rain = use_rain + self.use_dam = use_dam def fit_predict( self, X_tr, y_abs_tr, X_te @@ -74,6 +76,14 @@ class Variant: raise ValueError( f"{self.name} requires the rain series (run without --no-rain)" ) + if not self.use_dam: + drop = [c for c in features.DAM_FEATURES if c in X_tr.columns] + X_tr = X_tr.drop(columns=drop) + X_te = X_te.drop(columns=drop) + elif "dam_storage_pct" not in X_tr.columns: + raise ValueError( + f"{self.name} requires the dam series (rid_reservoir_daily backfilled)" + ) 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 @@ -103,8 +113,16 @@ VARIANTS: Dict[str, Variant] = { "rise_quantile": Variant("rise_quantile", target="rise", weighted=True, quantile=True), "rise_rain": Variant("rise_rain", target="rise", use_rain=True), + "rise_rain_dam": Variant("rise_rain_dam", target="rise", use_rain=True, + use_dam=True), + "rise_dam": Variant("rise_dam", target="rise", use_dam=True), } +# Dam variants are opt-in by name: they require dam columns that only exist +# for features.DAM_STATIONS and only when the reservoir series loaded, and +# the 2026-08-13 ablation concluded them a negative result. +DEFAULT_VARIANTS = [k for k, v in VARIANTS.items() if not v.use_dam] + def _find_events(observed: pd.Series, thr: float) -> List[dict]: """Contiguous >=thr episodes (gaps under EVENT_GAP_H merged).""" @@ -198,11 +216,12 @@ def evaluate_station( variants: Optional[List[str]] = None, seasons: Tuple[int, ...] = SEASONS, rain: Optional[pd.Series] = None, + dam: Optional[pd.DataFrame] = None, ) -> 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, rain=rain) + X_all = features.build_features(grid, station, rain=rain, dam=dam) observed = grid.observed[(station, "water_level")] keep = X_all["obs_age_h"].notna() @@ -211,7 +230,7 @@ def evaluate_station( keep &= X_all.index >= pd.Timestamp(train_start) X_all = X_all.loc[keep] - chosen = {k: VARIANTS[k] for k in (variants or VARIANTS)} + chosen = {k: VARIANTS[k] for k in (variants or DEFAULT_VARIANTS)} results: Dict = {"station": station, "warn_thr": warn_thr, "folds": []} for year in seasons: @@ -253,7 +272,14 @@ def evaluate_station( } for name, variant in chosen.items(): - pred_abs, sigma = variant.fit_predict(X_tr, y_tr, X_te) + try: + pred_abs, sigma = variant.fit_predict(X_tr, y_tr, X_te) + except ValueError as error: + # A variant whose required feature family is absent (e.g. a + # dam variant on a non-DAM_STATIONS target) skips this fold + # instead of killing the whole run and its finished results. + logger.warning(f"{station} {year} {name}: skipped ({error})") + continue 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 @@ -352,6 +378,8 @@ def main(argv=None) -> int: parser.add_argument("--out", default="models/eval_variants.json") parser.add_argument("--no-rain", action="store_true", help="skip loading the Open-Meteo rain series") + parser.add_argument("--no-dam", action="store_true", + help="skip loading the Mae Ngat reservoir series") args = parser.parse_args(argv) logging.basicConfig( @@ -375,12 +403,27 @@ def main(argv=None) -> int: f"{rain_series.index.max()}" ) + dam_frame = None + if not args.no_dam: + from . import dam as dam_mod + + dam_frame = dam_mod.load_history(db_url=args.db_url) + if dam_frame is None: + logger.warning("dam history unavailable; dam features will be absent") + else: + logger.info( + f"dam series loaded: {dam_frame.index.min()} .. " + f"{dam_frame.index.max()}" + ) + 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, rain=rain_series) + results = evaluate_station( + df, station, variant_names, rain=rain_series, dam=dam_frame + ) all_results.append(results) print(summarize(results)) diff --git a/src/ml/features.py b/src/ml/features.py index 871171f..65d04a8 100644 --- a/src/ml/features.py +++ b/src/ml/features.py @@ -202,9 +202,18 @@ def _hours_since_observed(mask_col: pd.Series) -> pd.Series: RAIN_FEATURES = ("rain_6h", "rain_24h", "rain_72h", "rain_fc24") +DAM_FEATURES = ("dam_storage_pct", "dam_storage_pct_d3", "dam_inflow", "dam_outflow") +# Stations hydrologically downstream of the Mae Ngat confluence (Ping mainstem +# at/below Mae Taeng) — the only ones where reservoir state is causal. West- +# tributary and upper-mainstem stations never receive dam columns. +DAM_STATIONS = frozenset({"P.1", "P.103", "P.67", "P.21", "P.5", "P.81"}) + def build_features( - grid: HourlyGrid, station: str, rain: Optional[pd.Series] = None + grid: HourlyGrid, + station: str, + rain: Optional[pd.Series] = None, + dam: Optional[pd.DataFrame] = None, ) -> pd.DataFrame: """Build the deterministic-order feature matrix for one target station. @@ -217,6 +226,11 @@ def build_features( bundles even when the live fetch fails. rain_fc24 is the forward 24 h sum: the archived forecast series at training time, a real weather forecast at serving time; it never contains river data. + + ``dam`` is the hourly Mae Ngat reservoir frame (src/ml/dam.py; columns + storage_pct/inflow_mcm/outflow_mcm, already leakage-shifted to 07:00 + report time). Same contract as rain: None omits the columns, an empty + frame yields NaN columns; only DAM_STATIONS receive them. """ idx = grid.observed.index cols: Dict[str, pd.Series] = {} @@ -280,6 +294,18 @@ def build_features( r.shift(-1).iloc[::-1].rolling(24, min_periods=1).sum().iloc[::-1] ) + if dam is not None and station in DAM_STATIONS: + d = dam.reindex(idx) + + def _dam_col(name: str) -> pd.Series: + return d[name] if name in d.columns else pd.Series(np.nan, index=idx) + + storage = _dam_col("storage_pct") + cols["dam_storage_pct"] = storage + cols["dam_storage_pct_d3"] = storage - storage.shift(72) + cols["dam_inflow"] = _dam_col("inflow_mcm") + cols["dam_outflow"] = _dam_col("outflow_mcm") + return pd.DataFrame(cols, index=idx) @@ -358,10 +384,11 @@ def build_matrix( horizons: Tuple[int, ...] = (6, 12, 24), stats_end: Optional[str] = None, rain: Optional[pd.Series] = None, + dam: Optional[pd.DataFrame] = 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, rain=rain) + X = build_features(grid, station, rain=rain, dam=dam) Y = build_labels(grid, station, horizons, stats_end=stats_end) keep = X["obs_age_h"].notna() diff --git a/src/ml/predict.py b/src/ml/predict.py index e1582ee..5b9549a 100644 --- a/src/ml/predict.py +++ b/src/ml/predict.py @@ -128,6 +128,7 @@ def _model_forecast( as_of: pd.Timestamp, current_level: float, rain: Optional[pd.Series] = None, + dam: Optional[pd.DataFrame] = None, ) -> List[dict]: warn_thr = bundle["thresholds"]["warning"] danger_thr = bundle["thresholds"]["danger"] @@ -146,7 +147,9 @@ def _model_forecast( ) warn_thr, danger_thr = cfg_warn, cfg_danger - feature_row = features.build_features(grid, station_code, rain=rain).loc[[as_of]] + feature_row = features.build_features(grid, station_code, rain=rain, dam=dam).loc[ + [as_of] + ] expected_columns = bundle["feature_names"] missing = [c for c in expected_columns if c not in feature_row.columns] if missing: @@ -231,6 +234,7 @@ def _forecast_station( now: pd.Timestamp, horizons: Tuple[int, ...], rain: Optional[pd.Series] = None, + dam: Optional[pd.DataFrame] = None, ) -> List[dict]: level_col = (station_code, "water_level") if level_col not in grid.observed.columns: @@ -267,7 +271,7 @@ def _forecast_station( bundle = _load_bundle(bundle_path) model_results = _model_forecast( - station_code, grid, bundle, as_of, current_level, rain=rain + station_code, grid, bundle, as_of, current_level, rain=rain, dam=dam ) if model_results is None: return _heuristic_forecast( @@ -306,6 +310,7 @@ def get_forecasts( models_dir: Union[str, Path] = DEFAULT_MODELS_DIR, now: Optional[Union[datetime.datetime, str]] = None, rain: Optional[pd.Series] = None, + dam: Optional[pd.DataFrame] = None, ) -> List[dict]: """Produce flood forecasts for every station present in `readings_by_station`. @@ -329,7 +334,13 @@ def get_forecasts( try: results.extend( _forecast_station( - station_code, grid, models_dir, now, DEFAULT_HORIZONS, rain=rain + station_code, + grid, + models_dir, + now, + DEFAULT_HORIZONS, + rain=rain, + dam=dam, ) ) except Exception as error: @@ -376,4 +387,19 @@ def get_latest_forecasts( logger.warning("live rain unavailable; rain features will be NaN") rain = pd.Series(dtype=float) - return get_forecasts(readings_by_station, models_dir=models_dir, rain=rain) + # Recent Mae Ngat reservoir state; same empty-not-None contract so + # dam-trained bundles keep their columns (NaN) when the DB read fails. + from . import dam as dam_mod + + try: + dam = dam_mod.serving_frame(db_url=db_url) + except Exception as error: + logger.warning(f"dam serving frame failed: {error}") + dam = None + if dam is None: + logger.warning("dam state unavailable; dam features will be NaN") + dam = pd.DataFrame() + + return get_forecasts( + readings_by_station, models_dir=models_dir, rain=rain, dam=dam + ) diff --git a/src/ml/train.py b/src/ml/train.py index caf8aa7..682a97e 100644 --- a/src/ml/train.py +++ b/src/ml/train.py @@ -204,9 +204,10 @@ def train_station( split_test_start: str = SPLIT_B_TEST_START, split_test_end: str = SPLIT_B_TEST_END, rain: Optional[pd.Series] = None, + dam: Optional[pd.DataFrame] = None, ) -> Tuple[Optional[dict], dict]: """Train every head for one station. Returns (bundle_or_None, station_metrics).""" - X, Y, meta = features.build_matrix(df_long, station, horizons, rain=rain) + X, Y, meta = features.build_matrix(df_long, station, horizons, rain=rain, dam=dam) if meta["n_rows"] < MIN_ROWS_TO_TRAIN: return None, { "status": "failed", @@ -412,8 +413,13 @@ def train_station( ] = f"only {n_pos} positives in train span (< {MIN_POSITIVES_FOR_CLASSIFIER})" final_heads[head_key] = None - # v3 = rise target + Open-Meteo rain features; v2 = rise target only - version_prefix = "hgb-v3" if "rain_24h" in feature_names else "hgb-v2" + # v4 = + Mae Ngat dam features; v3 = rise + rain; v2 = rise target only + if "dam_storage_pct" in feature_names: + version_prefix = "hgb-v4" + elif "rain_24h" in feature_names: + version_prefix = "hgb-v3" + else: + version_prefix = "hgb-v2" bundle = { "station_code": station, "model_version": f"{version_prefix}+{_git_short_sha()}", @@ -443,6 +449,8 @@ def train_all( skip_eval: bool = False, hgb_overrides: Optional[dict] = None, use_rain: bool = True, + use_dam: bool = False, + db_url: Optional[str] = None, ) -> dict: """Train and save every requested station's models. Returns the metrics.json payload.""" models_dir = Path(models_dir) @@ -462,7 +470,42 @@ def train_all( logger.info( f"rain series: {rain_series.index.min()} .. {rain_series.index.max()}" ) - version_prefix = "hgb-v3" if rain_series is not None else "hgb-v2" + + # Mae Ngat reservoir state (rid_reservoir_daily, 2018+). OFF by default: + # the 2026-08-13 backtest ablation showed every dam-feature subset COSTS + # 1-3 h of first-alert lead on the 2024 record flood (the daily report + # lags up to 31 h, so during fast onset the columns describe yesterday's + # benign reservoir and damp the alarm). Kept as an opt-in for post-monsoon + # re-evaluation once the 2026 season adds dam-era flood events. + dam_frame = None + if use_dam: + try: + from . import dam as dam_mod + + dam_frame = dam_mod.load_history(db_url=db_url) + except Exception as error: + logger.warning(f"dam history unavailable, training without it: {error}") + if dam_frame is None: + # load_history returns None (no raise) when both DB and cache + # miss — an explicitly requested experiment must say so loudly. + logger.warning( + "--dam requested but no dam history available; " + "training v3-style bundles WITHOUT dam features" + ) + if dam_frame is not None: + logger.info( + f"dam series: {dam_frame.index.min()} .. {dam_frame.index.max()}" + ) + + # Run-level version: v4 only if some requested station actually receives + # dam columns (they are gated to DAM_STATIONS; per-bundle versions are + # derived from each station's own feature_names and remain authoritative). + if dam_frame is not None and any(s in features.DAM_STATIONS for s in stations): + version_prefix = "hgb-v4" + elif rain_series is not None: + version_prefix = "hgb-v3" + else: + version_prefix = "hgb-v2" model_version = f"{version_prefix}+{_git_short_sha()}" station_results: Dict[str, dict] = {} @@ -480,6 +523,7 @@ def train_all( skip_eval=skip_eval, hgb_overrides=hgb_overrides, rain=rain_series, + dam=dam_frame, ) if bundle is None: logger.warning(f"{station}: failed ({station_metrics.get('reason')})") @@ -542,6 +586,12 @@ def main(argv: Optional[List[str]] = None) -> None: action="store_true", help="train without the Open-Meteo rain features (v2-style bundles)", ) + parser.add_argument( + "--dam", + action="store_true", + help="EXPERIMENTAL: include Mae Ngat reservoir features (v4 bundles); " + "the 2026-08 ablation showed they cost 1-3 h of alert lead", + ) args = parser.parse_args(argv) if args.stations == "all": @@ -570,6 +620,8 @@ def main(argv: Optional[List[str]] = None) -> None: models_dir=Path(args.models_dir), skip_eval=args.skip_eval, use_rain=not args.no_rain, + use_dam=args.dam, + db_url=resolve_db_url(args.db_url), ) trained = sum( 1 for s in metrics_payload["stations"].values() if s["status"] == "trained" diff --git a/tests/test_dam_features.py b/tests/test_dam_features.py new file mode 100644 index 0000000..9aaa3bd --- /dev/null +++ b/tests/test_dam_features.py @@ -0,0 +1,128 @@ +"""Tests for the Mae Ngat dam series (src/ml/dam.py) and its feature gating.""" + +import datetime + +import numpy as np +import pandas as pd + +from src.ml import features +from src.ml.dam import FFILL_LIMIT_H, REPORT_HOUR, hourly_frame + + +def _daily(days=5, start="2024-09-20"): + idx = pd.date_range(start, periods=days, freq="D") + return pd.DataFrame( + { + "storage_pct": np.linspace(90, 110, days), + "inflow_mcm": np.linspace(2, 20, days), + "outflow_mcm": np.linspace(0.5, 5, days), + }, + index=idx, + ) + + +class TestHourlyFrame: + def test_daily_value_visible_from_report_hour_only(self): + hourly = hourly_frame(_daily()) + day0 = pd.Timestamp("2024-09-20") + # Nothing before the first report hour + assert hourly.index.min() == day0 + pd.Timedelta(hours=REPORT_HOUR) + # The day's value holds from 07:00 through the next morning + assert hourly.loc[day0 + pd.Timedelta(hours=7), "storage_pct"] == 90.0 + assert hourly.loc[day0 + pd.Timedelta(hours=23), "storage_pct"] == 90.0 + next_6am = day0 + pd.Timedelta(days=1, hours=6) + next_7am = day0 + pd.Timedelta(days=1, hours=7) + assert hourly.loc[next_6am, "storage_pct"] == 90.0 # yesterday's value + assert hourly.loc[next_7am, "storage_pct"] == 95.0 # today's report + + def test_ffill_capped_after_missing_days(self): + daily = _daily(days=2).drop(index=pd.Timestamp("2024-09-21")) + # extend with a far-later row so the gap sits mid-frame + late = _daily(days=1, start="2024-09-28") + hourly = hourly_frame(pd.concat([daily, late])) + gap_ts = pd.Timestamp("2024-09-20") + pd.Timedelta( + hours=REPORT_HOUR + FFILL_LIMIT_H + 1 + ) + assert np.isnan(hourly.loc[gap_ts, "storage_pct"]) + + def test_none_and_empty(self): + assert hourly_frame(None) is None + assert hourly_frame(pd.DataFrame()) is None + + +class TestLoadDaily: + def test_empty_db_result_does_not_wipe_cache(self, tmp_path): + from sqlalchemy import create_engine, text + + from src.ml.dam import CACHE_FILE, load_daily + + # Good cache from a previous run + cache_path = tmp_path / CACHE_FILE + _daily(3).rename_axis("date").to_csv(cache_path, compression="gzip") + # Reachable DB whose table exists but is empty + db = f"sqlite:///{tmp_path}/empty_dam.db" + with create_engine(db).begin() as conn: + conn.execute( + text( + "CREATE TABLE rid_reservoir_daily (dam_id TEXT, date DATE, " + "storage_pct REAL, inflow_mcm REAL, outflow_mcm REAL)" + ) + ) + result = load_daily(db_url=db, cache_dir=tmp_path) + assert result.empty # honest empty result... + cached = pd.read_csv(cache_path, index_col=0) + assert len(cached) == 3 # ...but the good cache survives + + +def _grid(hours=400, start="2024-09-15"): + idx = pd.date_range(start, periods=hours, freq="h") + frames = [] + for code in ("P.1", "P.82"): + frames.append( + pd.DataFrame( + { + "timestamp": idx, + "station_code": code, + "water_level": 2.0, + "discharge": 100.0, + } + ) + ) + return features.make_hourly_grid(pd.concat(frames, ignore_index=True)) + + +class TestFeatureGating: + def test_dam_columns_only_for_dam_stations(self): + grid = _grid() + dam = hourly_frame(_daily(days=20, start="2024-09-10")) + X_p1 = features.build_features(grid, "P.1", dam=dam) + X_p82 = features.build_features(grid, "P.82", dam=dam) + for col in features.DAM_FEATURES: + assert col in X_p1.columns + assert col not in X_p82.columns + # values actually aligned, not all-NaN + assert X_p1["dam_storage_pct"].notna().any() + assert X_p1["dam_inflow"].notna().any() + + def test_none_dam_omits_columns(self): + X = features.build_features(_grid(), "P.1", dam=None) + for col in features.DAM_FEATURES: + assert col not in X.columns + + def test_empty_dam_frame_yields_nan_columns(self): + # Serving contract: empty frame -> columns exist as NaN so dam-trained + # bundles pass the feature guard when the DB read fails. + X = features.build_features(_grid(), "P.1", dam=pd.DataFrame()) + for col in features.DAM_FEATURES: + assert col in X.columns + assert X[col].isna().all() + + def test_storage_delta_72h(self): + grid = _grid(hours=24 * 12, start="2024-09-15") + dam = hourly_frame(_daily(days=20, start="2024-09-10")) + X = features.build_features(grid, "P.1", dam=dam) + ts = pd.Timestamp("2024-09-24 12:00") + expected = X.loc[ts, "dam_storage_pct"] - X.loc[ + ts - pd.Timedelta(hours=72), "dam_storage_pct" + ] + assert abs(X.loc[ts, "dam_storage_pct_d3"] - expected) < 1e-9 diff --git a/tests/test_flood_forecast.py b/tests/test_flood_forecast.py index 6e66cfa..4aa1aee 100644 --- a/tests/test_flood_forecast.py +++ b/tests/test_flood_forecast.py @@ -197,7 +197,7 @@ def test_train_smoke_and_roundtrip(tmp_path): df = make_synth(n, data_stations, seed=7, pulses=pulses) metrics = train.train_all( - df, target_stations, models_dir=tmp_path, skip_eval=True, hgb_overrides={"max_iter": 20}, use_rain=False + df, target_stations, models_dir=tmp_path, skip_eval=True, hgb_overrides={"max_iter": 20}, use_rain=False, use_dam=False ) assert metrics["stations"]["P.1"]["status"] == "trained" assert metrics["stations"]["P.20"]["status"] == "trained" @@ -237,7 +237,7 @@ def test_feature_name_stability(tmp_path): upstream = [code for code, _lead in features.UPSTREAM_LEADS["P.1"]] data_stations = ["P.1"] + upstream df = make_synth(300, data_stations, seed=11, pulses={"P.1": [(100, 20, 2.0)]}) - train.train_all(df, ["P.1"], models_dir=tmp_path, skip_eval=True, hgb_overrides={"max_iter": 10}, use_rain=False) + train.train_all(df, ["P.1"], models_dir=tmp_path, skip_eval=True, hgb_overrides={"max_iter": 10}, use_rain=False, use_dam=False) # Safe: loading the bundle this same test just wrote to tmp_path, not an external file. bundle = joblib.load(tmp_path / "flood_P.1.joblib")