Northern-Thailand-Ping-River-Monitor/docs/FLOOD_FORECASTING.md

# 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+<git short SHA>`), `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.

![Observed P.1 level vs the model's 24 h-ahead predicted peak through the
October 2024 flood, with the warning probability below](img/backtest-2024-p1.png)

The render above shows the whole event hour by hour. Top: the observed level
(blue) against the 24 h-ahead predicted peak the model issued at each hour
(amber, dashed) — the amber line leads the blue one into both flood waves,
which is the entire point of the system. Bottom: the model's probability of
flooding within 24 h; it fires its **first alert at 01:00 on 24 September, a
full day before the river crossed the 3.70 m flooding line**, stays pinned near
1.0 through both waves, and stands down between and after them. Also visible,
honestly: the predicted peak tops out ~0.4 m short of the actual 5.30 m record
(the extreme-compression limitation discussed below), and the prediction is
noisier on the recession limbs. The same model track drives the dashboard's
"Replay Oct 2024 flood" feature, so this chart can be watched live on the map.

- **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 per-station as of 2026-08-10.** `THRESHOLDS` now carries
calibrated (warning, danger) pairs for all 16 stations, derived from the DB's
`discharge_percent` (RID % of channel capacity): warning = median level at
75–85% capacity, danger = median level at 95–105%. P.1 instead uses the official
Chiang Mai inundation map (`P1_FLOOD_STAGES`): warning 3.70 m (city flooding
begins, stage 1) and danger 4.20 m (stage 5). The prior single default of
(3.0, 4.5) m made P.103 badly over-alert (its bank-full level is ~6.75 m) and
P.67 under-alert (overflow at ~2.9 m, 1.6 m below the old danger line).
**A retrain is required after any threshold change** — classifier labels depend
on them; until then, model rows report the thresholds baked into their bundle.
P.1 additionally reports `stages`: exceedance probability for each of the seven
official inundation stages (3.70–4.60 m), computed from the regression head and
its calibration sigma, so they need no retrain and no per-stage classifiers.

## 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'])"
```