diff --git a/docs/FLOOD_FORECASTING.md b/docs/FLOOD_FORECASTING.md index e798865..12470fb 100644 --- a/docs/FLOOD_FORECASTING.md +++ b/docs/FLOOD_FORECASTING.md @@ -203,16 +203,17 @@ section 6). 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. +1. **Belt-and-braces probability** *(since 2026-08-11 — see the re-examination + note in section 7)*: the sigmoid-of-regression probability + `p = 1/(1 + exp(−(predicted_max − threshold)/σ))` is always computed (σ = + the regressor's test-residual std, floor `MIN_SIGMA = 0.15` m), and when a + classifier head exists — trained only if the span had at least + `MIN_POSITIVES_FOR_CLASSIFIER = 30` positives — the served probability is + `max(classifier, sigmoid)`. The classifier can raise the alarm but never + silence it: on the gap-filled data a trained classifier stayed near zero + through the 2024 record crossing while the regression tracked it. +2. **Sigmoid only**, when the classifier head is absent or skipped + (recorded in `skipped_heads` with its reason). 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 @@ -255,47 +256,92 @@ invalidates the cache without a restart. ### 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. +Two evaluations exist and they differ sharply — the re-examination note in the +next section explains why (the hourly grid was gap-filled from ~56% to ~93% +between them, roughly doubling the test rows and collapsing the warning base +rates). -P.1 (Nawarat Bridge), the station that matters most: +**Current model** `hgb-v1+d2d0e65`, generated 2026-08-12 on the gap-filled DB +(~976k rows). Train ≤ 2024-12-31, test 2025-01-01 → 2026-08-12. P.1: -| 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% | +| Horizon | Warning PR-AUC | MAE | MAE above 2 m | Test rows | Base rate | +|---|---|---|---|---|---| +| 6 h | 0.783 | 5.5 cm | 7.9 cm | 14,034 | 0.12% | +| 12 h | 0.508 | 8.1 cm | 14.6 cm | 14,028 | 0.16% | +| 24 h | 0.288 | 10.5 cm | 24.0 cm | 14,034 | 0.25% | -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. +Level accuracy improved; standalone classifier discrimination did not survive +the data change (which is why serving is now `max(classifier, sigmoid)` — see +"Head gating"). Recall-at-FAR is null at all horizons on this run. Across +stations the 6 h warning PR-AUC now spans 0.987 (P.77) / 0.982 (P.5) / 0.956 +(P.85) / 0.950 (P.67) down to 0.436 (P.84), and danger heads are now evaluable +at nine stations — strongest P.5 (0.958/0.883/0.811 at 6/12/24 h) and P.77 +(0.942/0.863/0.786); P.103's danger metrics, previously the highlight, are null +on this span. -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. +**Historical evaluation** (`hgb-v1+49a3de0`, 2026-08-10, pre-gap-fill DB — +kept for the record; these numbers described the sparser 56%-filled grid and do +not reproduce on today's data): -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. +| Horizon | Warning PR-AUC | Recall @1% FAR | MAE | Test rows | Base rate | +|---|---|---|---|---|---| +| 6 h | 0.974 | 98.3% | 6.1 cm | 8,536 | 1.36% | +| 12 h | 0.904 | 93.8% | 9.0 cm | 7,932 | 1.61% | +| 24 h | 0.900 | 90.1% | 11.3 cm | 8,572 | 1.77% | + +The dramatic PR-AUC difference is mostly the base rate: the filled grid adds +~5,500 quiet test hours per horizon while the number of positive hours barely +changes, so the same ranking quality scores far lower — and the classifier's +genuine out-of-distribution weakness (see the backtest sections) does the rest. + +### 2026-08-11 re-examination: fuller data changes the backtest story + +> **Read this before the two backtest sections below.** On 2026-08-11 the +> backtests were codified into `scripts/backtest_render.py` (previously they +> were one-off runs) and re-run after the database grew from 592k to ~976k +> rows (a `--fill-gaps all` pass repaired most of the missing 44% of the +> hourly grid). Three things changed: +> +> 1. **The 2024 crossing was 8 hours earlier than documented.** The recovered +> hours show P.1 crossing 3.70 m at **17:00 on 24 September 2024**, not +> 01:00 on 25 September — confirmed independently by the HII sensor at +> Nawarat Bridge (hii_waterlevel, station 3226: 3.73 m at 17:00). The +> originally celebrated "24-hour warning" was therefore ~16 hours measured +> against the real river. +> 2. **Retraining on the fuller data improves level accuracy but degrades the +> warning classifiers.** P.1 24 h MAE improved (11.3 → 10.5 cm), but the +> warning-head PR-AUC collapsed (0.900 → 0.288 at 24 h): with the filled +> grid the classifier trains on many more dry-season rows and now stays +> silent through the September 2024 record crossing while the regression +> head tracks it. Serving was changed to belt-and-braces — +> `max(classifier, sigmoid(regression))` — so alerting can never be worse +> than the regression path. +> 3. **Honest current lead times, from the regenerated charts below:** the +> retrained configuration first alerts ~18 h *after* the true 24 Sep 2024 +> crossing and roughly *at* the 27 Sep 2025 crossing. The earlier, better +> numbers came from models trained and evaluated on the sparser data. The +> conclusion is not that the old system was better — it is that gauge-only +> features fundamentally lack lead time for fast rises, which is exactly +> the rainfall-input and rise-target work now queued (see "Honest limits"). +> +> `scripts/backtest_render.py` regenerates all three charts and fails its +> acceptance gate while the 2024 lead stays under 12 h — keeping this page +> honest is now automatic. ### The September 2025 flood, as the deployed configuration saw it ![Observed vs predicted through the September 2025 flood — model trained only through 2024](img/backtest-2025-p1.png) -This is the strongest single piece of evidence, because it uses the exact -deployed configuration (train ≤ 2024-12-31) on an event it never saw: the -model's **first alert came 26 September 2025 at 18:00, exactly 24 hours before -the river crossed 3.70 m**, and it predicted a 4.00 m peak against an actual -3.93 m — within 7 cm. Note the discrimination: the near-miss crest of 3.51 m on -26 September never triggered an alert, the probability fires only for the real -event, and it stands down as the water recedes. +This uses the deployed configuration (train ≤ 2024-12-31) on an event it never +saw. *(Chart regenerated 2026-08-11 on the gap-filled data — see the +re-examination note above; the original one-off render, trained on the sparser +data, alerted 24 h ahead and predicted the 3.93 m peak within 7 cm.)* On +today's fuller dataset the retrained equivalent first alerts at **18:00 on +27 September 2025 — as the river crosses 3.70 m**, not a day ahead. The +discrimination remains good: the near-miss 3.51 m crest on 26 September never +triggers, the probability fires only for the real event, and it stands down as +the water recedes. ### Headline validation: the October 2024 record flood @@ -307,27 +353,33 @@ 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. +The render above shows the whole event hour by hour *(regenerated 2026-08-11 +on the gap-filled data)*. Top: the observed level (blue) against the 24 h-ahead +predicted peak the model issued at each hour (amber, dashed) — the amber line +tracks both flood waves but no longer clearly leads the first one. Bottom: the +belt-and-braces probability of flooding within 24 h; on the fuller data its +**first alert comes at 11:00 on 25 September, ~18 hours after the true 17:00 +24 September crossing**, then stays correctly alarmed through the October +record wave. Also visible, honestly: the predicted peak tops out well short of +the actual 5.30 m record (the extreme-compression limitation discussed below). +The same historic model track drives the dashboard's "Replay Oct 2024 flood" +feature. ![Hour-by-hour detail of the detection window, 22–28 September 2024](img/backtest-2024-p1-detail.png) -The hour-by-hour detail of the detection window shows the sequence exactly: the -predicted 24 h peak (amber) starts pulling away from the observed level late on -23 September as upstream gauges rise, the warning probability snaps from ~0 to -1.0 at **01:00 on 24 September**, and the river crosses 3.70 m at **01:00 on -25 September** — a clean 24-hour warning, delivered while the river in town -still looked normal at 2.8 m. +The hour-by-hour detail of the detection window *(regenerated 2026-08-11)* +shows the corrected sequence: the river crosses 3.70 m at **17:00 on +24 September** (the hours recovered by gap-filling; independently confirmed by +the HII sensor at the same bridge), while the retrained model's probability +only crosses 0.5 at **11:00 on 25 September**. The original render — sparser +data, different trained model — alerted at 01:00 on 24 September against an +apparent 01:00 25 September crossing. Closing this real gap is what the +rainfall features and rise-target work are for. + +The event bullets below quote the original (pre-gap-fill) evaluation of the +deployed model and are kept for the historical record — see the re-examination +note above for why the lead times no longer reproduce: - **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 @@ -353,11 +405,14 @@ still looked normal at 2.8 m. 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. +is extrapolating basin state and season, not routing a wave. The regenerated +backtests bear this out — harder than first documented (see the re-examination +note in section 7): on the gap-filled data the retrained configuration alerts +the 27 September 2025 crossing *as it happens* and the 24 September 2024 +crossing ~18 h *late*. **Genuine gauge-only lead for P.1 is at best ~7–17 h, +and for fast rises can be zero.** Extending it requires rainfall inputs and +Mae Ngat/Mae Kuang dam release data, plus the rise-target/quantile modelling +work — rainfall collection began 2026-08-11 (see `docs/DATA_SOURCES.md`). **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 @@ -617,10 +672,12 @@ print({h: (d.get('pr_auc_warn'), d.get('mae')) for h, d in m['stations']['P.1'][ ``` 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. +`failed` names its reason in the same payload. Compare against the *previous +run's* `metrics.json`, not an absolute bar: after the 2026-08-11 gap-fill the +expected baseline is P.1 6 h warning PR-AUC ≈ 0.78 and MAE ≈ 5.5 cm (the +historical ~0.97 figure belonged to the sparse pre-fill grid — see section 5). +A *material drop from the previous run* 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): diff --git a/docs/img/backtest-2024-p1-detail.png b/docs/img/backtest-2024-p1-detail.png index 8b87f97..e57ff92 100644 Binary files a/docs/img/backtest-2024-p1-detail.png and b/docs/img/backtest-2024-p1-detail.png differ diff --git a/docs/img/backtest-2024-p1.png b/docs/img/backtest-2024-p1.png index 610544d..3bcf1a4 100644 Binary files a/docs/img/backtest-2024-p1.png and b/docs/img/backtest-2024-p1.png differ diff --git a/docs/img/backtest-2025-p1.png b/docs/img/backtest-2025-p1.png index e607a0a..ea7ab86 100644 Binary files a/docs/img/backtest-2025-p1.png and b/docs/img/backtest-2025-p1.png differ diff --git a/scripts/backtest_render.py b/scripts/backtest_render.py new file mode 100644 index 0000000..54a6b3d --- /dev/null +++ b/scripts/backtest_render.py @@ -0,0 +1,253 @@ +#!/usr/bin/env python3 +"""Regenerate the documented P.1 flood-backtest charts in docs/img/. + +For each chart an eval-only model (regression 24 h peak + warning classifier) +is trained on data STRICTLY BEFORE the event, then the event window is walked +hour by hour exactly as the live system would have seen it: + + backtest-2024-p1.png Oct 2024 record flood, trained < 1 Sep 2024 + backtest-2024-p1-detail.png 22-28 Sep 2024 zoom of the first crossing + backtest-2025-p1.png Sep 2025 flood, deployed config (trained <= 2024) + +This codifies the previously prose-only acceptance test: the run fails with a +non-zero exit if the model gives less than 12 h of warning before the first +3.70 m crossing of the 2024 event. + +Usage: + python scripts/backtest_render.py # uses FLOOD_ML_DB_URL/Config + python scripts/backtest_render.py --db-url postgresql://... +""" + +import argparse +import os +import sys + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..")) + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.dates as mdates +import matplotlib.pyplot as plt +import pandas as pd + +from src.ml import data, features +from src.ml.train import _make_classifier, _make_regressor + +STATION = "P.1" +STAGE1 = 3.70 # official Chiang Mai stage 1 - city flooding begins +STAGE7 = 4.60 # stage 7 - widespread +HORIZON = 24 + +INK = "#132b35" +BLUE = "#1c6ea4" +AMBER = "#c07d10" +RED = "#d9534f" + + +def fit_backtest_model(df_long: pd.DataFrame, train_end: str): + """Train the 24 h regression + warning heads on rows <= train_end only.""" + X, Y, _meta = features.build_matrix(df_long, STATION, (HORIZON,)) + train_mask = X.index <= pd.Timestamp(train_end) + X_train, Y_train = X.loc[train_mask], Y.loc[train_mask] + + max_col, warn_col = f"max_level_{HORIZON}", f"exceed_warn_{HORIZON}" + reg_rows = Y_train[max_col].notna() + reg = _make_regressor().fit(X_train.loc[reg_rows], Y_train.loc[reg_rows, max_col]) + warn_rows = Y_train[warn_col].notna() + clf = _make_classifier().fit( + X_train.loc[warn_rows], Y_train.loc[warn_rows, warn_col].astype(int) + ) + return X, reg, clf + + +def event_series(df_long, X, reg, clf, window_start: str, window_end: str): + """Observed level plus the forecasts the model would have issued hourly.""" + grid = features.make_hourly_grid(df_long) + # observed has MultiIndex columns (station_code, field) + observed = grid.observed[(STATION, "water_level")] + observed = observed.loc[window_start:window_end].dropna().astype(float) + + Xw = X.loc[window_start:window_end] + forecasts = pd.DataFrame(index=Xw.index) + forecasts["pred_max"] = reg.predict(Xw) + # Belt-and-braces probability: the classifier OR the regression-sigmoid, + # whichever is more alarmed. The classifier alone proved unreliable on + # out-of-distribution extremes (silent on the 2024 record flood). + import numpy as np + p_clf = clf.predict_proba(Xw)[:, 1] + p_sig = 1.0 / (1.0 + np.exp(-(forecasts["pred_max"] - STAGE1) / 0.15)) + forecasts["p_flood"] = np.maximum(p_clf, p_sig) + + flood_start = observed[observed >= STAGE1].index.min() + alerts = forecasts[forecasts["p_flood"] >= 0.5].index + first_alert = alerts.min() if len(alerts) else None + return observed, forecasts, flood_start, first_alert + + +def _style_axes(ax): + ax.spines[["top", "right"]].set_visible(False) + ax.tick_params(colors=INK, labelsize=11) + ax.grid(axis="y", color="#dfe9e7", linewidth=0.8) + ax.set_axisbelow(True) + + +def render(observed, forecasts, flood_start, first_alert, out_path, *, + title, subtitle, detail=False, show_stage7=False, peak_note=None): + fig, (ax, axp) = plt.subplots( + 2, 1, figsize=(12.6, 7.6), sharex=True, + gridspec_kw={"height_ratios": [2.2, 1], "hspace": 0.12}, + ) + fig.patch.set_facecolor("white") + + marker = dict(marker="o", markersize=3) if detail else {} + ax.plot(observed.index, observed.values, color=BLUE, linewidth=2.2, + label="Observed level" + (" (hourly)" if detail else ""), **marker) + marker = dict(marker="s", markersize=3) if detail else {} + ax.plot(forecasts.index, forecasts["pred_max"], color=AMBER, linewidth=2, + linestyle="--", label="Predicted 24 h peak (issued at that hour)", **marker) + + ax.axhline(STAGE1, color=RED, linewidth=1, alpha=0.65) + ax.annotate(f"{STAGE1:.2f} m · stage 1 · flooding begins", xy=(0.06, STAGE1), + xycoords=("axes fraction", "data"), xytext=(0, 5), + textcoords="offset points", color=RED, fontsize=10.5) + if show_stage7: + ax.axhline(STAGE7, color=RED, linewidth=1, alpha=0.65) + ax.annotate(f"{STAGE7:.2f} m · stage 7 · widespread", xy=(0.06, STAGE7), + xycoords=("axes fraction", "data"), xytext=(0, 5), + textcoords="offset points", color=RED, fontsize=10.5) + + if peak_note: + peak_ts = observed.idxmax() + ax.annotate(peak_note, xy=(peak_ts, observed.max()), + xytext=(12, 10), textcoords="offset points", + color=BLUE, fontsize=11.5, fontweight="bold") + + ax.set_ylabel("P.1 water level (m)", color=INK, fontsize=11.5) + ax.legend(loc="upper left", frameon=False, fontsize=10.5) + _style_axes(ax) + + axp.plot(forecasts.index, forecasts["p_flood"], color=AMBER, linewidth=1.8) + axp.fill_between(forecasts.index, 0, forecasts["p_flood"], + color=AMBER, alpha=0.28) + axp.axhline(0.5, color=INK, linewidth=0.9, linestyle=":", alpha=0.6) + axp.set_ylim(-0.02, 1.1) + axp.set_ylabel(f"P(flooding within {HORIZON} h)", color=INK, fontsize=11.5) + _style_axes(axp) + + if first_alert is not None: + lead_h = None if flood_start is None else \ + int((flood_start - first_alert).total_seconds() // 3600) + lead_txt = "" if lead_h is None else ( + f"\n({lead_h} h before flooding began)" if lead_h >= 0 + else f"\n({-lead_h} h after flooding began)" + ) + if detail and flood_start is not None: + for a in (ax, axp): + a.axvline(first_alert, color=AMBER, linewidth=1.4, alpha=0.85) + a.axvline(flood_start, color=BLUE, linewidth=1.4, alpha=0.85) + # Anchor labels away from each other in chronological order so a + # late alert (alert AFTER crossing) cannot overprint the labels. + events = sorted( + [(first_alert, "model alert", AMBER), (flood_start, "flooding begins", BLUE)] + ) + for (ts, label, color), (offset, align) in zip(events, ((-8, "right"), (8, "left"))): + ax.annotate(f"{label}\n{ts:%d %b %H:%M}", + xy=(ts, observed.min()), xytext=(offset, 18), + textcoords="offset points", ha=align, + color=color, fontsize=11, fontweight="bold") + mid_y = observed.min() + (observed.max() - observed.min()) * 0.28 + ax.annotate("", xy=(flood_start, mid_y), xytext=(first_alert, mid_y), + arrowprops=dict(arrowstyle="<->", color=INK, lw=1.3)) + arrow_label = ( + f"{lead_h} h warning" if lead_h >= 0 else f"alert {-lead_h} h late" + ) + ax.annotate(arrow_label, + xy=(first_alert + (flood_start - first_alert) / 2, mid_y), + xytext=(0, 8), textcoords="offset points", ha="center", + color=INK, fontsize=11.5, fontweight="bold") + else: + axp.annotate(f"first alert · {first_alert:%d %b %H:%M}{lead_txt}", + xy=(first_alert, 0.62), xytext=(10, 0), + textcoords="offset points", color=RED, fontsize=10.5, + bbox=dict(facecolor="white", alpha=0.75, edgecolor="none")) + + locator = mdates.DayLocator(interval=1 if detail else 3) + axp.xaxis.set_major_locator(locator) + axp.xaxis.set_major_formatter(mdates.DateFormatter("%d %b")) + fig.suptitle(f"{title}\n{subtitle}", x=0.07, y=0.985, ha="left", + fontsize=15, color=INK) + fig.subplots_adjust(top=0.885, left=0.07, right=0.97, bottom=0.07) + fig.savefig(out_path, dpi=110) + plt.close(fig) + print(f"wrote {out_path}") + + +def main(argv=None) -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--db-url", default=None) + parser.add_argument("--out-dir", default=os.path.join("docs", "img")) + args = parser.parse_args(argv) + + df = data.load_measurements(db_url=args.db_url) + if df.empty: + print("no measurement data available", file=sys.stderr) + return 1 + os.makedirs(args.out_dir, exist_ok=True) + + # --- October 2024 record flood: trained only on data before 1 Sep 2024 --- + X, reg, clf = fit_backtest_model(df, "2024-08-31") + obs, fc, flood_start, first_alert = event_series( + df, X, reg, clf, "2024-09-10", "2024-10-14 23:00") + peak = float(obs.max()) + render(obs, fc, flood_start, first_alert, + os.path.join(args.out_dir, "backtest-2024-p1.png"), + title="October 2024 flood: what the model saw coming", + subtitle="P.1 Nawarat Bridge — model trained only on data before 1 Sep 2024", + show_stage7=True, peak_note=f"record peak {peak:.2f} m") + + obs_d, fc_d, flood_d, alert_d = event_series( + df, X, reg, clf, "2024-09-21 18:00", "2024-09-28 06:00") + lead_h = None + if alert_d is not None and flood_d is not None: + lead_h = int((flood_d - alert_d).total_seconds() // 3600) + render(obs_d, fc_d, flood_d, alert_d, + os.path.join(args.out_dir, "backtest-2024-p1-detail.png"), + title="Detection in detail: 22–28 September 2024, hour by hour", + subtitle=( + f"the model alerts {lead_h} h before the river crosses the flooding line" + if lead_h is not None and lead_h > 0 + else "model alert vs the river crossing the flooding line" + ), + detail=True) + + # --- September 2025 flood: the deployed configuration (trained <= 2024) --- + X25, reg25, clf25 = fit_backtest_model(df, "2024-12-31") + obs25, fc25, flood25, alert25 = event_series( + df, X25, reg25, clf25, "2025-09-22", "2025-10-02 12:00") + pred_at_alert = float(fc25.loc[alert25:, "pred_max"].iloc[:24].max()) if alert25 is not None else None + note = f"peak {float(obs25.max()):.2f} m" + ( + f" (predicted {pred_at_alert:.2f} m)" if pred_at_alert is not None else "") + render(obs25, fc25, flood25, alert25, + os.path.join(args.out_dir, "backtest-2025-p1.png"), + title="The September 2025 flood — as forecast by the deployed configuration", + subtitle="model trained only on data through 2024; this event was never seen in training", + detail=True, peak_note=note) + + print(f"2024: flooding began {flood_start}, first alert {first_alert}") + print(f"2025: flooding began {flood25}, first alert {alert25}") + + # Acceptance gate: the flagship 2024 event must keep a >= 12 h warning + if first_alert is None or flood_start is None: + print("FAIL: 2024 event alert or crossing not found", file=sys.stderr) + return 1 + lead = (flood_start - first_alert).total_seconds() / 3600 + if lead < 12: + print(f"FAIL: 2024 first-alert lead {lead:.0f} h < 12 h", file=sys.stderr) + return 1 + print(f"PASS: 2024 first-alert lead {lead:.0f} h") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/scripts/load_test.py b/scripts/load_test.py new file mode 100644 index 0000000..e528af7 --- /dev/null +++ b/scripts/load_test.py @@ -0,0 +1,148 @@ +#!/usr/bin/env python3 +"""Staged load / client-stress test for the Ping River Monitor API + dashboard. + +Simulates a realistic traffic mix (dashboard page loads, the API calls the +dashboard itself makes, heavy history queries, external API consumers) at +increasing concurrency stages, and reports throughput, latency percentiles, +and errors per stage plus the slowest endpoints. + +Run against a LOCAL instance for full stress (never full-stress production — +it hosts live flood monitoring): + + python -m uvicorn src.web_api:app --port 8125 # separate shell + python scripts/load_test.py http://localhost:8125 + +A gentle production baseline (low, fixed concurrency): + + python scripts/load_test.py https://water.buildfor.life --gentle +""" + +import argparse +import random +import statistics +import threading +import time +from collections import Counter + +import requests + +# Weighted endpoint mix: dashboard session + API consumers +ENDPOINTS = [ + ("/", 10), + ("/measurements/latest?limit=500", 20), + ("/stations", 10), + ("/api/hii/rainfall/latest", 15), + ("/api/hii/waterlevel/latest", 15), + ("/forecast", 10), + ("/api/stats", 5), + ("/measurements/history/P.1?hours=168", 10), + ("/measurements/history/P.67?hours=720", 5), + ("/health", 5), +] +POOL = [endpoint for endpoint, weight in ENDPOINTS for _ in range(weight)] + +FULL_STAGES = [(10, 20), (50, 20), (200, 25)] # (clients, seconds) +GENTLE_STAGES = [(3, 15), (8, 15)] + + +def _worker(base, stop_at, results, errors): + session = requests.Session() + while time.time() < stop_at: + path = random.choice(POOL) + start = time.perf_counter() + try: + response = session.get(f"{base}{path}", timeout=30) + elapsed = time.perf_counter() - start + if response.status_code == 200: + results.append((path, elapsed)) + else: + errors.append((path, response.status_code)) + except Exception as error: + errors.append((path, type(error).__name__)) + + +def _pct(values, p): + if len(values) >= 100: + return statistics.quantiles(values, n=100)[p - 1] + return max(values) + + +def run_stage(base, clients, seconds): + results, errors = [], [] + stop_at = time.time() + seconds + threads = [ + threading.Thread( + target=_worker, args=(base, stop_at, results, errors), daemon=True + ) + for _ in range(clients) + ] + for thread in threads: + thread.start() + for thread in threads: + thread.join(timeout=seconds + 35) + + latencies = [elapsed for _, elapsed in results] + total = len(results) + len(errors) + print(f"\n== {clients} clients x {seconds}s ==") + print( + f"requests: {total} ok: {len(results)} errors: {len(errors)} " + f"rps: {total / seconds:.1f}" + ) + if latencies: + print( + f"latency ms p50: {statistics.median(latencies) * 1000:.0f} " + f"p95: {_pct(latencies, 95) * 1000:.0f} " + f"p99: {_pct(latencies, 99) * 1000:.0f} " + f"max: {max(latencies) * 1000:.0f}" + ) + by_endpoint = {} + for path, elapsed in results: + by_endpoint.setdefault(path, []).append(elapsed) + slowest = sorted( + by_endpoint.items(), key=lambda kv: -statistics.median(kv[1]) + )[:4] + for path, values in slowest: + print( + f" slow: {path:45} n={len(values):5} " + f"p50={statistics.median(values) * 1000:6.0f}ms " + f"max={max(values) * 1000:7.0f}ms" + ) + if errors: + top = Counter(f"{path} {code}" for path, code in errors).most_common(5) + print(f" errors: {top}") + return {"clients": clients, "total": total, "errors": len(errors)} + + +def main(argv=None) -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("base", nargs="?", default="http://localhost:8125") + parser.add_argument( + "--gentle", + action="store_true", + help="low fixed concurrency (safe for the production instance)", + ) + args = parser.parse_args(argv) + base = args.base.rstrip("/") + + # Warm caches first so stage 1 doesn't measure cold-start work + for path in ("/forecast", "/api/stats", "/measurements/latest?limit=500"): + try: + requests.get(f"{base}{path}", timeout=60) + except Exception: + pass + + print(f"target: {base} mode: {'gentle' if args.gentle else 'full'}") + stages = GENTLE_STAGES if args.gentle else FULL_STAGES + summary = [run_stage(base, clients, seconds) for clients, seconds in stages] + worst = max( + (stage["errors"] / stage["total"] for stage in summary if stage["total"]), + default=1.0, + ) + print(f"\nworst-stage error rate: {worst:.1%}") + return 0 if worst < 0.05 else 1 + + +if __name__ == "__main__": + import sys + + sys.exit(main()) diff --git a/scripts/locustfile.py b/scripts/locustfile.py new file mode 100644 index 0000000..20a0bb1 --- /dev/null +++ b/scripts/locustfile.py @@ -0,0 +1,93 @@ +"""Locust load profile for the Ping River Monitor API + dashboard. + +Two user types mirror real traffic: dashboard visitors (page + the API calls +the page makes, polling like the auto-refresh does) and API consumers +(direct endpoint hits, including heavy history queries). + +Full stress against a LOCAL instance (never full-stress production — it hosts +live flood monitoring): + + # separate shell: python -m uvicorn src.web_api:app --port 8125 + .venv/Scripts/python.exe -m locust -f scripts/locustfile.py \ + --host http://localhost:8125 --headless \ + --users 200 --spawn-rate 20 --run-time 2m \ + --html load-report.html + +Interactive UI instead: drop --headless and open http://localhost:8089. +""" + +import random + +from locust import FastHttpUser, between, task + + +class DashboardVisitor(FastHttpUser): + """A browser session: initial page load, then periodic refresh polling.""" + + weight = 3 + wait_time = between(2, 6) + + def on_start(self): + # What one real page load requests + self.client.get("/") + self.client.get("/stations") + self.client.get("/measurements/latest?limit=500") + self.client.get("/api/hii/waterlevel/latest") + self.client.get("/api/hii/rainfall/latest") + + @task(4) + def poll_latest(self): + self.client.get("/measurements/latest?limit=500") + + @task(2) + def poll_forecast(self): + self.client.get("/forecast") + + @task(2) + def poll_rain(self): + self.client.get("/api/hii/rainfall/latest") + + @task(1) + def view_history(self): + station = random.choice(["P.1", "P.67", "P.103", "P.75", "P.20"]) + hours = random.choice([24, 168, 720]) + self.client.get( + f"/measurements/history/{station}?hours={hours}", + name="/measurements/history/[station]", + ) + + @task(1) + def stats(self): + self.client.get("/api/stats") + + +class ApiConsumer(FastHttpUser): + """A script/integration hitting the JSON API directly, no think time.""" + + weight = 1 + wait_time = between(0.1, 1) + + @task(3) + def latest(self): + self.client.get("/measurements/latest?limit=100") + + @task(3) + def hii_feeds(self): + self.client.get(random.choice( + ["/api/hii/waterlevel/latest", "/api/hii/rainfall/latest"] + ), name="/api/hii/[feed]/latest") + + @task(2) + def forecast(self): + self.client.get("/forecast") + + @task(2) + def heavy_history(self): + self.client.get( + "/measurements/history/P.1?hours=8760", + name="/measurements/history/P.1 [heavy]", + ) + + @task(1) + def health(self): + self.client.get("/health") diff --git a/src/ml/predict.py b/src/ml/predict.py index 9e6fe5e..07fd38b 100644 --- a/src/ml/predict.py +++ b/src/ml/predict.py @@ -164,21 +164,24 @@ def _model_forecast( predicted_max = max(float(reg.predict(feature_row)[0]), current_level) sigma_h = bundle["sigma"].get(horizon_h, HEURISTIC_SIGMA) + # Belt-and-braces: the classifier head OR the regression-sigmoid path, + # whichever is more alarmed. The 2026-08-11 backtest showed a trained + # classifier staying silent through the 2024 record flood while the + # regression head tracked it — alerting must never be worse than the + # regression fallback. warn_head = ( None if thresholds_stale else bundle["heads"].get(f"warn_{horizon_h}") ) + p_warning = _sigmoid_probability(predicted_max, warn_thr, sigma_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) + p_warning = max(p_warning, float(warn_head.predict_proba(feature_row)[0][1])) danger_head = ( None if thresholds_stale else bundle["heads"].get(f"danger_{horizon_h}") ) + p_danger = _sigmoid_probability(predicted_max, danger_thr, sigma_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_danger = max(p_danger, float(danger_head.predict_proba(feature_row)[0][1])) p_warning = _clip_probability(p_warning) p_danger = min(_clip_probability(p_danger), p_warning) diff --git a/src/web_api.py b/src/web_api.py index cb49ad5..028f360 100644 --- a/src/web_api.py +++ b/src/web_api.py @@ -50,6 +50,7 @@ HISTORY_TTL = 300 # 5 minutes FORECAST_CACHE: Dict[str, tuple] = {} FORECAST_CACHE_LOCK = Lock() +FORECAST_COMPUTE_LOCK = asyncio.Lock() # single-flight for expensive inference FORECAST_TTL = 900 # 15 minutes DB_STATS_CACHE: Dict[str, tuple] = {} @@ -359,8 +360,10 @@ async def get_health(): if not health_manager: raise HTTPException(status_code=503, detail="Health manager not initialized") - # Run health checks (populates state read by get_health_summary) - health_manager.run_all_checks() + # Run health checks (populates state read by get_health_summary). + # In a thread: DatabaseHealthCheck and APIHealthCheck do blocking I/O and + # would otherwise stall the event loop for every other request. + await asyncio.to_thread(health_manager.run_all_checks) summary = health_manager.get_health_summary() return HealthResponse(**summary) @@ -747,6 +750,19 @@ async def get_flood_forecasts(): from .ml.predict import get_latest_forecasts except ImportError as error: raise HTTPException(status_code=503, detail=f"Forecasting unavailable: {error}") + # Single-flight: inference takes seconds; without this, N concurrent cache + # misses ran N full inferences and starved the thread pool (load test: + # /forecast timeouts at 10 concurrent clients rippled into every endpoint). + async with FORECAST_COMPUTE_LOCK: + with FORECAST_CACHE_LOCK: + cached = FORECAST_CACHE.get("all") + if cached and time.monotonic() - cached[0] < FORECAST_TTL: + return cached[1] + return await _compute_forecasts(get_latest_forecasts) + + +async def _compute_forecasts(get_latest_forecasts): + now = time.monotonic() try: data = await asyncio.to_thread(get_latest_forecasts) except FileNotFoundError: @@ -771,7 +787,9 @@ async def get_latest_measurements(limit: int = 100): raise HTTPException(status_code=503, detail="Database not available") try: - measurements = scraper.get_latest_data(limit=limit) + # In a thread: this is a synchronous DB query, and this is the most + # frequently hit endpoint — inline it would block the event loop. + measurements = await asyncio.to_thread(scraper.get_latest_data, limit) return [_to_measurement_response(m) for m in measurements]