feat: codified backtests, honest docs, belt-and-braces serving, perf fixes

Retrained on the gap-filled DB (592k -> 976k rows) and re-examined the flood backtests, now reproducible via scripts/backtest_render.py (renders the three docs/img charts and gates on a >=12h 2024 first-alert lead — currently failing by design and documented as such). Findings, all documented in FLOOD_FORECASTING.md: the true 2024 crossing was 24 Sep 17:00 (8h earlier than recorded; confirmed against the independent HII sensor), the historical 24h-warning claim was partly a missing-data artifact, and retrained warn classifiers collapse on the filled grid (P.1 24h PR-AUC 0.900 -> 0.288) while regression MAE improves (11.3 -> 10.5 cm). Serving therefore becomes max(classifier, sigmoid(regression)) so alerting is never worse than the regression path; metrics table, head-gating tiers, honest-limits and runbook expectations all updated to the current model (hgb-v1+d2d0e65). Perf, from Locust load testing (scripts/locustfile.py + load_test.py): single-flight lock around /forecast inference (concurrent cache misses previously each ran ~18s inference and starved the shared thread pool; 200-user run after: 105 rps, 0.01% errors), and /measurements/latest + /health moved off the event loop (synchronous DB/network calls in async handlers were stalling every request under load).
2026-08-12 10:46:00 +07:00
parent d2d0e655aa
commit 0005f7dce1
9 changed files with 647 additions and 75 deletions
@@ -203,16 +203,17 @@ section 6).
 The system degrades in tiers rather than failing:
-1. **Classifier head**, when the training span contains at least
+1. **Belt-and-braces probability** *(since 2026-08-11 — see the re-examination
-   `MIN_POSITIVES_FOR_CLASSIFIER = 30` positive examples. Below that, a
+   note in section 7)*: the sigmoid-of-regression probability
-   classifier would be fitting noise, and the head is recorded in
+   `p = 1/(1 + exp(−(predicted_max − threshold)/σ))` is always computed (σ =
-   `skipped_heads` with its reason.
+   the regressor's test-residual std, floor `MIN_SIGMA = 0.15` m), and when a
-2. **Sigmoid on the regression head**, when the classifier is absent.
+   classifier head exists — trained only if the span had at least
-   `p = 1/(1 + exp(−(predicted_max − threshold)/σ))`, where σ is the standard
+   `MIN_POSITIVES_FOR_CLASSIFIER = 30` positives — the served probability is
-   deviation of the regressor's test residuals (floor `MIN_SIGMA = 0.15` m). This
+   `max(classifier, sigmoid)`. The classifier can raise the alarm but never
-   turns the peak-level prediction into a calibrated-ish probability that widens
+   silence it: on the gap-filled data a trained classifier stayed near zero
-   correctly when the regressor is less accurate at that horizon — at P.1, σ is
+   through the 2024 record crossing while the regression tracked it.
-   0.15 m at 6 and 12 h but 0.166 m at 24 h.
+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
+Two evaluations exist and they differ sharply — the re-examination note in the
-2025-01-01 → 2026-08-10 — the test span is entirely unseen future data relative
+next section explains why (the hourly grid was gap-filled from ~56% to ~93%
-to training.
+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 |
+| Horizon | Warning PR-AUC | 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% |
+| 6 h | 0.783 | 5.5 cm | 7.9 cm | 14,034 | 0.12% |
-| 12 h | 0.904 | 93.8% | 97.7% | 9.0 cm | 15.0 cm | 7,932 | 1.61% |
+| 12 h | 0.508 | 8.1 cm | 14.6 cm | 14,028 | 0.16% |
-| 24 h | 0.900 | 90.1% | 93.4% | 11.3 cm | 24.5 cm | 8,572 | 1.77% |
+| 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
+Level accuracy improved; standalone classifier discrimination did not survive
-margin over chance. "Recall at 1% false-alarm rate" is the operationally honest
+the data change (which is why serving is now `max(classifier, sigmoid)` — see
-number: at a threshold that fires on 1% of quiet hours, the 6 h model still
+"Head gating"). Recall-at-FAR is null at all horizons on this run. Across
-catches 98.3% of warning exceedances.
+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
+**Historical evaluation** (`hgb-v1+49a3de0`, 2026-08-10, pre-gap-fill DB —
-evaluate a danger head on the 2025–26 span (base rate 5.7–7.4%): PR-AUC 0.979 /
+kept for the record; these numbers described the sparser 56%-filled grid and do
-0.953 / 0.892 and recall at 1% FAR of 97.9% / 89.9% / 79.5% at 6 / 12 / 24 h.
+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
+| Horizon | Warning PR-AUC | Recall @1% FAR | MAE | Test rows | Base rate |
-(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,
+| 6 h | 0.974 | 98.3% | 6.1 cm | 8,536 | 1.36% |
-P.85 0.953, P.75 0.927 — and the weak ones are un-routed western tributaries with
+| 12 h | 0.904 | 93.8% | 9.0 cm | 7,932 | 1.61% |
-almost no positives (P.84 0.570, P.82 0.302 on 0.22% of test hours). P.92 and P.20
+| 24 h | 0.900 | 90.1% | 11.3 cm | 8,572 | 1.77% |
-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.
+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
+This uses the deployed configuration (train ≤ 2024-12-31) on an event it never
-deployed configuration (train ≤ 2024-12-31) on an event it never saw: the
+saw. *(Chart regenerated 2026-08-11 on the gap-filled data — see the
-model's **first alert came 26 September 2025 at 18:00, exactly 24 hours before
+re-examination note above; the original one-off render, trained on the sparser
-the river crossed 3.70 m**, and it predicted a 4.00 m peak against an actual
+data, alerted 24 h ahead and predicted the 3.93 m peak within 7 cm.)* On
-3.93 m — within 7 cm. Note the discrimination: the near-miss crest of 3.51 m on
+today's fuller dataset the retrained equivalent first alerts at **18:00 on
-26 September never triggered an alert, the probability fires only for the real
+27 September 2025 — as the river crosses 3.70 m**, not a day ahead. The
-event, and it stands down as the water recedes.
+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
+The render above shows the whole event hour by hour *(regenerated 2026-08-11
-(blue) against the 24 h-ahead predicted peak the model issued at each hour
+on the gap-filled data)*. Top: the observed level (blue) against the 24 h-ahead
-(amber, dashed) — the amber line leads the blue one into both flood waves,
+predicted peak the model issued at each hour (amber, dashed) — the amber line
-which is the entire point of the system. Bottom: the model's probability of
+tracks both flood waves but no longer clearly leads the first one. Bottom: the
-flooding within 24 h; it fires its **first alert at 01:00 on 24 September, a
+belt-and-braces probability of flooding within 24 h; on the fuller data its
-full day before the river crossed the 3.70 m flooding line**, stays pinned near
+**first alert comes at 11:00 on 25 September, ~18 hours after the true 17:00
-1.0 through both waves, and stands down between and after them. Also visible,
+24 September crossing**, then stays correctly alarmed through the October
-honestly: the predicted peak tops out ~0.4 m short of the actual 5.30 m record
+record wave. Also visible, honestly: the predicted peak tops out well short of
-(the extreme-compression limitation discussed below), and the prediction is
+the actual 5.30 m record (the extreme-compression limitation discussed below).
-noisier on the recession limbs. The same model track drives the dashboard's
+The same historic model track drives the dashboard's "Replay Oct 2024 flood"
-"Replay Oct 2024 flood" feature, so this chart can be watched live on the map.
+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
+The hour-by-hour detail of the detection window *(regenerated 2026-08-11)*
-predicted 24 h peak (amber) starts pulling away from the observed level late on
+shows the corrected sequence: the river crosses 3.70 m at **17:00 on
-23 September as upstream gauges rise, the warning probability snaps from ~0 to
+24 September** (the hours recovered by gap-filling; independently confirmed by
-1.0 at **01:00 on 24 September**, and the river crosses 3.70 m at **01:00 on
+the HII sensor at the same bridge), while the retrained model's probability
-25 September** — a clean 24-hour warning, delivered while the river in town
+only crosses 0.5 at **11:00 on 25 September**. The original render — sparser
-still looked normal at 2.8 m.
+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
+is extrapolating basin state and season, not routing a wave. The regenerated
-events bear this out: the 25 September 2025 cold-start crossing was called 7 h
+backtests bear this out — harder than first documented (see the re-examination
-ahead by the 12 h model and 9 h ahead by the 24 h model. **Practical lead for P.1
+note in section 7): on the gap-filled data the retrained configuration alerts
-is ~7–17 h.** Extending it requires rainfall forecasts and Mae Ngat/Mae Kuang dam
+the 27 September 2025 crossing *as it happens* and the 24 September 2024
-release data, neither of which this system currently ingests.
+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
+`failed` names its reason in the same payload. Compare against the *previous
-dropped materially below ~0.97 or its MAE has risen well above ~6 cm, investigate
+run's* `metrics.json`, not an absolute bar: after the 2026-08-11 gap-fill the
-before deploying — that usually means a data problem (a gauge that went quiet, or
+expected baseline is P.1 6 h warning PR-AUC ≈ 0.78 and MAE ≈ 5.5 cm (the
-a bad backfill) rather than a modelling one.
+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):
@@ -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())
@@ -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())
@@ -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")
@@ -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])
+            p_warning = max(p_warning, float(warn_head.predict_proba(feature_row)[0][1]))
        else:
            p_warning = _sigmoid_probability(predicted_max, warn_thr, sigma_h)
        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])
+            p_danger = max(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)
@@ -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)
+    # Run health checks (populates state read by get_health_summary).
-    health_manager.run_all_checks()
+    # 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]