grabowski
7b31d4d0dd
CI / Test suite (push) Successful in 22s
CI / Format & lint (push) Successful in 16s
Docs / Validate documentation (push) Successful in 10s
Security / Dependency vulnerabilities (push) Successful in 1m34s
Security / Static analysis (push) Successful in 10s
Security / License report (push) Successful in 13s
src/ml/skill.py joins forecast_history (what each deployed version predicted for the 24 h peak, hourly) to water_measurements (what the river did) and reports per version: verified hours, peak MAE, bias, the persistence baseline (peak = current level), skill = 1 - MAE/persistence, and the same MAE restricted to observed peaks >= 2 m. Only forecasts whose window has elapsed with >= 75 % of hours observed count; a version needs 24 verified hours before it is compared. GET /api/forecast/skill?station_code=P.1&horizon=24 returns it (SWR cached, 15 min). The dashboard's forecast card gains a panel with a one-line verdict (current vs previous version), the per-version table, and a caveat that quiet weeks measure quiet-river accuracy only: the model is judged on flood-onset lead, which the backtests cover. EN + TH. On today's production data: hgb-v3+28b62e5 (369 h, Aug 13 - Sep 1) MAE 15.2 cm, skill -0.05; hgb-v2+f6570ac (224 h, Sep 1 - 11) MAE 12.3 cm, skill 0.36 - the "worse" v2 model scores better on a quieter fortnight, which is exactly why the panel shows the >= 2 m column and the caveat. Tests: 3, sqlite, synthetic. scripts/dev_proxy.py: DEV_PROXY_LOCAL lets a not-yet-deployed endpoint be answered from a local JSON file while everything else goes to prod.
94 lines
3.9 KiB
Python
94 lines
3.9 KiB
Python
"""Forecast skill verification: issued forecasts vs observed peaks (sqlite)."""
|
|
|
|
import datetime
|
|
|
|
import pytest
|
|
from sqlalchemy import create_engine, text
|
|
|
|
from src.ml import skill
|
|
|
|
|
|
@pytest.fixture
|
|
def engine(tmp_path):
|
|
eng = create_engine(f"sqlite:///{tmp_path / 'skill.db'}")
|
|
with eng.begin() as c:
|
|
c.execute(text("CREATE TABLE stations (id INTEGER PRIMARY KEY, station_code TEXT)"))
|
|
c.execute(text("INSERT INTO stations VALUES (1, 'P.1')"))
|
|
c.execute(
|
|
text(
|
|
"CREATE TABLE water_measurements (timestamp DATETIME, station_id INTEGER, water_level REAL)"
|
|
)
|
|
)
|
|
c.execute(
|
|
text(
|
|
"CREATE TABLE forecast_history (as_of TIMESTAMP, station_code TEXT, horizon_hours INTEGER, "
|
|
"predicted_max_level REAL, p_warning REAL, p_danger REAL, current_level REAL, "
|
|
"model_version TEXT, source TEXT)"
|
|
)
|
|
)
|
|
return eng
|
|
|
|
|
|
def _fill(engine, start, hours, level_fn, forecasts):
|
|
"""hours of hourly observations from `start`, plus (as_of_offset_h, version, pred) rows."""
|
|
with engine.begin() as c:
|
|
for h in range(hours):
|
|
ts = start + datetime.timedelta(hours=h)
|
|
c.execute(
|
|
text("INSERT INTO water_measurements VALUES (:t, 1, :l)"),
|
|
{"t": ts, "l": level_fn(h)},
|
|
)
|
|
for off, version, pred in forecasts:
|
|
ts = start + datetime.timedelta(hours=off)
|
|
c.execute(
|
|
text(
|
|
"INSERT INTO forecast_history VALUES (:t, 'P.1', 24, :p, 0, 0, :cur, :v, 'model')"
|
|
),
|
|
{"t": ts, "p": pred, "cur": level_fn(off), "v": version},
|
|
)
|
|
|
|
|
|
def test_skill_per_version_and_trend(engine):
|
|
start = datetime.datetime(2026, 8, 1)
|
|
# river: flat 1.5 m, with a bump to 2.4 m around hour 100
|
|
level = lambda h: 2.4 if 96 <= h <= 104 else 1.5
|
|
forecasts = []
|
|
# old version: always predicts 1.5 (persistence-like, misses the bump)
|
|
for off in range(0, 60):
|
|
forecasts.append((off, "hgb-v2+aaaaaaa", 1.5))
|
|
# new version: predicts 1.5 normally and 2.3 ahead of the bump
|
|
for off in range(60, 200):
|
|
pred = 2.3 if 72 <= off <= 104 else 1.5
|
|
forecasts.append((off, "hgb-v3+bbbbbbb", pred))
|
|
_fill(engine, start, 260, level, forecasts)
|
|
|
|
out = skill.compute_skill(engine, "sqlite", "P.1", 24, now=start + datetime.timedelta(hours=300))
|
|
assert [v["model_version"] for v in out["versions"]] == ["hgb-v2+aaaaaaa", "hgb-v3+bbbbbbb"]
|
|
old, new = out["versions"]
|
|
assert old["n"] == 60 and old["enough_data"]
|
|
assert new["n"] == 140 and new["enough_data"]
|
|
# the old version issued only on flat hours: perfect there, no bump rows
|
|
assert old["mae_m"] == 0.0 and old["above_2m_n"] == 0
|
|
# the new version saw the bump: nonzero MAE but positive skill vs persistence
|
|
assert new["above_2m_n"] > 0
|
|
assert new["skill"] is not None and new["skill"] > 0
|
|
assert out["current"]["model_version"] == "hgb-v3+bbbbbbb"
|
|
assert out["trend"]["previous_version"] == "hgb-v2+aaaaaaa"
|
|
assert out["trend"]["better"] is False # honest: old had an easier period
|
|
|
|
|
|
def test_skill_requires_full_window(engine):
|
|
start = datetime.datetime(2026, 8, 1)
|
|
# forecasts issued at the very end have no observed window yet
|
|
_fill(engine, start, 30, lambda h: 1.5, [(o, "hgb-v3+ccccccc", 1.5) for o in range(0, 30)])
|
|
out = skill.compute_skill(engine, "sqlite", "P.1", 24, now=start + datetime.timedelta(hours=30))
|
|
# only as_of <= now-24h AND with >= 18 observed hours in the window count
|
|
assert out["versions"] and out["versions"][0]["n"] == 7 # as_of 0..6 h: <= now-24h with >= 18 observed hours
|
|
assert out["versions"][0]["enough_data"] is False
|
|
assert out["trend"] is None
|
|
|
|
|
|
def test_skill_empty(engine):
|
|
out = skill.compute_skill(engine, "sqlite", "P.1", 24)
|
|
assert out["versions"] == [] and out["current"] is None and out["trend"] is None
|