From 4358d52d5583b9af3dcb323dd41b77367c2428cc Mon Sep 17 00:00:00 2001 From: grabowski Date: Mon, 10 Aug 2026 12:49:47 +0700 Subject: [PATCH] feat: ML flood-event forecasting from 8 years of gauge history Add src/ml/ package predicting, per station and per 6/12/24 h horizon, the probability of exceeding warning (3.0 m) and danger (4.5 m) levels plus expected peak level, trained on the 592k-row PostgreSQL history: - features.py: hourly grid with coverage gating and no future leakage; upstream stations enter at empirically measured travel-time lags (P.20 +17h ... P.103 +1h vs P.1); hour-of-day deliberately excluded (it encodes the scrape schedule, not hydrology) - train.py: HistGradientBoosting regression + warn/danger classifier heads per station x horizon, >=30-positives gate with calibrated sigmoid-on-regression fallback, strict temporal splits, per-event lead-time evaluation; guards against sklearn 1.9.0 crash on degenerate feature columns - predict.py: bundle loading with feature-name checks, heuristic fallback tier, get_latest_forecasts() for the API; raises when no models are trained so the endpoint 503s instead of serving persistence output as forecasts - data.py: Postgres-first loader (FLOOD_ML_DB_URL override), HTTP API fallback (flagged: that path backfills synthetic discharge), csv.gz cache - /forecast endpoint (15-min TTL cache) + dashboard flood-risk panel (hidden until models exist) - docs/FLOOD_FORECASTING.md: full system doc with measured deployment numbers (~335 MB RSS, CPU negligible, ~6 min full retrain) and retraining policy Validation: out-of-sample backtest of the record 2024 flood season (train <= Aug 2024) alerted 24-48 h ahead of the Oct 5 peak; 2025-26 test split: P.1 6h PR-AUC 0.974, recall 98.3% at 1% false-alarm rate. Also: fix P.81 station coordinates (was Ban Pong/Ratchaburi, 493 km out of basin; now 18.6936 N 99.0819 E per RID station page), pin scikit-learn==1.9.0 and numpy<2, gitignore model artifacts (~100 MB, train on the server via scripts/train_flood_model.py). --- .gitignore | 11 + docs/FLOOD_FORECASTING.md | 621 +++++++++++++++++++++++++++++++++++ models/.gitkeep | 0 pyproject.toml | 3 + requirements.txt | 4 + scripts/train_flood_model.py | 17 + src/data/stations.json | 4 +- src/ml/__init__.py | 1 + src/ml/data.py | 223 +++++++++++++ src/ml/features.py | 272 +++++++++++++++ src/ml/predict.py | 284 ++++++++++++++++ src/ml/train.py | 414 +++++++++++++++++++++++ src/static/dashboard.html | 67 ++++ src/web_api.py | 31 ++ tests/test_flood_forecast.py | 238 ++++++++++++++ 15 files changed, 2188 insertions(+), 2 deletions(-) create mode 100644 docs/FLOOD_FORECASTING.md create mode 100644 models/.gitkeep create mode 100644 scripts/train_flood_model.py create mode 100644 src/ml/__init__.py create mode 100644 src/ml/data.py create mode 100644 src/ml/features.py create mode 100644 src/ml/predict.py create mode 100644 src/ml/train.py create mode 100644 tests/test_flood_forecast.py diff --git a/.gitignore b/.gitignore index 22b9916..d66c2f5 100644 --- a/.gitignore +++ b/.gitignore @@ -137,3 +137,14 @@ vm_data/ grafana_data/ # Runtime station config (persisted CRUD); bundled default lives in src/data/ /stations.json + +# Ruflo local secrets and runtime data +.env.*.local +.claude-flow/data/ +.claude-flow/logs/ +.claude-flow/sessions/ + +# Trained flood-forecast model artifacts (produced on the server, ~100 MB; see docs/FLOOD_FORECASTING.md) +models/*.joblib +models/cache/ +models/metrics.json diff --git a/docs/FLOOD_FORECASTING.md b/docs/FLOOD_FORECASTING.md new file mode 100644 index 0000000..d0838f5 --- /dev/null +++ b/docs/FLOOD_FORECASTING.md @@ -0,0 +1,621 @@ +# Flood forecasting + +Short-range flood-risk forecasts for the Ping River gauge network, trained on the +monitor's own PostgreSQL history. This document covers what the system predicts, +what it is built from, how well it actually performs, how to run it on the server, +and when to retrain. + +Code lives in `src/ml/` (`data.py`, `features.py`, `train.py`, `predict.py`), the +training entry point is `scripts/train_flood_model.py`, tests are in +`tests/test_flood_forecast.py`, and trained artifacts land in `models/`. + +## 1. Overview + +For every station the system answers three questions at three lead times (6, 12 +and 24 hours): + +- **`p_warning`** — probability the water level reaches or exceeds the warning + threshold (3.0 m) at any point within the horizon. +- **`p_danger`** — same for the danger threshold (4.5 m). +- **`predicted_max_level`** — the expected peak level within the horizon, in + metres on the station's own datum. + +The window is open-ended forward: `exceed_warn_6` at 09:00 asks whether the level +touches 3.0 m anywhere in (09:00, 15:00], not what it will be at 15:00 — the +question an operator actually has. + +Fifteen of the sixteen stations have trained models. P.4A (Ban Mae Taeng) is +excluded by `features.NOT_TRAINABLE` (17.2% hourly fill, effectively dead +2019–2024: 289 rows in 2019, 769 in 2024) and is served by the persistence +heuristic instead. It still feeds downstream stations as an *input*, where +HistGradientBoosting's native NaN handling copes with the gaps. + +Every forecast row carries `source` (`model` or `heuristic`), `model_version` +and `trained_at`, so a stale or degraded forecast is visible in the payload +rather than silently indistinguishable from a good one. + +## 2. Data + +**Source of record.** PostgreSQL table `water_measurements` joined to `stations`. +As verified on 2026-08-10 (`inventory.json`, `db_cross_check`): **592,240 rows, +16 stations, 2018-08-01 through 2026-08-10**, zero mismatches against the HTTP +API, and a `status` column that is uniformly `active` (there are no quality flags +to filter on — bad readings must be caught by the feature pipeline, not by the +database). + +**Coverage is the dominant data constraint.** Readings are nominally hourly +(modal interval 1 h, ~95% of gaps), but only about **56% of hours on the complete +hourly grid have a reading**: P.1, P.67, P.76 and P.84 at 56.1%, P.103 at 55.8%, +P.21 at 55.4%, P.20 at 53.7%, P.87 at 53.3%, P.5 at 50.7%, P.4A at 17.2%. + +The missingness is **systematic, not random**. Measured over the 587 days from +2025-01-01 in `models/cache/P.1.csv.gz`, the fraction of days with a reading at +each hour is roughly 0.80 for 01:00–12:00, 0.68–0.69 for 13:00–16:00, 0.46–0.49 +for 17:00–21:00, 0.40–0.42 for 22:00–23:00, and 0.32 at midnight. That is a +scrape-schedule fingerprint, not hydrology — and it is why hour-of-day is +deliberately *not* a feature (see section 3). + +Long outages would poison training if used naively: P.87 lost 3,961 hours +(165 days) in 2023, P.20 lost 2,681 hours in 2021, P.77 lost 2,522 hours in early +2022, P.5 lost 2,429 hours over the 2020–21 turn. `features.TRAIN_START` excludes +P.5 before 2022-01-01; the rest are handled by the per-row coverage gates. + +**How `src/ml/data.py` loads it.** `resolve_db_url()` picks a connection string in +priority order: an explicit `--db-url` argument, then the `FLOOD_ML_DB_URL` +environment variable, then `Config.get_database_config()` when `DB_TYPE` is +`postgresql`, else `None`. `load_measurements()` then tries three tiers: + +1. **PostgreSQL** (`_fetch_from_db`) — the primary path. NULL discharge stays + NULL, which matters because the models must learn from the real missingness + pattern. +2. **HTTP API** (`_fetch_from_api`, default `http://100.81.167.42:8000`) — a + fallback for running off-server. **Caveat:** the public history endpoint + backfills missing discharge with a synthetic rating-curve estimate, so this + path is not equivalent to the DB path. It is flagged as + `discharge_maybe_synthetic: true` in the cache metadata. +3. **On-disk cache** (`models/cache/{station}.csv.gz` plus `meta.json`) — last + resort only. A successful DB or API fetch refreshes the cache; the cache is + never treated as a source of fresh data. + +`load_latest()` (used by the API) pulls the trailing 336 hours and never writes +the cache. + +## 3. Physics and features + +### Upstream routing + +The Ping mainstem gives real forecast skill for free: a flood wave takes hours to +travel downstream, so an upstream gauge reading *now* is information about a +downstream gauge *later*. `data-scout` measured these travel times by +cross-correlating water-level anomalies against the basin anchor P.1. The peak +correlation lags, which are hard-coded in `features.UPSTREAM_LEADS`: + +| Station | Lead vs P.1 | Peak anomaly correlation | Distance to P.1 (km) | +|---|---|---|---| +| P.20 (Ban Chiang Dao) | 17 h | 0.59 | 84.5 | +| P.92 (Ban Muang Aut) | 15 h | 0.66 | 63.8 | +| P.75 (Ban Chai Lat) | 12 h | 0.61 | 45.0 | +| P.4A (Ban Mae Taeng) | 12 h | 0.73 | 37.9 | +| P.67 (Ban Tae) | 7 h | 0.74 | 25.3 | +| P.21 (Ban Rim Tai) | 9 h | 0.56 | 15.0 | +| P.103 (Ring Bridge 3) | 1 h | 0.88 | 9.3 | + +(Distances are cumulative straight-line gauge-to-gauge, from the inventory's +`spatial_order_north_to_south`, not channel length — the real river is longer.) + +The lags are broadly consistent with distance, with one exception worth knowing +about: P.21 is 10 km closer to P.1 than P.67 yet lags by 9 h rather than 7 h, and +it has the weakest correlation of the mainstem set (0.56). Whatever the cause, +the table encodes the measured lag rather than the one distance would predict — +which is the point of measuring instead of assuming. + +Two stations are *downstream* of P.1 (P.5 at −12 h, P.81 at −4 h). For those, +`UPSTREAM_LEADS` routes P.1 and P.103 forward as their inputs, which is the same +physics running the other direction. + +Six western-tributary stations — P.82, P.84, P.87, P.77, P.85, P.76 — have empty +`UPSTREAM_LEADS` and are **un-routed**. Their anomaly correlations with P.1 are +0.22–0.32, low enough that routing them would inject noise rather than signal. +They are forecast from their own history plus the P.1 basin-state features. This +is a known gap: those catchments have no upstream gauge of their own in this +network. + +### Feature set (`features.build_features`) + +Everything is computed on an hourly grid built by `make_hourly_grid()`, which +keeps three aligned frames: `observed` (raw, NaN where nothing was recorded), +`filled` (forward-filled with `FFILL_LIMIT_H = 3`), and `mask` (True where a real +reading exists). Features read `filled`; labels read `observed` only. + +Per target station: + +- **Self level**: current level, lags at 1/2/3/6/12/24/48/72 h. +- **Rate of rise**: level minus its own value 1/3/6/12/24 h ago — a river at + 2.5 m and falling is a different situation from one at 2.5 m rising 30 cm/h. +- **Rolling statistics**: 6/24 h means, 6/24/72 h maxima, 24 h minimum. +- **Discharge**: current, lags at 6/24 h, 6 h rise (read from `observed`, so NULL + discharge stays NULL). +- **Observation health**: `obs_age_h` (hours since the last real reading, capped + at `FFILL_LIMIT_H`) and `cov_24h` (fraction of the last 24 h actually observed), + so the model can learn to hedge when a gauge is going quiet. +- **Routed upstream**, per `(upstream, lead)` pair: the upstream level at + `lead−3`, `lead`, and `lead+3` hours ago, its 6 h rise at `lead`, and its 24 h + rolling max at `lead−3`. The three-point bracket absorbs error in the measured + travel time rather than depending on it being exact. +- **Basin state** (non-P.1 stations only): P.1 level, its 24 h rolling max, and + its 24 h rise. +- **Seasonality**: `doy_sin`, `doy_cos` and an `is_monsoon` flag for June–October. + +**Hour-of-day is deliberately excluded.** Given the availability profile in +section 2, an hour-of-day feature would let the model learn "readings at 03:00 +are more likely to exist" and route that through to the label — an artefact of +when the scraper runs, with no hydrological content, that would evaporate the +moment the scrape schedule changed. + +### No-leakage guarantees + +- Only `shift()`, backward `rolling()` and forward-fill are used — nothing + interpolates, and no row can read a value timestamped after itself. +- `test_no_future_leakage` enforces this empirically: it adds +50 m to every + reading after time *t*, rebuilds the features, and asserts the rows at or + before *t* are bit-identical. +- Labels come from `observed`, never `filled`, so a forward-filled value can + never become its own target. +- The split is strictly temporal, and `early_stopping` is disabled in + `HGB_PARAMS` specifically because scikit-learn's internal validation split is + random and would leak across time. + +### Coverage gating + +A label is only trusted if enough of its forward window was actually observed. +`build_labels` requires `MIN_WINDOW_COVERAGE = 0.5` — at least half the horizon's +hours present — otherwise the label is NaN and the row is dropped from that head's +training set. The one exception is deliberate: **an observed exceedance always +produces a positive label regardless of coverage**, because a confirmed 3.5 m +reading inside a sparse window is not ambiguous. Rows whose own features are +stale (`obs_age_h` is NaN, i.e. the last real reading is more than 3 h old) are +dropped entirely in `build_matrix`. + +## 4. Models + +### Architecture + +One `HistGradientBoosting` model per **station × horizon × head**: + +| Head | Type | Target | +|---|---|---| +| `max_{h}` | `HistGradientBoostingRegressor` (squared error) | max observed level in (t, t+h] | +| `warn_{h}` | `HistGradientBoostingClassifier` | level ≥ 3.0 m anywhere in (t, t+h] | +| `danger_{h}` | `HistGradientBoostingClassifier` | level ≥ 4.5 m anywhere in (t, t+h] | + +Nine heads per station, three horizons (6/12/24 h), fifteen trained stations. +Hyperparameters are fixed (`HGB_PARAMS`: 300 iterations, learning rate 0.06, 31 +leaf nodes, minimum 50 samples per leaf, L2 1.0, `random_state=42`), chosen in an +earlier sweep and not re-searched per run — training is deterministic and +repeatable. + +HistGradientBoosting was chosen for three concrete reasons: it handles NaN +natively (essential given ~44% missing hours), it needs no feature scaling, and +it trains on CPU alone — no GPU anywhere in this pipeline (measured cost in +section 6). + +### Head gating and fallbacks + +The system degrades in tiers rather than failing: + +1. **Classifier head**, when the training span contains at least + `MIN_POSITIVES_FOR_CLASSIFIER = 30` positive examples. Below that, a + classifier would be fitting noise, and the head is recorded in + `skipped_heads` with its reason. +2. **Sigmoid on the regression head**, when the classifier is absent. + `p = 1/(1 + exp(−(predicted_max − threshold)/σ))`, where σ is the standard + deviation of the regressor's test residuals (floor `MIN_SIGMA = 0.15` m). This + turns the peak-level prediction into a calibrated-ish probability that widens + correctly when the regressor is less accurate at that horizon — at P.1, σ is + 0.15 m at 6 and 12 h but 0.166 m at 24 h. +3. **Persistence heuristic** (`predict._heuristic_forecast`), when there is no + model file at all, or the station's newest reading is more than + `STALE_AFTER_H = 6` hours old. It extrapolates the last 3 h rate of rise + forward with a 0.7 damping factor and a fixed σ of 0.3 m. It is not skilful; it + exists so the endpoint always returns something structurally valid. + +A station is skipped entirely if it has fewer than `MIN_ROWS_TO_TRAIN = 200` +usable rows; an individual head is skipped below `MIN_ROWS_FOR_HEAD = 50` labeled +rows. `_safe_fit` converts any fit failure (typically HistGradientBoosting's +binning step rejecting an all-NaN or constant column) into a recorded skip rather +than a station-killing exception. + +### Training procedure + +`train_station` runs two passes. First it evaluates on the strict temporal +holdout (train ≤ 2024-12-31, test 2025-01-01 → 2026-08-10) to produce the metrics +and the σ calibration. Then it **refits every head on the entire record** for the +deployed artifact, so the shipped model has seen the most recent data. Because +the full record has more labeled rows than the training half, the head-gating +decisions can differ between the two passes — `skipped_heads` is therefore +re-derived during the refit so it always describes what is actually in the saved +bundle, not what the evaluation pass decided. + +### Bundle format + +`models/flood_{station}.joblib` contains: `station_code`, `model_version` +(`hgb-v1+`), `trained_at`, `sklearn_version`, `feature_names`, +`horizons`, `thresholds`, `heads`, `sigma`, `skipped_heads`, `train_span`, and +`n_train_rows`. + +`feature_names` is the important one. At prediction time `_model_forecast` +rebuilds the feature row from live data and checks it against the bundle's stored +list; if any expected column is missing it logs an error and falls back to the +heuristic rather than feeding scikit-learn silently misaligned columns. +`test_feature_name_stability` guards the same invariant at build time. Bundles are +cached in memory keyed by `(path, mtime)`, so dropping in a retrained file +invalidates the cache without a restart. + +## 5. Measured performance + +### Holdout metrics (`models/metrics.json`) + +Model version `hgb-v1+49a3de0`, generated 2026-08-10. Train ≤ 2024-12-31, test +2025-01-01 → 2026-08-10 — the test span is entirely unseen future data relative +to training. + +P.1 (Nawarat Bridge), the station that matters most: + +| Horizon | Warning PR-AUC | Recall @1% FAR | Recall @5% FAR | MAE | MAE above 2 m | Test rows | Base rate | +|---|---|---|---|---|---|---|---| +| 6 h | 0.974 | 98.3% | 100% | 6.1 cm | 9.2 cm | 8,536 | 1.36% | +| 12 h | 0.904 | 93.8% | 97.7% | 9.0 cm | 15.0 cm | 7,932 | 1.61% | +| 24 h | 0.900 | 90.1% | 93.4% | 11.3 cm | 24.5 cm | 8,572 | 1.77% | + +Read PR-AUC against the base rate — 0.974 versus a 1.36% positive rate is a wide +margin over chance. "Recall at 1% false-alarm rate" is the operationally honest +number: at a threshold that fires on 1% of quiet hours, the 6 h model still +catches 98.3% of warning exceedances. + +P.103 (Ring Bridge 3) is the only station with enough danger-level events to +evaluate a danger head on the 2025–26 span (base rate 5.7–7.4%): PR-AUC 0.979 / +0.953 / 0.892 and recall at 1% FAR of 97.9% / 89.9% / 79.5% at 6 / 12 / 24 h. + +Across the other stations the 6 h warning PR-AUC spans 0.996 (P.5) down to 0.302 +(P.82), and tracks almost exactly with how many exceedances that station saw. The +strong ones are the frequently-flooded gauges — P.5 0.996, P.81 0.992, P.77 0.968, +P.85 0.953, P.75 0.927 — and the weak ones are un-routed western tributaries with +almost no positives (P.84 0.570, P.82 0.302 on 0.22% of test hours). P.92 and P.20 +have no evaluable warning metric at all: neither crossed 3.0 m often enough in the +test span (P.92 not once, P.20 in 0.09% of hours) to score. + +### Headline validation: the October 2024 record flood + +The holdout above never sees a true extreme, because the 2024 flood is in the +training half. So the model was retrained on data **ending 2024-08-31** and asked +to forecast September–November 2024 cold, with no knowledge of the event that +followed. This is the closest thing to a real operational test available. + +- **25 September cold start.** P.1's first warning crossing of the episode was + alerted **24–26 hours ahead**. This is the genuinely impressive case: the river + was in normal state, and the alert came from upstream routing alone. +- **5 October record peak** (P.1 5.30 m, P.103 9.93 m — the highest levels in the + eight-year record). Alerted **48 hours ahead**. Read this one carefully: the + river was already in sustained flood by then, so "48 hours" is the + `_first_alert_at` lookback window (`lookback_h = 48`) saturating, not a + measurement of true lead time. The model was correctly alarmed throughout; + the metric simply cannot express how much earlier than 48 h that started. +- **P.103 danger head** over the same window: PR-AUC 0.98–0.99, recall at 1% FAR + 87–95%. It called the danger-level crossings, not just the warning ones. +- **8 November re-flood.** Caught **26–31 hours ahead** by the 12 and 24 h models + — a second, independent event in the same test window. +- **P.103's 1 September "miss"** is a test-boundary artefact: the event begins in + the first hours of the test span, before the feature window has enough test-side + history to have produced a sustained alert. It is not a model failure, but it is + also not evidence of skill. + +### Honest limits + +**Genuine lead time is capped by gauge-only physics.** The longest upstream travel +time into P.1 is 17 h (P.20), and the strongest predictors are much closer: +P.103 at 1 h, P.67 at 7 h, P.21 at 9 h. Once a 24 h forecast reaches past roughly +17 h, there is no observation that has "already happened" to inform it — the model +is extrapolating basin state and season, not routing a wave. The 2025–26 test +events bear this out: the 25 September 2025 cold-start crossing was called 7 h +ahead by the 12 h model and 9 h ahead by the 24 h model. **Practical lead for P.1 +is ~7–17 h.** Extending it requires rainfall forecasts and Mae Ngat/Mae Kuang dam +release data, neither of which this system currently ingests. + +**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 +metric is `null`). The danger head exists and is trained on the full record — the +river has spent 57 hours above 4.5 m historically, 0.144% of all hours — but no +out-of-sample number backs it. Treat `p_danger` at P.1 as indicative, not +validated. + +**Thresholds are P.1-datum-specific but applied network-wide.** `THRESHOLDS` +carries a single default of (3.0, 4.5) m for every station, and gauge datums +differ enormously: P.5 spends 27.8% of all recorded hours above 3.0 m, P.103 +11.1%, P.81 8.2%, P.77 7.6% — versus 0.95% at P.1. On those stations "warning" is +close to a normal wet-season level, so **P.103 currently over-alerts on the +default thresholds**. The models are calibrated to the labels they were given, so +the metrics are internally valid; it is the operational meaning of the label that +is wrong. Per-station calibration against RID's published flood stages is the +pending follow-up — it changes only `features.THRESHOLDS`, plus a retrain. + +## 6. Deployment + +### API + +`GET /forecast` (`src/web_api.py`) returns one JSON row per station × horizon with +the fields listed in section 1. Results are cached in-process for +`FORECAST_TTL = 900` seconds (15 minutes), which matches the data cadence — the +underlying readings do not update faster than hourly. Inference runs in a thread +via `asyncio.to_thread` so it never blocks the event loop. + +Failure modes: **503** if the `src.ml` package cannot be imported (missing +scikit-learn, say), **502** on any other exception. + +One behaviour worth knowing, because the code comments suggest otherwise: the +endpoint's `FileNotFoundError` ("No trained flood models found") and `RuntimeError` +handlers are unreachable — nothing in `src/ml/` raises either, and +`predict._forecast_station` checks `bundle_path.exists()` and falls back to the +heuristic instead. So **before the first training run `/forecast` returns 200 with +an all-heuristic payload**, not a 503, provided there is recent gauge data; you +get an empty `200 []` only when there is no recent data at all. Judge deployment +state by the `source` field, not the status code. + +### Dashboard + +The "Flood risk outlook" panel (`src/static/dashboard.html`, `loadForecasts()`) +loads non-blocking after the map renders and **stays hidden unless `/forecast` +returns a non-empty array** — a non-OK response, an empty array, or a thrown +fetch all just leave the panel hidden, and the rest of the dashboard is +unaffected. Per the note above, this means the panel appears with heuristic-only +content once data is flowing but before any model is trained; the per-chip +tooltip is what tells you so. Stations are sorted worst-risk first, each showing +three chips (6/12/24 h) coloured by risk band, with the tooltip carrying the exact +warning and danger percentages, the predicted peak level, and a "heuristic +fallback" note when the row did not come from a model. The panel is labelled +*experimental*. + +### Training on the server + +The server already has the PostgreSQL connection configured, so no host override +is needed: + +```bash +cd /path/to/Northern-Thailand-Ping-River-Monitor +python scripts/train_flood_model.py --stations all +``` + +`resolve_db_url()` picks up `Config.get_database_config()` automatically when +`DB_TYPE=postgresql`. The run writes fifteen `models/flood_{station}.joblib` +bundles plus `models/metrics.json`. + +### Artifacts and dependencies + +The fifteen bundles total **101.4 MB** — mean 6.76 MB, from 4.04 MB (P.20) to +9.35 MB (P.103 and P.87, with P.5 next at 8.99 MB) — plus `metrics.json` at +0.34 MB and a 2.1 MB `models/cache/`. **These are not in +git**, and they should stay that way — artifacts are produced on the server, not +shipped. `.gitignore` excludes `models/*.joblib`, `models/cache/` and +`models/metrics.json` for exactly this reason. + +Two pins matter and are already in `requirements.txt` / `pyproject.toml`: +`scikit-learn==1.9.0` and `numpy>=1.24,<2` (pandas 2.0.3 wheels are ABI +incompatible with numpy 2.x). Bundles record `sklearn_version`; unpickling a +bundle under a different scikit-learn version is not guaranteed to work, so +retrain after any scikit-learn upgrade rather than assuming the artifacts carry +over. + +### Measured resource use + +All figures below were measured on 2026-08-10 on a development workstation — +**24 physical / 32 logical cores at 2.20 GHz, 32 GiB RAM** (Python 3.11.9, +scikit-learn 1.9.0, joblib 1.5.3, numpy 1.26.4, pandas 2.0.3) — **not** on the +production server. They come from two independent benchmark runs on that same +machine, which is why a couple of figures below are quoted as narrow ranges. +Treat the CPU times as a floor and the memory figures as representative, since +RSS barely depends on core count. Training read the `models/cache/` csv.gz files +(592,240 rows load in 0.5 s); loading the same history from PostgreSQL was not +measured and will be slower. + +**Training** (`train_all`, all 15 stations, evaluation pass plus full refit): + +| Measurement | Value | +|---|---| +| Full 15-station run, unrestricted threads | **199 s (3.3 min)** | +| Peak RSS during the full run | **209 MB** | +| Single station, unrestricted (P.1 / P.103) | 16.5 s / 18.7 s | + +HistGradientBoosting threads through OpenMP, and it scales only modestly. Timing +P.1 alone under `OMP_NUM_THREADS`: + +| Threads | 1 | 2 | 4 | unrestricted (32) | +|---|---|---|---|---| +| P.1 train time | 36.5 s | 23.5 s | 14.7 s | 16.5 s | + +Two things follow. **Four threads is the sweet spot** — 32 threads was marginally +*slower* than 4, so oversubscription costs you a little. And **even one core is +enough**: at 36.5 s per station, a single-core box retrains all fifteen in roughly +9 minutes (extrapolated, not measured end-to-end). + +Per station the fit costs **8–17 s**, and P.1 is the worst case at 16.9 s — it is +the basin anchor, so it carries 64 features against 32 for stations with fewer +upstream inputs (P.85 9.6 s, P.20 8.2 s). Two things are *not* the cost driver. +Evaluation isn't: P.1 with `skip_eval=True` took 17.3 s, no faster than the full +path. Nor is feature engineering — `build_matrix` over P.1's whole 8-year history +is 256 ms against 8–17 s of fitting. **The fit is the cost.** + +One honest caveat about the run that produced the current artifacts. By file +mtime it wrote all fifteen models between 11:55:29 and 12:01:15 — **5 min 46 s**, +averaging 25 s/station including joblib serialization, which lines up with the +measured fits. But `models/cache/meta.json` records the data fetch finishing at +11:45:49, so end to end that run spanned about 15.5 minutes, and the 9 min 40 s +gap between fetch and first model could not be reconstructed from the surviving +artifacts. Do not attribute it to per-station training cost. Either way the +conclusion holds: **retraining is minutes, not tens of minutes.** + +**Inference** (15 bundles, 16 stations × 3 horizons = 48 rows): + +| Measurement | Value | +|---|---| +| Cold call — every bundle unpickled from disk | **6.7 s** | +| Warm call — bundles in `_MODEL_CACHE` | **0.72 s** median (0.63–0.84 s) | +| RSS after imports, before any model | 71 MB | +| RSS with all 15 bundles resident | **288 MB** | + +The 101.4 MB of on-disk pickles expand to roughly **203–211 MB resident** — about +2× — and they stay there: `_MODEL_CACHE` replaces an entry when the file's mtime +changes but never drops one to reclaim memory. That is the single largest memory +cost of the whole feature. + +Where the time goes: cold start is 5.30 s, of which 0.46 s is the import and +4.84 s is unpickling, and 2.6 s of *that* is the first bundle alone paying a +one-time lazy `sklearn.ensemble` import — the remaining fourteen average 159 ms. +Of the ~640 ms warm compute, model prediction is ~525 ms, the hourly grid 47 ms, +and feature building 71 ms across all sixteen stations. + +**Live-endpoint measurements** (a second, independent benchmark run against a real +uvicorn instance of the app, same day, same workstation, RSS summed over the +process tree): + +| Measurement | Value | +|---|---| +| `/forecast` cache hit (15-min TTL) | **2.4 ms** median | +| `/forecast` cache miss, default threads | 15.2 s (≈7 s of that was the HTTP data fallback; a local DB replaces it) | +| `/forecast` cache miss, `OMP_NUM_THREADS=1` | 10.9 s | +| API process RSS, idle → models resident | 76 MB → **335 MB** | + +The endpoint-level RSS (335 MB) is higher than the models-only figure above +because the live process also retains the pandas frames from the data pull and +the HTTP/JSON machinery — use 335 MB as the sizing number. + +One threading subtlety cuts the other way in serving: inference is ~135 +single-row predicts, and at one row OpenMP thread dispatch costs more than the +math — `OMP_NUM_THREADS=1` makes the warm compute 2.6× faster (642 ms → 252 ms). +Training shows the opposite (2.3× slower single-threaded), so set the variable +per process, never globally. + +**Server sizing, in plain terms:** this is a small workload and almost any server +runs it. **RAM is the binding constraint, not CPU.** Budget about **1 GB for the +API process** so the ~335 MB steady state has headroom on top of the rest of the +app; training peaks at only ~210–315 MB and can share the same box. No GPU +anywhere. Pin thread counts per process — `OMP_NUM_THREADS=1` in the serving +unit, `OMP_NUM_THREADS=4` for retraining so it cannot monopolise every core while +the API is serving. And since a cold call costs seconds against a 2.4 ms cache +hit, consider warming `/forecast` once at startup rather than letting a user +absorb it. + +## 7. Retraining policy + +**Why it matters here specifically.** This is not a generic "models go stale" +argument: + +- **Channel geometry changes after every major flood.** Scour, deposition and + bank failure shift the level-to-discharge relationship at a gauge, and RID + revises rating curves after big events. A model trained on the pre-2024 channel + is predicting levels for a cross-section that no longer exists. +- **Extreme events extend the label range.** The highest P.103 reading before the + 2024 season was 7.54 m (October 2022); the 2024 event pushed it to 8.27 m on + 26 September and 9.93 m on 5 October. Gradient boosting cannot extrapolate past + its training range — predictions saturate at the largest value it has seen — so + every new record is what makes the next one predictable. +- **Station outages change feature availability.** P.87's 165-day gap in 2023 and + P.4A's five dead years mean the set of populated features drifts over time. + Retraining lets head gating and NaN handling re-adapt to the current sensors. + +**Recommended schedule:** + +| When | Why | +|---|---| +| **Every year, May–June (pre-monsoon)** | The minimum. Ensures the model entering the flood season has seen last season in full. | +| **Monthly, July–November** | Cheap insurance during the season — a full retrain costs minutes, not hours (section 6), so `nice` it and forget it. | +| **After any major flood event** | Non-negotiable. Channel geometry and rating curves have changed, and the new extreme extends the trainable label range. | + +Staleness is auditable without guesswork: `model_version` embeds the git short SHA +of the code that trained the bundle (`hgb-v1+49a3de0`), and `trained_at` is a +timestamp in every bundle. Both are echoed in every `/forecast` row, so you can +tell from the API response alone which code produced a forecast and how old the +model is. + +## 8. Operations runbook + +All commands assume the project virtualenv is active (`.venv` locally). + +**Train (all stations, with evaluation):** + +```bash +python scripts/train_flood_model.py --stations all +``` + +**Train a subset, refit-only (skips the holdout evaluation — much faster, but +produces no metrics and leaves σ at the `MIN_SIGMA` floor):** + +```bash +python scripts/train_flood_model.py --stations P.1,P.103 --skip-eval +``` + +**Train from a workstation against the server's database:** + +```bash +export FLOOD_ML_DB_URL='postgresql://user:pass@host:5432/dbname' +python scripts/train_flood_model.py --stations all +``` + +Do not commit that URL anywhere. If the DB is unreachable the loader silently +falls back to the HTTP API, whose discharge values are partly synthetic — check +the log line `PostgreSQL fetch failed, falling back to HTTP API` before trusting a +run. + +**Verify before promoting.** Training writes `models/metrics.json` alongside the +bundles. Check it before treating a run as good: + +```bash +python -c "import json; m=json.load(open('models/metrics.json')); \ +print(m['model_version'], m['split']); \ +print({s: v['status'] for s, v in m['stations'].items()}); \ +print({h: (d.get('pr_auc_warn'), d.get('mae')) for h, d in m['stations']['P.1']['per_horizon'].items()})" +``` + +Expect fifteen `trained` and one `heuristic` (P.4A). A station that reports +`failed` names its reason in the same payload. If P.1's 6 h warning PR-AUC has +dropped materially below ~0.97 or its MAE has risen well above ~6 cm, investigate +before deploying — that usually means a data problem (a gauge that went quiet, or +a bad backfill) rather than a modelling one. + +**Run the tests** (synthetic data only, no database or network required): + +```bash +python -m pytest tests/test_flood_forecast.py -v +``` + +Seven tests covering leakage, label alignment, the coverage gate, forward-fill and +staleness, a train/predict round trip, the heuristic fallback, and feature-name +stability. The whole suite runs in about 8 seconds, so there is no excuse for +skipping it before a deploy. + +**Understanding graceful degradation.** Three things can make a forecast row +non-model-backed, and all of them are visible in the payload: + +- `source: "heuristic"`, `model_version: "heuristic-v1"` — either no bundle exists + for that station (P.4A always, every station before the first training run), or + the station's newest reading is more than 6 hours old. +- A single horizon coming back heuristic while others are model-backed — that + horizon's head is in the bundle's `skipped_heads`, almost always because the + station had fewer than 30 positive examples for that threshold. +- A whole station flipping to heuristic after a code change — the feature-name + check in `_model_forecast` caught a mismatch between the live feature builder + and the stored `feature_names`. The fix is to retrain; the log line names the + missing columns. + +A live example from the 2026-08-10 cache: of 48 forecast rows, 42 came from +models and 6 were heuristic — three for P.4A, which has no bundle by design, and +three for P.92, whose newest reading was 02:00 while the basin's newest was 09:00. +That 7 hours of staleness crossed `STALE_AFTER_H = 6`, so P.92 correctly dropped +to persistence. Both fallback triggers, working as intended, in one ordinary call. + +Inspect a bundle's skipped heads directly (`joblib.load` unpickles, so only ever +point it at a bundle this pipeline's own `train.py` wrote — never a file from +elsewhere): + +```bash +python -c "import joblib; b=joblib.load('models/flood_P.1.joblib'); \ +print(b['model_version'], b['trained_at'], b['n_train_rows']); print(b['skipped_heads'])" +``` diff --git a/models/.gitkeep b/models/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/pyproject.toml b/pyproject.toml index 729f1aa..cb8f714 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -42,6 +42,9 @@ dependencies = [ "requests==2.31.0", "schedule==1.2.0", "pandas==2.0.3", + "numpy>=1.24,<2", + # Flood forecasting (ML) + "scikit-learn==1.9.0", # Web API framework "fastapi==0.104.1", "uvicorn[standard]==0.24.0", diff --git a/requirements.txt b/requirements.txt index 527c3b4..a52459b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,6 +2,10 @@ requests==2.31.0 schedule==1.2.0 pandas==2.0.3 +numpy>=1.24,<2 # pandas 2.0.3 wheels are ABI-incompatible with numpy 2.x + +# Flood forecasting (ML) +scikit-learn==1.9.0 # Web API framework fastapi==0.104.1 diff --git a/scripts/train_flood_model.py b/scripts/train_flood_model.py new file mode 100644 index 0000000..9e9c212 --- /dev/null +++ b/scripts/train_flood_model.py @@ -0,0 +1,17 @@ +#!/usr/bin/env python3 +"""CLI entry point for training the Ping River flood forecast models. + +Usage: + python scripts/train_flood_model.py --stations all + python scripts/train_flood_model.py --stations P.1,P.103 --skip-eval +""" + +import os +import sys + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..")) + +from src.ml.train import main + +if __name__ == "__main__": + main() diff --git a/src/data/stations.json b/src/data/stations.json index c849b1c..655a56d 100644 --- a/src/data/stations.json +++ b/src/data/stations.json @@ -83,8 +83,8 @@ "code": "P.81", "thai_name": "บ้านโป่ง", "english_name": "Ban Pong", - "latitude": 13.805661820610888, - "longitude": 99.87174946122846, + "latitude": 18.693611, + "longitude": 99.081944, "geohash": null }, "12": { diff --git a/src/ml/__init__.py b/src/ml/__init__.py new file mode 100644 index 0000000..5211975 --- /dev/null +++ b/src/ml/__init__.py @@ -0,0 +1 @@ +"""Flood forecasting ML package: data loading, feature/label engineering, training, and prediction.""" diff --git a/src/ml/data.py b/src/ml/data.py new file mode 100644 index 0000000..6684e79 --- /dev/null +++ b/src/ml/data.py @@ -0,0 +1,223 @@ +"""Loaders for flood-model training/inference data. + +Primary path reads raw measurements straight from PostgreSQL (keeping NULL +discharge as NULL). HTTP fallback goes through the public API's history +endpoint, which backfills missing discharge with a synthetic rating-curve +estimate -- callers are told about that via the `discharge_maybe_synthetic` +cache metadata flag. +""" + +import datetime +import gzip +import json +import logging +import os +from pathlib import Path +from typing import Dict, List, Optional + +import pandas as pd +from sqlalchemy import create_engine, text + +from ..config import Config +from .features import UPSTREAM_LEADS + +logger = logging.getLogger(__name__) + +DEFAULT_API_URL = "http://100.81.167.42:8000" +CACHE_DIR = Path("models/cache") +_MEASUREMENT_COLUMNS = ["timestamp", "station_code", "water_level", "discharge"] + + +def resolve_db_url(db_url: Optional[str] = None) -> Optional[str]: + """Resolve a Postgres connection string: explicit param > FLOOD_ML_DB_URL env > + Config's postgresql connection string > None (caller should fall back to HTTP).""" + if db_url: + return db_url + + env_url = os.getenv("FLOOD_ML_DB_URL") + if env_url: + return env_url + + try: + db_config = Config.get_database_config() + except Exception as error: + logger.warning(f"Could not resolve database config: {error}") + return None + + if db_config.get("type") == "postgresql": + return db_config.get("connection_string") + return None + + +def _default_stations() -> List[str]: + return list(UPSTREAM_LEADS.keys()) + + +def _normalize_long(df: pd.DataFrame) -> pd.DataFrame: + if df.empty: + return pd.DataFrame(columns=_MEASUREMENT_COLUMNS) + df = df.copy() + df["timestamp"] = pd.to_datetime(df["timestamp"]).dt.floor("h") + df["water_level"] = pd.to_numeric(df["water_level"], errors="coerce") + df["discharge"] = pd.to_numeric(df["discharge"], errors="coerce") + df = df.drop_duplicates(subset=["station_code", "timestamp"], keep="last") + df = df.sort_values("timestamp").reset_index(drop=True) + return df[_MEASUREMENT_COLUMNS] + + +def _fetch_from_db( + db_url: str, + stations: Optional[List[str]], + start: Optional[datetime.datetime], + end: Optional[datetime.datetime], +) -> pd.DataFrame: + engine = create_engine(db_url, pool_pre_ping=True) + query = ( + "SELECT m.timestamp, s.station_code, m.water_level, m.discharge " + "FROM water_measurements m JOIN stations s ON m.station_id = s.id WHERE 1=1" + ) + params: Dict = {} + if start is not None: + query += " AND m.timestamp >= :start_time" + params["start_time"] = start + if end is not None: + query += " AND m.timestamp <= :end_time" + params["end_time"] = end + if stations: + placeholders = ", ".join(f":station_{i}" for i in range(len(stations))) + query += f" AND s.station_code IN ({placeholders})" + for i, code in enumerate(stations): + params[f"station_{i}"] = code + query += " ORDER BY m.timestamp" + + with engine.connect() as connection: + df = pd.read_sql(text(query), connection, params=params) + return _normalize_long(df) + + +def _fetch_station_from_api(api_url: str, station_code: str, hours: int, limit: int = 100000) -> pd.DataFrame: + import requests + + response = requests.get( + f"{api_url}/measurements/history/{station_code}", + params={"hours": hours, "limit": limit}, + timeout=30, + ) + response.raise_for_status() + rows = response.json() + for row in rows: + row["station_code"] = station_code + return pd.DataFrame(rows, columns=_MEASUREMENT_COLUMNS + ["discharge_percent"]) + + +def _fetch_from_api( + api_url: str, + stations: List[str], + start: Optional[datetime.datetime], + end: Optional[datetime.datetime], +) -> pd.DataFrame: + now = datetime.datetime.now() + reference_end = end or now + reference_start = start or (reference_end - datetime.timedelta(days=365 * 8)) + hours = max(1, int((reference_end - reference_start).total_seconds() // 3600) + 1) + + frames = [] + for code in stations: + try: + frames.append(_fetch_station_from_api(api_url, code, hours)) + except Exception as error: + logger.warning(f"HTTP fallback failed for station {code}: {error}") + if not frames: + return pd.DataFrame(columns=_MEASUREMENT_COLUMNS) + df = pd.concat(frames, ignore_index=True) + return _normalize_long(df) + + +def _write_cache(df: pd.DataFrame, cache_dir: Path, source: str, discharge_maybe_synthetic: bool) -> None: + cache_dir.mkdir(parents=True, exist_ok=True) + for code, group in df.groupby("station_code"): + path = cache_dir / f"{code}.csv.gz" + with gzip.open(path, "wt", encoding="utf-8", newline="") as handle: + group.to_csv(handle, index=False) + meta = { + "fetched_at": datetime.datetime.now().isoformat(), + "source": source, + "discharge_maybe_synthetic": discharge_maybe_synthetic, + } + with open(cache_dir / "meta.json", "w", encoding="utf-8") as handle: + json.dump(meta, handle) + + +def _read_cache(cache_dir: Path, stations: Optional[List[str]]) -> pd.DataFrame: + if not cache_dir.exists(): + return pd.DataFrame(columns=_MEASUREMENT_COLUMNS) + frames = [] + for path in sorted(cache_dir.glob("*.csv.gz")): + code = path.name[: -len(".csv.gz")] + if stations and code not in stations: + continue + with gzip.open(path, "rt", encoding="utf-8") as handle: + frames.append(pd.read_csv(handle, parse_dates=["timestamp"])) + if not frames: + return pd.DataFrame(columns=_MEASUREMENT_COLUMNS) + return _normalize_long(pd.concat(frames, ignore_index=True)) + + +def load_measurements( + db_url: Optional[str] = None, + stations: Optional[List[str]] = None, + start: Optional[datetime.datetime] = None, + end: Optional[datetime.datetime] = None, + use_cache: bool = True, + cache_dir: Path = CACHE_DIR, + api_url: str = DEFAULT_API_URL, +) -> pd.DataFrame: + """Load the long-format [timestamp, station_code, water_level, discharge] history. + + Tries PostgreSQL first, then the HTTP API, then the on-disk cache as a last + resort. A successful DB/API fetch refreshes the cache; the cache itself is + never treated as a source of fresh data. + """ + resolved_db_url = resolve_db_url(db_url) + + if resolved_db_url: + try: + df = _fetch_from_db(resolved_db_url, stations, start, end) + if use_cache: + _write_cache(df, cache_dir, source="postgres", discharge_maybe_synthetic=False) + return df + except Exception as error: + logger.warning(f"PostgreSQL fetch failed, falling back to HTTP API: {error}") + + try: + api_stations = stations or _default_stations() + df = _fetch_from_api(api_url, api_stations, start, end) + if not df.empty: + if use_cache: + _write_cache(df, cache_dir, source="api", discharge_maybe_synthetic=True) + return df + except Exception as error: + logger.warning(f"HTTP API fetch failed: {error}") + + if use_cache: + logger.warning("Falling back to on-disk cache for measurement history") + return _read_cache(cache_dir, stations) + + return pd.DataFrame(columns=_MEASUREMENT_COLUMNS) + + +def load_latest( + db_url: Optional[str] = None, + hours: int = 336, + stations: Optional[List[str]] = None, +) -> pd.DataFrame: + """Load the last `hours` of history for all (or given) stations. Never cached to disk.""" + end = datetime.datetime.now() + start = end - datetime.timedelta(hours=hours) + return load_measurements( + db_url=db_url, + stations=stations, + start=start, + end=end, + use_cache=False, + ) diff --git a/src/ml/features.py b/src/ml/features.py new file mode 100644 index 0000000..66ab020 --- /dev/null +++ b/src/ml/features.py @@ -0,0 +1,272 @@ +"""Static config and feature/label engineering for the Ping River flood forecast models. + +All feature computation is strictly causal (no row uses information timestamped after +itself) so it is safe to run identically at training time and at prediction time. +""" + +import datetime +import logging +from dataclasses import dataclass +from typing import Dict, List, Optional, Tuple + +import numpy as np +import pandas as pd + +logger = logging.getLogger(__name__) + +# --------------------------------------------------------------------------- +# Static configuration +# --------------------------------------------------------------------------- + +# Per-station (warning, danger) level thresholds in meters. "*" is the default +# applied to any station without an explicit override. +THRESHOLDS: Dict[str, Tuple[float, float]] = { + "*": (3.0, 4.5), +} + +MONSOON_MONTHS = {6, 7, 8, 9, 10} +FFILL_LIMIT_H = 3 +MIN_WINDOW_COVERAGE = 0.5 + +BASIN_ANCHOR = "P.1" + +# Empirical hours a station's water-level anomaly leads the basin anchor (P.1), +# derived from data-scout cross-correlation analysis. UPSTREAM_LEADS[station] +# lists, for each station, the (upstream_code, lead_hours) pairs to use as +# routed-upstream input features when forecasting `station`. +UPSTREAM_LEADS: Dict[str, List[Tuple[str, int]]] = { + "P.1": [("P.103", 1), ("P.67", 7), ("P.21", 9), ("P.75", 12), ("P.4A", 12), ("P.92", 15), ("P.20", 17)], + "P.103": [("P.67", 6), ("P.21", 8), ("P.75", 11), ("P.4A", 11), ("P.92", 14), ("P.20", 16)], + "P.21": [("P.67", 1), ("P.75", 3), ("P.4A", 3), ("P.92", 6), ("P.20", 8)], + "P.67": [("P.75", 5), ("P.4A", 5), ("P.92", 8), ("P.20", 10)], + "P.75": [("P.92", 3), ("P.20", 5)], + "P.4A": [("P.92", 3), ("P.20", 5)], + "P.92": [("P.20", 2)], + "P.20": [], + "P.5": [("P.1", 12), ("P.103", 13)], + "P.81": [("P.1", 4), ("P.103", 5)], + "P.82": [], + "P.84": [], + "P.87": [], + "P.77": [], + "P.85": [], + "P.76": [], +} + +# Per-station usable-from dates: data before this cutoff is excluded from training +# because of known data-quality holes (see data-scout inventory). +TRAIN_START: Dict[str, str] = {"P.5": "2022-01-01"} + +# Stations with data too sparse/broken to ever be a regression/classification +# target. They are still usable as upstream *input* features (HGB tolerates NaN). +NOT_TRAINABLE: Dict[str, str] = {"P.4A": "17% fill, dead 2019-2024"} + + +def get_thresholds(station_code: str) -> Tuple[float, float]: + """Return (warning, danger) level thresholds for a station, falling back to the default.""" + return THRESHOLDS.get(station_code, THRESHOLDS["*"]) + + +# --------------------------------------------------------------------------- +# Hourly grid +# --------------------------------------------------------------------------- + + +@dataclass +class HourlyGrid: + """A complete hourly time grid pivoted wide across stations. + + observed: raw values, NaN where nothing was recorded that hour (pristine; used for labels). + filled: observed forward-filled per column with limit=FFILL_LIMIT_H (causal; used for features). + mask: boolean, True where `observed` has a real reading. + """ + + observed: pd.DataFrame + filled: pd.DataFrame + mask: pd.DataFrame + + +def make_hourly_grid(df_long: pd.DataFrame) -> HourlyGrid: + """Pivot a long station/timestamp measurement frame onto a complete hourly grid. + + df_long columns: timestamp, station_code, water_level, discharge. + """ + if df_long.empty: + empty = pd.DataFrame( + index=pd.DatetimeIndex([], name="timestamp"), + columns=pd.MultiIndex.from_tuples([], names=["station_code", "field"]), + ) + return HourlyGrid(observed=empty, filled=empty.copy(), mask=empty.copy()) + + df = df_long.copy() + df["timestamp"] = pd.to_datetime(df["timestamp"]).dt.floor("h") + df = df.drop_duplicates(subset=["station_code", "timestamp"], keep="last") + + full_index = pd.date_range(df["timestamp"].min(), df["timestamp"].max(), freq="h", name="timestamp") + + wide = df.pivot(index="timestamp", columns="station_code", values=["water_level", "discharge"]) + wide = wide.reorder_levels([1, 0], axis=1).sort_index(axis=1) + wide = wide.reindex(full_index) + + observed = wide + mask = observed.notna() + # Forward-fill only — never interpolate — so no row ever depends on a future value. + filled = observed.ffill(limit=FFILL_LIMIT_H) + + return HourlyGrid(observed=observed, filled=filled, mask=mask) + + +def _series(grid_frame: pd.DataFrame, station: str, field: str, index: pd.Index) -> pd.Series: + """Fetch a (station, field) column, or an all-NaN series if the station is absent.""" + if (station, field) in grid_frame.columns: + return grid_frame[(station, field)] + return pd.Series(np.nan, index=index) + + +def _hours_since_observed(mask_col: pd.Series) -> pd.Series: + """Hours since the last True in `mask_col` (0 at an observed hour; NaN if never observed yet).""" + idx = mask_col.index + obs_time = pd.Series(idx, index=idx).where(mask_col.to_numpy()) + last_obs_time = obs_time.ffill() + age_hours = (idx.to_series() - last_obs_time).dt.total_seconds() / 3600.0 + return age_hours + + +# --------------------------------------------------------------------------- +# Features +# --------------------------------------------------------------------------- + + +def build_features(grid: HourlyGrid, station: str) -> pd.DataFrame: + """Build the deterministic-order feature matrix for one target station.""" + idx = grid.observed.index + cols: Dict[str, pd.Series] = {} + + level = _series(grid.filled, station, "water_level", idx) + discharge = _series(grid.observed, station, "discharge", idx) + obs_mask = _series(grid.mask, station, "water_level", idx).fillna(False) + + cols["level"] = level + for k in (1, 2, 3, 6, 12, 24, 48, 72): + cols[f"level_lag_{k}"] = level.shift(k) + for k in (1, 3, 6, 12, 24): + cols[f"rise_{k}"] = level - level.shift(k) + cols["roll_mean_6"] = level.rolling(6, min_periods=1).mean() + cols["roll_mean_24"] = level.rolling(24, min_periods=1).mean() + cols["roll_max_6"] = level.rolling(6, min_periods=1).max() + cols["roll_max_24"] = level.rolling(24, min_periods=1).max() + cols["roll_max_72"] = level.rolling(72, min_periods=1).max() + cols["roll_min_24"] = level.rolling(24, min_periods=1).min() + + cols["discharge"] = discharge + cols["discharge_lag_6"] = discharge.shift(6) + cols["discharge_lag_24"] = discharge.shift(24) + cols["discharge_rise_6"] = discharge - discharge.shift(6) + + obs_age_h = _hours_since_observed(obs_mask) + cols["obs_age_h"] = obs_age_h.where(obs_age_h <= FFILL_LIMIT_H) + cols["cov_24h"] = obs_mask.rolling(24, min_periods=1).mean() + + for upstream_code, lead_h in UPSTREAM_LEADS.get(station, []): + u_level = _series(grid.filled, upstream_code, "water_level", idx) + u_rise_6 = u_level - u_level.shift(6) + u_rollmax_24 = u_level.rolling(24, min_periods=1).max() + near_lag = max(0, lead_h - 3) + cols[f"{upstream_code}_level_lag_{near_lag}"] = u_level.shift(near_lag) + cols[f"{upstream_code}_level_lag_{lead_h}"] = u_level.shift(lead_h) + cols[f"{upstream_code}_level_lag_{lead_h + 3}"] = u_level.shift(lead_h + 3) + cols[f"{upstream_code}_rise_6_lag_{lead_h}"] = u_rise_6.shift(lead_h) + cols[f"{upstream_code}_rollmax_24_lag_{near_lag}"] = u_rollmax_24.shift(near_lag) + + if station != BASIN_ANCHOR: + p1_level = _series(grid.filled, BASIN_ANCHOR, "water_level", idx) + cols["P1_level"] = p1_level + cols["P1_rollmax_24"] = p1_level.rolling(24, min_periods=1).max() + cols["P1_rise_24"] = p1_level - p1_level.shift(24) + + doy = idx.to_series().dt.dayofyear.astype(float) + cols["doy_sin"] = np.sin(2 * np.pi * doy / 365.25) + cols["doy_cos"] = np.cos(2 * np.pi * doy / 365.25) + cols["is_monsoon"] = idx.to_series().dt.month.isin(MONSOON_MONTHS).astype(float) + + return pd.DataFrame(cols, index=idx) + + +# --------------------------------------------------------------------------- +# Labels +# --------------------------------------------------------------------------- + + +def _future_window_stats(col: pd.Series, horizon_h: int) -> Tuple[pd.Series, pd.Series]: + """For every t, (max, count) of observed values in the OPEN window (t, t+horizon_h].""" + reversed_col = col.iloc[::-1] + shifted = reversed_col.shift(1) # excludes t itself + fut_max = shifted.rolling(horizon_h, min_periods=1).max().iloc[::-1] + fut_count = shifted.rolling(horizon_h, min_periods=1).count().iloc[::-1] + return fut_max, fut_count + + +def build_labels(grid: HourlyGrid, station: str, horizons: Tuple[int, ...] = (6, 12, 24)) -> pd.DataFrame: + """Build max-level and threshold-exceedance labels for one target station.""" + idx = grid.observed.index + observed_level = _series(grid.observed, station, "water_level", idx) + warn_thr, danger_thr = get_thresholds(station) + + out: Dict[str, pd.Series] = {} + for horizon_h in horizons: + fut_max, fut_count = _future_window_stats(observed_level, horizon_h) + cov = fut_count / horizon_h + enough_cov = cov >= MIN_WINDOW_COVERAGE + + exceed_warn = pd.Series(np.nan, index=idx) + exceed_warn[fut_max >= warn_thr] = 1.0 + exceed_warn[enough_cov & exceed_warn.isna()] = 0.0 + + exceed_danger = pd.Series(np.nan, index=idx) + exceed_danger[fut_max >= danger_thr] = 1.0 + exceed_danger[enough_cov & exceed_danger.isna()] = 0.0 + + max_level_valid = fut_max.where(enough_cov | (fut_max >= warn_thr)) + + out[f"max_level_{horizon_h}"] = max_level_valid + out[f"exceed_warn_{horizon_h}"] = exceed_warn + out[f"exceed_danger_{horizon_h}"] = exceed_danger + + return pd.DataFrame(out, index=idx) + + +# --------------------------------------------------------------------------- +# Glue +# --------------------------------------------------------------------------- + + +def build_matrix( + df_long: pd.DataFrame, + station: str, + horizons: Tuple[int, ...] = (6, 12, 24), +) -> 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) + + keep = X["obs_age_h"].notna() + train_start = TRAIN_START.get(station) + if train_start: + keep &= X.index >= pd.Timestamp(train_start) + + X = X.loc[keep] + Y = Y.loc[keep] + + positive_counts = { + col: int(Y[col].sum()) for col in Y.columns if col.startswith("exceed_") and Y[col].notna().any() + } + meta = { + "station_code": station, + "n_rows": int(len(X)), + "span": ( + (X.index.min().isoformat(), X.index.max().isoformat()) if len(X) else (None, None) + ), + "positive_counts": positive_counts, + } + return X, Y, meta diff --git a/src/ml/predict.py b/src/ml/predict.py new file mode 100644 index 0000000..f4f1cd6 --- /dev/null +++ b/src/ml/predict.py @@ -0,0 +1,284 @@ +"""Flood forecast inference. + +Integration contract (see get_forecasts / get_latest_forecasts): callers pass +raw station readings, get back one forecast dict per station x horizon. A +station with a stale, missing, or version-mismatched model transparently +falls back to a simple persistence heuristic instead of raising -- this +module must never crash the caller (e.g. the web API). +""" + +import datetime +import logging +from pathlib import Path +from typing import Dict, List, Optional, Tuple, Union + +import joblib +import numpy as np +import pandas as pd + +from . import features + +logger = logging.getLogger(__name__) + +DEFAULT_HORIZONS: Tuple[int, ...] = (6, 12, 24) +STALE_AFTER_H = 6.0 +HEURISTIC_SIGMA = 0.3 +HEURISTIC_VERSION = "heuristic-v1" + +# Keyed by (path, mtime) so a retrained model (new mtime) invalidates the old entry. +_MODEL_CACHE: Dict[Tuple[str, float], dict] = {} + + +def _load_bundle(path: Path) -> dict: + # joblib.load runs arbitrary pickle code; safe here because `path` is always + # models/flood_{station}.joblib, an artifact this pipeline's own train.py wrote -- + # never a user- or network-supplied file. + key = (str(path), path.stat().st_mtime) + cached = _MODEL_CACHE.get(key) + if cached is not None: + return cached + bundle = joblib.load(path) + for stale_key in [k for k in _MODEL_CACHE if k[0] == str(path)]: + del _MODEL_CACHE[stale_key] + _MODEL_CACHE[key] = bundle + return bundle + + +def _readings_to_long_df(readings_by_station: Dict[str, List[dict]]) -> pd.DataFrame: + rows = [] + for station_code, readings in readings_by_station.items(): + for reading in readings: + timestamp = reading.get("timestamp") + if isinstance(timestamp, str): + timestamp = pd.to_datetime(timestamp) + rows.append( + { + "timestamp": timestamp, + "station_code": station_code, + "water_level": reading.get("water_level"), + "discharge": reading.get("discharge"), + } + ) + if not rows: + return pd.DataFrame(columns=["timestamp", "station_code", "water_level", "discharge"]) + df = pd.DataFrame(rows) + return df.dropna(subset=["timestamp"]) + + +def _clip_probability(value: float) -> float: + return float(min(max(value, 0.0), 1.0)) + + +def _sigmoid_probability(predicted_max: float, threshold: float, sigma: float) -> float: + return 1.0 / (1.0 + np.exp(-(predicted_max - threshold) / sigma)) + + +def _heuristic_forecast( + station_code: str, + as_of: pd.Timestamp, + current_level: float, + level_t_minus_3: Optional[float], + warn_thr: float, + danger_thr: float, + horizons: Tuple[int, ...], +) -> List[dict]: + if level_t_minus_3 is None: + rate = 0.0 + else: + rate = max(0.0, (current_level - level_t_minus_3) / 3.0) + + results = [] + for horizon_h in horizons: + predicted_max = max(current_level + rate * horizon_h * 0.7, current_level) + p_warning = _clip_probability(_sigmoid_probability(predicted_max, warn_thr, HEURISTIC_SIGMA)) + p_danger = _clip_probability(_sigmoid_probability(predicted_max, danger_thr, HEURISTIC_SIGMA)) + p_danger = min(p_danger, p_warning) + results.append( + { + "station_code": station_code, + "horizon_hours": horizon_h, + "p_warning": p_warning, + "p_danger": p_danger, + "predicted_max_level": predicted_max, + "current_level": current_level, + "as_of": as_of.isoformat(), + "model_version": HEURISTIC_VERSION, + "trained_at": None, + "source": "heuristic", + "threshold_warning": warn_thr, + "threshold_danger": danger_thr, + } + ) + return results + + +def _model_forecast( + station_code: str, + grid: features.HourlyGrid, + bundle: dict, + as_of: pd.Timestamp, + current_level: float, +) -> List[dict]: + warn_thr = bundle["thresholds"]["warning"] + danger_thr = bundle["thresholds"]["danger"] + + feature_row = features.build_features(grid, station_code).loc[[as_of]] + expected_columns = bundle["feature_names"] + missing = [c for c in expected_columns if c not in feature_row.columns] + if missing: + logger.error(f"Feature mismatch for {station_code} (missing {missing}); falling back to heuristic") + return None + feature_row = feature_row[expected_columns] + + results = [] + for horizon_h in bundle["horizons"]: + reg = bundle["heads"].get(f"max_{horizon_h}") + if reg is None: + results.append(None) + continue + predicted_max = max(float(reg.predict(feature_row)[0]), current_level) + sigma_h = bundle["sigma"].get(horizon_h, HEURISTIC_SIGMA) + + warn_head = bundle["heads"].get(f"warn_{horizon_h}") + if warn_head is not None: + p_warning = float(warn_head.predict_proba(feature_row)[0][1]) + else: + p_warning = _sigmoid_probability(predicted_max, warn_thr, sigma_h) + + danger_head = bundle["heads"].get(f"danger_{horizon_h}") + if danger_head is not None: + p_danger = float(danger_head.predict_proba(feature_row)[0][1]) + else: + p_danger = _sigmoid_probability(predicted_max, danger_thr, sigma_h) + + p_warning = _clip_probability(p_warning) + p_danger = min(_clip_probability(p_danger), p_warning) + + results.append( + { + "station_code": station_code, + "horizon_hours": horizon_h, + "p_warning": p_warning, + "p_danger": p_danger, + "predicted_max_level": predicted_max, + "current_level": current_level, + "as_of": as_of.isoformat(), + "model_version": bundle["model_version"], + "trained_at": bundle["trained_at"], + "source": "model", + "threshold_warning": warn_thr, + "threshold_danger": danger_thr, + } + ) + return results + + +def _forecast_station( + station_code: str, + grid: features.HourlyGrid, + models_dir: Path, + now: pd.Timestamp, + horizons: Tuple[int, ...], +) -> List[dict]: + level_col = (station_code, "water_level") + if level_col not in grid.observed.columns: + logger.warning(f"No data for station {station_code}; omitting") + return [] + observed_level = grid.observed[level_col].dropna() + if observed_level.empty: + logger.warning(f"No observed readings for station {station_code}; omitting") + return [] + + as_of = observed_level.index.max() + current_level = float(observed_level.loc[as_of]) + staleness_h = (pd.Timestamp(now) - as_of).total_seconds() / 3600.0 + + warn_thr, danger_thr = features.get_thresholds(station_code) + t_minus_3 = as_of - pd.Timedelta(hours=3) + level_t_minus_3 = float(observed_level.loc[t_minus_3]) if t_minus_3 in observed_level.index else None + + bundle_path = models_dir / f"flood_{station_code}.joblib" + if not bundle_path.exists() or staleness_h > STALE_AFTER_H: + return _heuristic_forecast( + station_code, as_of, current_level, level_t_minus_3, warn_thr, danger_thr, horizons + ) + + bundle = _load_bundle(bundle_path) + model_results = _model_forecast(station_code, grid, bundle, as_of, current_level) + if model_results is None: + return _heuristic_forecast( + station_code, as_of, current_level, level_t_minus_3, warn_thr, danger_thr, horizons + ) + + # Per-horizon heads that were skipped at train time (e.g. too few positives) still + # need a forecast row -- fall back to the single-horizon heuristic for just that row. + filled = [] + for horizon_h, row in zip(bundle["horizons"], model_results): + if row is not None: + filled.append(row) + else: + filled.extend( + _heuristic_forecast( + station_code, as_of, current_level, level_t_minus_3, warn_thr, danger_thr, (horizon_h,) + ) + ) + return filled + + +def get_forecasts( + readings_by_station: Dict[str, List[dict]], + models_dir: Union[str, Path] = "models", + now: Optional[Union[datetime.datetime, str]] = None, +) -> List[dict]: + """Produce flood forecasts for every station present in `readings_by_station`. + + Each reading dict needs at least {timestamp, water_level, discharge}; extra + keys are ignored so raw API/DB rows can be passed straight through. At + least 96 hours of span is required to populate every feature; 336 hours + (14 days) is recommended. + """ + models_dir = Path(models_dir) + if now is None: + now = datetime.datetime.now() + now = pd.Timestamp(now) + + df_long = _readings_to_long_df(readings_by_station) + if df_long.empty: + return [] + + grid = features.make_hourly_grid(df_long) + results: List[dict] = [] + for station_code in readings_by_station.keys(): + try: + results.extend(_forecast_station(station_code, grid, models_dir, now, DEFAULT_HORIZONS)) + except Exception as error: + logger.error(f"Forecast failed for station {station_code}: {error}") + return results + + +def get_latest_forecasts( + db_url: Optional[str] = None, + models_dir: Union[str, Path] = "models", + hours: int = 336, +) -> List[dict]: + """Convenience wrapper for web_api: load the latest window from the DB/API and forecast. + + Raises FileNotFoundError when no trained model bundle exists at all, so the + API can 503 instead of serving purely heuristic output as if it were a forecast. + """ + from .data import load_latest + + if not sorted(Path(models_dir).glob("flood_*.joblib")): + raise FileNotFoundError(f"no trained model bundles in {models_dir}") + + df_long = load_latest(db_url=db_url, hours=hours) + readings_by_station: Dict[str, List[dict]] = {} + if not df_long.empty: + for station_code, group in df_long.groupby("station_code"): + readings_by_station[station_code] = group[["timestamp", "water_level", "discharge"]].to_dict("records") + + expected_stations = set(features.UPSTREAM_LEADS.keys()) + for missing_station in expected_stations - set(readings_by_station.keys()): + logger.warning(f"No recent data for station {missing_station}; omitting from forecasts") + + return get_forecasts(readings_by_station, models_dir=models_dir) diff --git a/src/ml/train.py b/src/ml/train.py new file mode 100644 index 0000000..49331c6 --- /dev/null +++ b/src/ml/train.py @@ -0,0 +1,414 @@ +"""Training CLI for the Ping River flood forecast models. + +Per station: build the feature/label matrix once, evaluate with a strict +temporal holdout (Split B), then refit each head on the full record for the +deployed artifact. Hyperparameters are fixed (chosen via an earlier Split A +sweep, not repeated here) -- no random search, no shuffling, no sklearn +early_stopping (its internal validation split is random and would leak +across time). +""" + +import argparse +import datetime +import json +import logging +import subprocess +from pathlib import Path +from typing import Dict, List, Optional, Tuple + +import joblib +import numpy as np +import pandas as pd +import sklearn +from sklearn.ensemble import HistGradientBoostingClassifier, HistGradientBoostingRegressor +from sklearn.metrics import average_precision_score, brier_score_loss, mean_absolute_error, mean_squared_error + +from . import features +from .data import DEFAULT_API_URL, load_measurements, resolve_db_url + +logger = logging.getLogger(__name__) + +HORIZONS: Tuple[int, ...] = (6, 12, 24) +SPLIT_B_TRAIN_END = "2024-12-31" +SPLIT_B_TEST_START = "2025-01-01" +SPLIT_B_TEST_END = "2026-08-10" +MIN_POSITIVES_FOR_CLASSIFIER = 30 +MIN_SIGMA = 0.15 +MIN_ROWS_TO_TRAIN = 200 +MIN_ROWS_FOR_HEAD = 50 + +HGB_PARAMS = { + "max_iter": 300, + "learning_rate": 0.06, + "max_leaf_nodes": 31, + "min_samples_leaf": 50, + "l2_regularization": 1.0, + "early_stopping": False, + "random_state": 42, +} + + +def _git_short_sha() -> str: + try: + result = subprocess.run( + ["git", "rev-parse", "--short", "HEAD"], capture_output=True, text=True, timeout=5, check=True + ) + sha = result.stdout.strip() + return sha or "nogit" + except Exception: + return "nogit" + + +def _make_regressor(overrides: Optional[dict] = None) -> HistGradientBoostingRegressor: + params = {**HGB_PARAMS, **(overrides or {})} + return HistGradientBoostingRegressor(loss="squared_error", **params) + + +def _make_classifier(overrides: Optional[dict] = None) -> HistGradientBoostingClassifier: + params = {**HGB_PARAMS, **(overrides or {})} + return HistGradientBoostingClassifier(**params) + + +def _safe_fit(estimator, X: pd.DataFrame, y: pd.Series, head_key: str, skipped_heads: Dict[str, str]): + """Fit an estimator, converting any failure (e.g. HistGradientBoosting's binning + step rejecting an all-NaN/constant feature column) into a recorded skip rather + than a station-killing exception.""" + try: + estimator.fit(X, y) + return estimator + except Exception as error: + skipped_heads[head_key] = f"fit failed: {error}" + logger.warning(f"{head_key}: fit failed, skipping ({error})") + return None + + +def _recall_at_far(y_true: np.ndarray, y_score: np.ndarray, target_far: float) -> Optional[float]: + """Recall at the score threshold whose false-positive rate over true negatives is <= target_far.""" + y_true = np.asarray(y_true) + y_score = np.asarray(y_score) + neg_scores = np.sort(y_score[y_true == 0])[::-1] + n_pos = int((y_true == 1).sum()) + n_neg = len(neg_scores) + if n_pos == 0 or n_neg == 0: + return None + k = int(np.floor(target_far * n_neg)) + threshold = neg_scores[k - 1] if k > 0 else neg_scores[0] + 1e-9 + predicted_positive = y_score >= threshold + tp = int(np.sum(predicted_positive & (y_true == 1))) + return tp / n_pos + + +def _p_warning_series(head, reg, X: pd.DataFrame, threshold: float, sigma: float) -> pd.Series: + """Model score if a classifier head exists, else the sigmoid-derived fallback probability.""" + if head is not None: + return pd.Series(head.predict_proba(X)[:, 1], index=X.index) + predicted_max = pd.Series(reg.predict(X), index=X.index) + return 1.0 / (1.0 + np.exp(-(predicted_max - threshold) / sigma)) + + +def _find_events(observed_level: pd.Series, warn_thr: float) -> List[dict]: + """Group contiguous observed hours >= warn_thr into flood events.""" + above = observed_level >= warn_thr + events: List[dict] = [] + start = None + prev_t = None + for t, is_above in above.items(): + if is_above and start is None: + start = t + elif not is_above and start is not None: + window = observed_level.loc[start:prev_t] + events.append({"crossed_warn_at": start, "peak_time": window.idxmax(), "peak_level": float(window.max())}) + start = None + prev_t = t + if start is not None: + window = observed_level.loc[start:] + events.append({"crossed_warn_at": start, "peak_time": window.idxmax(), "peak_level": float(window.max())}) + return events + + +def _first_alert_at(p_series: pd.Series, crossed_at, lookback_h: int = 48): + """Earliest time p_warning was sustained (>=0.5 for 2 consecutive hours) within the prior lookback_h.""" + window = p_series.loc[crossed_at - pd.Timedelta(hours=lookback_h) : crossed_at] + sustained = (window >= 0.5) & (window.shift(1) >= 0.5) + hits = sustained[sustained].index + if len(hits) == 0: + return None + return hits.min() - pd.Timedelta(hours=1) + + +def _events_with_lead_time( + observed_level_test: pd.Series, warn_thr: float, p_warning_test: pd.Series +) -> List[dict]: + events = _find_events(observed_level_test, warn_thr) + for event in events: + first_alert_at = _first_alert_at(p_warning_test, event["crossed_warn_at"]) + event["first_alert_at"] = first_alert_at.isoformat() if first_alert_at is not None else None + if first_alert_at is not None: + lead_hours = (event["crossed_warn_at"] - first_alert_at).total_seconds() / 3600.0 + else: + lead_hours = None + event["lead_hours"] = lead_hours + event["crossed_warn_at"] = event["crossed_warn_at"].isoformat() + event["peak_time"] = event["peak_time"].isoformat() + return events + + +def train_station( + df_long: pd.DataFrame, + station: str, + horizons: Tuple[int, ...] = HORIZONS, + skip_eval: bool = False, + hgb_overrides: Optional[dict] = None, + split_train_end: str = SPLIT_B_TRAIN_END, + split_test_start: str = SPLIT_B_TEST_START, + split_test_end: str = SPLIT_B_TEST_END, +) -> 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) + if meta["n_rows"] < MIN_ROWS_TO_TRAIN: + return None, {"status": "failed", "reason": f"only {meta['n_rows']} usable rows (< {MIN_ROWS_TO_TRAIN})"} + + warn_thr, danger_thr = features.get_thresholds(station) + feature_names = list(X.columns) + + if skip_eval: + train_mask = pd.Series(True, index=X.index) + test_mask = pd.Series(False, index=X.index) + else: + train_mask = X.index <= pd.Timestamp(split_train_end) + test_mask = (X.index >= pd.Timestamp(split_test_start)) & (X.index <= pd.Timestamp(split_test_end)) + X_train, Y_train = X.loc[train_mask], Y.loc[train_mask] + X_test, Y_test = X.loc[test_mask], Y.loc[test_mask] + eval_X, eval_Y = (X, Y) if skip_eval else (X_train, Y_train) + + heads: Dict[str, object] = {} + sigma: Dict[int, float] = {} + skipped_heads: Dict[str, str] = {} + per_horizon: Dict[int, dict] = {} + observed_grid = features.make_hourly_grid(df_long).observed + + for h in horizons: + max_col, warn_col, danger_col = f"max_level_{h}", f"exceed_warn_{h}", f"exceed_danger_{h}" + horizon_metrics: dict = {} + + # --- regression head (max level) --- + reg_labeled = eval_Y[max_col].notna() + reg = None + if reg_labeled.sum() >= MIN_ROWS_FOR_HEAD: + reg = _safe_fit( + _make_regressor(hgb_overrides), + eval_X.loc[reg_labeled], + eval_Y.loc[reg_labeled, max_col], + f"max_{h}", + skipped_heads, + ) + else: + skipped_heads[f"max_{h}"] = f"only {int(reg_labeled.sum())} labeled rows" + + sigma_h = MIN_SIGMA + if reg is not None and not skip_eval: + test_labeled = Y_test[max_col].notna() + if test_labeled.sum() > 0: + y_true = Y_test.loc[test_labeled, max_col] + y_pred = reg.predict(X_test.loc[test_labeled]) + residuals = y_true.to_numpy() - y_pred + sigma_h = max(float(np.std(residuals)), MIN_SIGMA) + horizon_metrics["n_test"] = int(test_labeled.sum()) + horizon_metrics["mae"] = float(mean_absolute_error(y_true, y_pred)) + horizon_metrics["rmse"] = float(np.sqrt(mean_squared_error(y_true, y_pred))) + above_2m = y_true >= 2.0 + horizon_metrics["mae_above_2m"] = ( + float(mean_absolute_error(y_true[above_2m], y_pred[above_2m])) if above_2m.any() else None + ) + sigma[h] = sigma_h + horizon_metrics["sigma"] = sigma_h + + # --- classification heads (warn / danger) --- + p_warning_test = None + for label_name, col, thr in (("warn", warn_col, warn_thr), ("danger", danger_col, danger_thr)): + train_labeled = eval_Y[col].notna() + n_pos = int(eval_Y.loc[train_labeled, col].sum()) if train_labeled.any() else 0 + head_key = f"{label_name}_{h}" + clf = None + if n_pos >= MIN_POSITIVES_FOR_CLASSIFIER: + clf = _safe_fit( + _make_classifier(hgb_overrides), + eval_X.loc[train_labeled], + eval_Y.loc[train_labeled, col], + head_key, + skipped_heads, + ) + else: + skipped_heads[head_key] = f"only {n_pos} positives in train span (< {MIN_POSITIVES_FOR_CLASSIFIER})" + heads[head_key] = clf + + if not skip_eval: + test_labeled = Y_test[col].notna() + horizon_metrics[f"base_rate_{label_name}"] = ( + float(Y_test.loc[test_labeled, col].mean()) if test_labeled.any() else None + ) + if clf is not None and test_labeled.sum() > 0 and Y_test.loc[test_labeled, col].nunique() > 1: + y_true = Y_test.loc[test_labeled, col] + y_score = clf.predict_proba(X_test.loc[test_labeled])[:, 1] + horizon_metrics[f"pr_auc_{label_name}"] = float(average_precision_score(y_true, y_score)) + horizon_metrics[f"brier_{label_name}"] = float(brier_score_loss(y_true, y_score)) + horizon_metrics[f"recall_{label_name}_at_far1pct"] = _recall_at_far(y_true, y_score, 0.01) + horizon_metrics[f"recall_{label_name}_at_far5pct"] = _recall_at_far(y_true, y_score, 0.05) + else: + horizon_metrics[f"pr_auc_{label_name}"] = None + horizon_metrics[f"brier_{label_name}"] = None + horizon_metrics[f"recall_{label_name}_at_far1pct"] = None + horizon_metrics[f"recall_{label_name}_at_far5pct"] = None + + if label_name == "warn" and not skip_eval and reg is not None: + p_warning_test = _p_warning_series(clf, reg, X_test, thr, sigma_h) + + per_horizon[h] = horizon_metrics + heads[f"max_{h}"] = reg + + if not skip_eval and reg is not None and p_warning_test is not None: + observed_test_level = observed_grid.get((station, "water_level")) + if observed_test_level is not None: + observed_test_level = observed_test_level.loc[observed_test_level.index.isin(X_test.index)] + per_horizon[h]["events"] = _events_with_lead_time(observed_test_level, warn_thr, p_warning_test) + + # --- full refit on the ENTIRE record for the deployed artifact --- + # This may include/exclude different heads than the eval-phase gate above (the + # full record has more labeled rows), so skip reasons are re-derived here -- + # skipped_heads must reflect what actually ends up in the saved bundle. + final_heads: Dict[str, object] = {} + for h in horizons: + max_col, warn_col, danger_col = f"max_level_{h}", f"exceed_warn_{h}", f"exceed_danger_{h}" + head_key = f"max_{h}" + labeled = Y[max_col].notna() + if labeled.sum() >= MIN_ROWS_FOR_HEAD: + reg = _safe_fit(_make_regressor(hgb_overrides), X.loc[labeled], Y.loc[labeled, max_col], head_key, skipped_heads) + final_heads[head_key] = reg + if reg is not None: + skipped_heads.pop(head_key, None) + else: + skipped_heads[head_key] = f"only {int(labeled.sum())} labeled rows" + final_heads[head_key] = None + + for label_name, col in (("warn", warn_col), ("danger", danger_col)): + head_key = f"{label_name}_{h}" + train_labeled = Y[col].notna() + n_pos = int(Y.loc[train_labeled, col].sum()) if train_labeled.any() else 0 + if n_pos >= MIN_POSITIVES_FOR_CLASSIFIER: + clf = _safe_fit( + _make_classifier(hgb_overrides), X.loc[train_labeled], Y.loc[train_labeled, col], head_key, skipped_heads + ) + final_heads[head_key] = clf + if clf is not None: + skipped_heads.pop(head_key, None) + else: + skipped_heads[head_key] = f"only {n_pos} positives in train span (< {MIN_POSITIVES_FOR_CLASSIFIER})" + final_heads[head_key] = None + + bundle = { + "station_code": station, + "model_version": f"hgb-v1+{_git_short_sha()}", + "trained_at": datetime.datetime.now().isoformat(), + "sklearn_version": sklearn.__version__, + "feature_names": feature_names, + "horizons": list(horizons), + "thresholds": {"warning": warn_thr, "danger": danger_thr}, + "heads": final_heads, + "sigma": sigma, + "skipped_heads": skipped_heads, + "train_span": meta["span"], + "n_train_rows": meta["n_rows"], + } + station_metrics = {"status": "trained", "per_horizon": per_horizon} + return bundle, station_metrics + + +def train_all( + df_long: pd.DataFrame, + stations: List[str], + horizons: Tuple[int, ...] = HORIZONS, + models_dir: Path = Path("models"), + skip_eval: bool = False, + hgb_overrides: Optional[dict] = None, +) -> dict: + """Train and save every requested station's models. Returns the metrics.json payload.""" + models_dir = Path(models_dir) + models_dir.mkdir(parents=True, exist_ok=True) + model_version = f"hgb-v1+{_git_short_sha()}" + + station_results: Dict[str, dict] = {} + for station in stations: + if station in features.NOT_TRAINABLE: + reason = features.NOT_TRAINABLE[station] + logger.info(f"{station}: heuristic ({reason})") + station_results[station] = {"status": "heuristic", "reason": reason} + continue + try: + bundle, station_metrics = train_station( + df_long, station, horizons, skip_eval=skip_eval, hgb_overrides=hgb_overrides + ) + if bundle is None: + logger.warning(f"{station}: failed ({station_metrics.get('reason')})") + station_results[station] = station_metrics + continue + joblib.dump(bundle, models_dir / f"flood_{station}.joblib") + logger.info( + f"{station}: trained, {bundle['n_train_rows']} rows, " + f"{len(bundle['skipped_heads'])} heads skipped" + ) + station_results[station] = station_metrics + except Exception as error: + logger.error(f"{station}: failed with exception: {error}") + station_results[station] = {"status": "failed", "reason": str(error)} + + metrics_payload = { + "generated_at": datetime.datetime.now().isoformat(), + "model_version": model_version, + "split": { + "train_end": SPLIT_B_TRAIN_END, + "test_start": SPLIT_B_TEST_START, + "test_end": SPLIT_B_TEST_END, + }, + "stations": station_results, + } + with open(models_dir / "metrics.json", "w", encoding="utf-8") as handle: + json.dump(metrics_payload, handle, indent=2, default=str) + return metrics_payload + + +def main(argv: Optional[List[str]] = None) -> None: + logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(name)s: %(message)s") + + parser = argparse.ArgumentParser(description="Train Ping River flood forecast models") + parser.add_argument("--stations", default="all", help="'all' or a comma-separated list of station codes") + parser.add_argument("--models-dir", default="models") + parser.add_argument("--db-url", default=None) + parser.add_argument("--api-url", default=DEFAULT_API_URL) + parser.add_argument("--skip-eval", action="store_true", help="Refit-only fast path; skip Split B evaluation") + parser.add_argument("--start", default=None, help="ISO date; earliest measurement to load") + parser.add_argument("--end", default=None, help="ISO date; latest measurement to load") + args = parser.parse_args(argv) + + if args.stations == "all": + stations = list(features.UPSTREAM_LEADS.keys()) + else: + stations = [s.strip() for s in args.stations.split(",") if s.strip()] + + start = datetime.datetime.fromisoformat(args.start) if args.start else None + end = datetime.datetime.fromisoformat(args.end) if args.end else None + + logger.info(f"Loading measurements for {len(stations)} stations...") + df_long = load_measurements( + db_url=resolve_db_url(args.db_url), stations=None, start=start, end=end, api_url=args.api_url + ) + logger.info(f"Loaded {len(df_long)} rows spanning {df_long['timestamp'].min()} .. {df_long['timestamp'].max()}") + + metrics_payload = train_all( + df_long, stations, models_dir=Path(args.models_dir), skip_eval=args.skip_eval + ) + trained = sum(1 for s in metrics_payload["stations"].values() if s["status"] == "trained") + logger.info(f"Done: {trained}/{len(stations)} stations trained. metrics.json written to {args.models_dir}") + + +if __name__ == "__main__": + main() diff --git a/src/static/dashboard.html b/src/static/dashboard.html index d65e291..ee00b64 100644 --- a/src/static/dashboard.html +++ b/src/static/dashboard.html @@ -123,6 +123,13 @@ @keyframes riverMove { to { stroke-dashoffset: -40; } } @media (prefers-reduced-motion: reduce) { .flow-line, .flow-marker::before { animation: none; } } .line-swatch { width: 24px; height: 4px; border-radius: 2px; flex: none; } + .forecast-grid { display: grid; grid-template-columns: repeat(auto-fill, minmax(215px, 1fr)); gap: 10px; margin-top: 14px; } + .forecast-station { border: 1px solid var(--border); border-radius: 12px; padding: 10px 12px; } + .forecast-station strong { font-size: .8rem; } + .forecast-station .fc-name { color: var(--muted); font-size: .68rem; margin: 2px 0 8px; overflow: hidden; text-overflow: ellipsis; white-space: nowrap; } + .risk-chips { display: flex; gap: 6px; } + .risk-chip { flex: 1; text-align: center; border-radius: 8px; padding: 5px 4px; font-size: .64rem; font-weight: 800; color: white; } + .risk-chip span { display: block; font-weight: 650; font-size: .58rem; opacity: .85; } .leaflet-popup-content-wrapper { border-radius: 14px; box-shadow: 0 12px 35px rgba(14,45,54,.2); } .popup { min-width: 190px; } .popup-code { font-size: .7rem; color: var(--river); font-weight: 850; text-transform: uppercase; letter-spacing: .08em; } @@ -201,6 +208,14 @@ + +

PostgreSQL history

Select a RID flow station to load the last 7 days

@@ -556,6 +571,7 @@ renderThaiWaterSensors(thaiWaterSensors, new Set(stations.map((station) => station.station_code))); renderSummary(stations, readings); $('loading').style.display = 'none'; + loadForecasts(); // non-blocking; the card stays hidden until models are deployed } catch (error) { $('loading').style.display = 'none'; $('error').style.display = 'grid'; @@ -566,6 +582,57 @@ } } + function riskColor(pWarning, pDanger) { + if (pDanger >= .5) return '#cc4b37'; + if (pWarning >= .5 || pDanger >= .2) return '#d99018'; + if (pWarning >= .2) return '#0a91b9'; + return '#1e8b60'; + } + + async function loadForecasts() { + const card = $('forecast-card'); + try { + const response = await fetch('/forecast'); + if (!response.ok) { card.style.display = 'none'; return; } + const rows = await response.json(); + if (!Array.isArray(rows) || !rows.length) { card.style.display = 'none'; return; } + const byStation = new Map(); + rows.forEach((row) => { + if (!byStation.has(row.station_code)) byStation.set(row.station_code, []); + byStation.get(row.station_code).push(row); + }); + const grid = $('forecast-grid'); + grid.replaceChildren(); + const stations = [...byStation.entries()].map(([code, list]) => ({ + code, list: list.sort((a, b) => a.horizon_hours - b.horizon_hours), + worst: Math.max(...list.map((r) => Math.max(r.p_danger ?? 0, (r.p_warning ?? 0) * .5))) + })).sort((a, b) => b.worst - a.worst); + stations.forEach(({ code, list }) => { + const first = list[0]; + const cardEl = document.createElement('div'); + cardEl.className = 'forecast-station'; + const chips = list.map((r) => { + const pw = r.p_warning ?? 0, pd = r.p_danger ?? 0; + const pct = Math.round(Math.max(pw, pd) * 100); + const title = `+${r.horizon_hours}h · warning ${Math.round(pw * 100)}% · danger ${Math.round(pd * 100)}%` + + (r.predicted_max_level == null ? '' : ` · peak ~${Number(r.predicted_max_level).toFixed(2)} m`) + + (r.source === 'heuristic' ? ' · heuristic fallback' : ''); + return `
${pct}%${r.horizon_hours}h
`; + }).join(''); + cardEl.innerHTML = `${escapeHtml(code)}` + + `
${first.current_level == null ? '' : `now ${Number(first.current_level).toFixed(2)} m · `}peak risk next 24h
` + + `
${chips}
`; + grid.appendChild(cardEl); + }); + const asOf = rows[0].as_of ? new Date(rows[0].as_of).toLocaleString([], { dateStyle: 'medium', timeStyle: 'short' }) : null; + $('forecast-status').textContent = `Probability of exceeding warning / danger level within 6, 12 and 24 h` + + (asOf ? ` · based on readings up to ${asOf}` : ''); + card.style.display = 'block'; + } catch (error) { + card.style.display = 'none'; + } + } + $('refresh-button').addEventListener('click', loadDashboard); $('history-range').addEventListener('change', () => { if (state.selectedStation) loadHistory(state.selectedStation); }); loadDashboard(); diff --git a/src/web_api.py b/src/web_api.py index 8e08bde..339961b 100644 --- a/src/web_api.py +++ b/src/web_api.py @@ -41,6 +41,10 @@ HISTORY_CACHE: Dict[str, tuple] = {} HISTORY_CACHE_LOCK = Lock() HISTORY_TTL = 300 # 5 minutes +FORECAST_CACHE: Dict[str, tuple] = {} +FORECAST_CACHE_LOCK = Lock() +FORECAST_TTL = 900 # 15 minutes + # Dashboard HTML is loaded once at import from src/static/dashboard.html. _DASHBOARD_HTML_PATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), "static", "dashboard.html") try: @@ -495,6 +499,33 @@ async def get_postgres_history( raise HTTPException(status_code=502, detail="Measurement history unavailable") +@app.get("/forecast") +async def get_flood_forecasts(): + """Flood-risk forecasts per station for the 6/12/24 h horizons.""" + increment_counter("api_requests", labels={"endpoint": "forecast"}) + now = time.monotonic() + with FORECAST_CACHE_LOCK: + cached = FORECAST_CACHE.get("all") + if cached and now - cached[0] < FORECAST_TTL: + return cached[1] + try: + from .ml.predict import get_latest_forecasts + except ImportError as error: + raise HTTPException(status_code=503, detail=f"Forecasting unavailable: {error}") + try: + data = await asyncio.to_thread(get_latest_forecasts) + except FileNotFoundError: + raise HTTPException(status_code=503, detail="No trained flood models found") + except RuntimeError as error: + raise HTTPException(status_code=503, detail=str(error)) + except Exception as error: + logger.error(f"Error computing flood forecasts: {error}") + raise HTTPException(status_code=502, detail="Flood forecast unavailable") + with FORECAST_CACHE_LOCK: + FORECAST_CACHE["all"] = (now, data) + return data + + @app.get("/measurements/latest", response_model=List[MeasurementResponse]) async def get_latest_measurements(limit: int = 100): """Get latest measurements from all stations""" diff --git a/tests/test_flood_forecast.py b/tests/test_flood_forecast.py new file mode 100644 index 0000000..367098b --- /dev/null +++ b/tests/test_flood_forecast.py @@ -0,0 +1,238 @@ +"""Tests for the flood forecast ML package. Synthetic data only -- no DB/network.""" + +import datetime +from typing import Dict, List, Optional + +import joblib +import numpy as np +import pandas as pd +import pytest + +from src.ml import features, predict, train + + +def make_synth( + n_hours: int, + stations: List[str], + seed: int = 0, + start: str = "2020-01-01", + pulses: Optional[Dict[str, List[tuple]]] = None, + missing_patches: Optional[Dict[str, List[tuple]]] = None, +) -> pd.DataFrame: + """Generate a synthetic long measurement frame with smooth levels, flood pulses, + and optional missing patches, for `stations` over `n_hours` hourly steps. + + pulses: {station: [(start_hour, width_hours, peak_add), ...]} + missing_patches: {station: [(start_hour, length_hours), ...]} + """ + rng = np.random.default_rng(seed) + idx = pd.date_range(start, periods=n_hours, freq="h") + rows = [] + for station in stations: + base = 1.5 + 0.1 * np.sin(np.linspace(0, 6 * np.pi, n_hours)) + noise = rng.normal(0, 0.02, n_hours) + level = base + noise + for pulse_start, width, peak_add in (pulses or {}).get(station, []): + t = np.arange(n_hours) + bump = peak_add * np.exp(-0.5 * ((t - (pulse_start + width / 2)) / (width / 4)) ** 2) + level = level + bump + discharge = 20.0 * level + rng.normal(0, 1.0, n_hours) + + missing = np.zeros(n_hours, dtype=bool) + for patch_start, length in (missing_patches or {}).get(station, []): + missing[patch_start : patch_start + length] = True + + for i in range(n_hours): + if missing[i]: + continue + rows.append( + { + "timestamp": idx[i], + "station_code": station, + "water_level": round(float(level[i]), 3), + "discharge": round(float(discharge[i]), 2), + } + ) + return pd.DataFrame(rows) + + +def test_no_future_leakage(): + stations = ["P.1", "P.20"] + df_a = make_synth(200, stations, seed=1, pulses={"P.1": [(150, 10, 3.0)]}) + grid_a = features.make_hourly_grid(df_a) + feat_a = features.build_features(grid_a, "P.1") + + t0 = grid_a.observed.index[120] + + df_b = df_a.copy() + future_mask = df_b["timestamp"] > t0 + df_b.loc[future_mask, "water_level"] = df_b.loc[future_mask, "water_level"] + 50.0 + df_b.loc[future_mask, "discharge"] = df_b.loc[future_mask, "discharge"] + 500.0 + grid_b = features.make_hourly_grid(df_b) + feat_b = features.build_features(grid_b, "P.1") + + past_a = feat_a.loc[feat_a.index <= t0] + past_b = feat_b.loc[feat_b.index <= t0] + pd.testing.assert_frame_equal(past_a, past_b) + + +def test_label_alignment(): + idx = pd.date_range("2020-01-01", periods=12, freq="h") + levels = [1.0, 1.0, 1.0, 1.0, 1.0, 3.5, 3.5, 1.0, 1.0, 1.0, 1.0, 1.0] + df = pd.DataFrame( + { + "timestamp": idx, + "station_code": "P.1", + "water_level": levels, + "discharge": [20.0 * lvl for lvl in levels], + } + ) + grid = features.make_hourly_grid(df) + labels = features.build_labels(grid, "P.1", horizons=(6,)) + + # Level crosses warn (3.0) at t=5. A 6h forward window (t, t+6] first + # includes t=5 for t=0 .. t=4 (inclusive), so exceed_warn_6 should be 1 + # for t=0..4 and not (necessarily) for later rows in this hand-built series. + for t in range(5): + assert labels["exceed_warn_6"].iloc[t] == 1.0, f"t={t} expected warn exceedance" + + # max_level_6 at t=0 covers hours 1..6 -> includes the 3.5 peak. + assert labels["max_level_6"].iloc[0] == pytest.approx(3.5) + + +def test_label_coverage_gate(): + n = 40 + idx = pd.date_range("2020-01-01", periods=n, freq="h") + levels = [1.0] * n + df = pd.DataFrame( + {"timestamp": idx, "station_code": "P.1", "water_level": levels, "discharge": [20.0] * n} + ) + # Drop 70% of a future window (hours 21..26) for the row at t=20, no exceedance in it. + df_missing = df[~df["timestamp"].isin(idx[21:26])].copy() + grid = features.make_hourly_grid(df_missing) + labels = features.build_labels(grid, "P.1", horizons=(6,)) + t20 = idx[20] + assert pd.isna(labels.loc[t20, "exceed_warn_6"]) + + # Same sparse window, but WITH an observed exceedance inside it -> must be 1, not NaN. + df_with_peak = df_missing.copy() + peak_row = pd.DataFrame( + [{"timestamp": idx[22], "station_code": "P.1", "water_level": 5.0, "discharge": 100.0}] + ) + df_with_peak = pd.concat([df_with_peak, peak_row], ignore_index=True) + grid2 = features.make_hourly_grid(df_with_peak) + labels2 = features.build_labels(grid2, "P.1", horizons=(6,)) + assert labels2.loc[t20, "exceed_warn_6"] == 1.0 + + +def test_ffill_and_staleness(): + n = 20 + idx = pd.date_range("2020-01-01", periods=n, freq="h") + df = pd.DataFrame( + { + "timestamp": idx, + "station_code": "P.1", + "water_level": [1.0 + 0.01 * i for i in range(n)], + "discharge": [20.0] * n, + } + ) + # Small gap: drop hours 5,6 (2h gap). + df_small_gap = df[~df["timestamp"].isin(idx[5:7])].copy() + grid = features.make_hourly_grid(df_small_gap) + feat = features.build_features(grid, "P.1") + assert feat.loc[idx[5], "obs_age_h"] == pytest.approx(1.0) + assert feat.loc[idx[6], "obs_age_h"] == pytest.approx(2.0) + + # Large gap: drop hours 5..9 (5h gap) -> rows with age>3 dropped (NaN). + df_big_gap = df[~df["timestamp"].isin(idx[5:10])].copy() + grid2 = features.make_hourly_grid(df_big_gap) + feat2 = features.build_features(grid2, "P.1") + assert feat2.loc[idx[8], "obs_age_h"] != feat2.loc[idx[8], "obs_age_h"] # NaN + assert feat2.loc[idx[9], "obs_age_h"] != feat2.loc[idx[9], "obs_age_h"] # NaN + assert feat2.loc[idx[7], "obs_age_h"] == pytest.approx(3.0) + + +_FORECAST_KEYS = { + "station_code", + "horizon_hours", + "p_warning", + "p_danger", + "predicted_max_level", + "current_level", + "as_of", + "model_version", + "trained_at", + "source", + "threshold_warning", + "threshold_danger", +} + + +def _assert_valid_forecast_row(row: dict) -> None: + assert set(row.keys()) == _FORECAST_KEYS + assert 0.0 <= row["p_warning"] <= 1.0 + assert 0.0 <= row["p_danger"] <= 1.0 + assert row["p_danger"] <= row["p_warning"] + assert row["predicted_max_level"] >= row["current_level"] + assert row["source"] in ("model", "heuristic") + + +def test_train_smoke_and_roundtrip(tmp_path): + # Include every station P.1's feature set actually references (its UPSTREAM_LEADS) + # so no upstream column is entirely NaN -- HistGradientBoosting's binning step + # cannot fit a fully-degenerate column (see train._safe_fit). + upstream = [code for code, _lead in features.UPSTREAM_LEADS["P.1"]] + data_stations = ["P.1"] + upstream + target_stations = ["P.1", "P.20"] + n = 700 + pulses = {station: [(start, 20, 2.0) for start in range(50, n - 50, 110)] for station in data_stations} + 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} + ) + assert metrics["stations"]["P.1"]["status"] == "trained" + assert metrics["stations"]["P.20"]["status"] == "trained" + assert (tmp_path / "flood_P.1.joblib").exists() + assert (tmp_path / "metrics.json").exists() + + readings_by_station = { + code: group[["timestamp", "water_level", "discharge"]].to_dict("records") + for code, group in df.groupby("station_code") + if code in target_stations + } + now = df["timestamp"].max() + forecasts = predict.get_forecasts(readings_by_station, models_dir=tmp_path, now=now) + + assert len(forecasts) > 0 + for row in forecasts: + _assert_valid_forecast_row(row) + assert any(row["source"] == "model" for row in forecasts) + + +def test_heuristic_fallback(tmp_path): + df = make_synth(50, ["P.1"], seed=3) + readings_by_station = {"P.1": df[["timestamp", "water_level", "discharge"]].to_dict("records")} + now = df["timestamp"].max() + + forecasts = predict.get_forecasts(readings_by_station, models_dir=tmp_path, now=now) + + assert len(forecasts) == len(predict.DEFAULT_HORIZONS) + for row in forecasts: + _assert_valid_forecast_row(row) + assert row["source"] == "heuristic" + assert row["model_version"] == "heuristic-v1" + assert row["trained_at"] is None + + +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}) + # 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") + + grid = features.make_hourly_grid(df) + fresh_columns = list(features.build_features(grid, "P.1").columns) + assert fresh_columns == bundle["feature_names"]