From 6d4beae65a23bfafbb05c5f0bda7e9be11b6c837 Mon Sep 17 00:00:00 2001 From: grabowski Date: Sun, 18 Oct 2026 14:07:56 +0700 Subject: [PATCH] fix(ml): serve the forecast rain the model was trained with build_features cut the rain series to the gauge grid before the forward sum, and the grid ends at the newest reading, so the row every live forecast is made from had rain_fc24 = NaN. Training and all evaluations had the real Open-Meteo forecast there. The forward sum now runs over the grid extended by 24 h; earlier rows are unchanged. Harness (new opt-in variant rise_rain_fc24_missing = rise_rain with rain_fc24 blanked on test rows, five folds, --from-cache): P.1 and P.103 nearly unaffected, P.21 false alarms 10 -> 5 and MAE 0.232 -> 0.215 with the fix. Write-up in docs/FLOOD_FORECASTING.md section 5; results in models/eval_2026-10-18_fc24.json. Also: HII mirror rows in one hour now resolve to the latest reading via a stable sort on the raw time (was an unstable sort on the floored hour, i.e. row-order dependent). Needed for an exact Rust port. --- docs/FLOOD_FORECASTING.md | 28 ++ models/eval_2026-10-18_fc24.json | 737 +++++++++++++++++++++++++++++++ src/ml/data.py | 15 +- src/ml/evaluate.py | 17 + src/ml/features.py | 21 +- tests/test_rain_serving.py | 64 +++ 6 files changed, 876 insertions(+), 6 deletions(-) create mode 100644 models/eval_2026-10-18_fc24.json create mode 100644 tests/test_rain_serving.py diff --git a/docs/FLOOD_FORECASTING.md b/docs/FLOOD_FORECASTING.md index 6de304d..5c38013 100644 --- a/docs/FLOOD_FORECASTING.md +++ b/docs/FLOOD_FORECASTING.md @@ -584,6 +584,34 @@ span. `src/ml/hii_rain.py` builds the catchment mean and a 24 h-sum bias/MAE/correlation, so the two sources' relationship is on record by the time the 2027 fold (train ≤ 2027-04-30, test Jun–Nov 2027) can test it. +### 2026-10-18: forecast rain was missing at serving time — fixed + +`build_features` cut the rain series to the gauge grid before taking the +forward sum, and the grid ends at the newest reading. So the row serving +predicts from always had `rain_fc24 = NaN`, although training and every +evaluation above had the real value (the Open-Meteo forecast, already fetched +with `forecast_days=2`). The model was answering as if no weather forecast +existed. The forward sum now runs over the grid extended by 24 h; earlier rows +are unchanged, so training is unaffected. + +Cost of the bug, measured with the harness: the `rise_rain_fc24_missing` +variant is `rise_rain` with `rain_fc24` blanked on the test rows (results in +`models/eval_2026-10-18_fc24.json`, five monsoon folds, `--from-cache`): + +| Station | Variant | Events | Warned early | Median lead | False alarms | MAE | +|---|---|---|---|---|---|---| +| P.1 | fixed (`rise_rain`) | 4 | 4 | 7.5 h | 0 | 0.087 | +| P.1 | as served | 4 | 4 | 7.0 h | 0 | 0.087 | +| P.103 | fixed | 7 | 7 | 9.0 h | 2 | 0.152 | +| P.103 | as served | 7 | 7 | 9.0 h | 2 | 0.153 | +| P.21 | fixed | 10 | 8 | 2.5 h | 5 | 0.215 | +| P.21 | as served | 10 | 7 | 2.0 h | 10 | 0.232 | + +P.1 and P.103 lost little; P.21 doubled its false alarms. The two very long +"as served" leads (P.21 2025-08-26 +69 h, P.103 2025-10-03 +35 h) are alarm +episodes that happened to precede an event, not real early warnings. Since the +fix, the served features match what this section evaluates. + ## 6. Deployment ### API diff --git a/models/eval_2026-10-18_fc24.json b/models/eval_2026-10-18_fc24.json new file mode 100644 index 0000000..e6f487d --- /dev/null +++ b/models/eval_2026-10-18_fc24.json @@ -0,0 +1,737 @@ +[ + { + "station": "P.1", + "warn_thr": 3.7, + "folds": [ + { + "year": 2021, + "n_train": 20024, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.07375926701460789, + "mae_above_2p5": null, + "brier_warn": 0.0, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.07349391292098781, + "mae_above_2p5": null, + "brier_warn": 0.0, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2022, + "n_train": 28784, + "n_test": 4392, + "events": [ + { + "crossing": "2022-10-02T19:00:00", + "peak_ts": "2022-10-03T15:00:00", + "peak_level": 4.65 + } + ], + "variants": { + "rise_rain": { + "mae": 0.08452847354970178, + "mae_above_2p5": 0.25082252888260664, + "brier_warn": 0.004300908725927739, + "events": [ + { + "crossing": "2022-10-02T19:00:00", + "lead_h": 5.0, + "peak_level": 4.65, + "peak_pred_24h_before": 3.8173954245046406 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.08511964468019045, + "mae_above_2p5": 0.2527703432441232, + "brier_warn": 0.004315341268867166, + "events": [ + { + "crossing": "2022-10-02T19:00:00", + "lead_h": 5.0, + "peak_level": 4.65, + "peak_pred_24h_before": 3.8102898827310137 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2023, + "n_train": 37539, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.07696637416933236, + "mae_above_2p5": 0.10643655855133666, + "brier_warn": 1.919860722404625e-15, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.07281928830426289, + "mae_above_2p5": 0.12255431002563086, + "brier_warn": 1.859416740909973e-15, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2024, + "n_train": 46323, + "n_test": 4392, + "events": [ + { + "crossing": "2024-09-24T17:00:00", + "peak_ts": "2024-09-26T02:00:00", + "peak_level": 4.93 + }, + { + "crossing": "2024-10-03T09:00:00", + "peak_ts": "2024-10-05T12:00:00", + "peak_level": 5.3 + } + ], + "variants": { + "rise_rain": { + "mae": 0.0890019157543212, + "mae_above_2p5": 0.2093245927883516, + "brier_warn": 0.005826169840715695, + "events": [ + { + "crossing": "2024-09-24T17:00:00", + "lead_h": 10.0, + "peak_level": 4.93, + "peak_pred_24h_before": 4.817519939833057 + }, + { + "crossing": "2024-10-03T09:00:00", + "lead_h": 21.0, + "peak_level": 5.3, + "peak_pred_24h_before": 5.546588884631041 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.0896559193295349, + "mae_above_2p5": 0.2184968282693981, + "brier_warn": 0.00611626529479383, + "events": [ + { + "crossing": "2024-09-24T17:00:00", + "lead_h": 9.0, + "peak_level": 4.93, + "peak_pred_24h_before": 4.7051072221287855 + }, + { + "crossing": "2024-10-03T09:00:00", + "lead_h": 21.0, + "peak_level": 5.3, + "peak_pred_24h_before": 5.563518687086929 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2025, + "n_train": 55083, + "n_test": 4392, + "events": [ + { + "crossing": "2025-09-27T18:00:00", + "peak_ts": "2025-09-27T22:00:00", + "peak_level": 3.93 + } + ], + "variants": { + "rise_rain": { + "mae": 0.11202835461695825, + "mae_above_2p5": 0.2456782496767171, + "brier_warn": 0.005178356039745929, + "events": [ + { + "crossing": "2025-09-27T18:00:00", + "lead_h": 2.0, + "peak_level": 3.93, + "peak_pred_24h_before": 3.23 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.11232147044143151, + "mae_above_2p5": 0.2457513740974484, + "brier_warn": 0.004988582804599792, + "events": [ + { + "crossing": "2025-09-27T18:00:00", + "lead_h": 4.0, + "peak_level": 3.93, + "peak_pred_24h_before": 3.3264027753671597 + } + ], + "false_alarm_episodes": 0 + } + } + } + ] + }, + { + "station": "P.103", + "warn_thr": 5.95, + "folds": [ + { + "year": 2021, + "n_train": 20009, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.14438683486790602, + "mae_above_2p5": null, + "brier_warn": 0.0, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.14435605468178428, + "mae_above_2p5": null, + "brier_warn": 0.0, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2022, + "n_train": 28769, + "n_test": 4392, + "events": [ + { + "crossing": "2022-08-14T04:00:00", + "peak_ts": "2022-08-14T08:00:00", + "peak_level": 6.09 + }, + { + "crossing": "2022-10-02T16:00:00", + "peak_ts": "2022-10-03T16:00:00", + "peak_level": 7.54 + } + ], + "variants": { + "rise_rain": { + "mae": 0.14392334606606189, + "mae_above_2p5": 0.41075382386412596, + "brier_warn": 0.00764679988213082, + "events": [ + { + "crossing": "2022-08-14T04:00:00", + "lead_h": 6.0, + "peak_level": 6.09, + "peak_pred_24h_before": 4.825270553204425 + }, + { + "crossing": "2022-10-02T16:00:00", + "lead_h": 9.0, + "peak_level": 7.54, + "peak_pred_24h_before": 7.194064117976157 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.1438805977972288, + "mae_above_2p5": 0.41655333933205874, + "brier_warn": 0.007670614570637225, + "events": [ + { + "crossing": "2022-08-14T04:00:00", + "lead_h": 5.0, + "peak_level": 6.09, + "peak_pred_24h_before": 4.822220911656332 + }, + { + "crossing": "2022-10-02T16:00:00", + "lead_h": 9.0, + "peak_level": 7.54, + "peak_pred_24h_before": 7.187888665288963 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2023, + "n_train": 37524, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.16195301512514512, + "mae_above_2p5": null, + "brier_warn": 3.0208441147560167e-07, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.15117839869986002, + "mae_above_2p5": null, + "brier_warn": 6.508469941612771e-08, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2024, + "n_train": 46308, + "n_test": 4058, + "events": [ + { + "crossing": "2024-09-24T10:00:00", + "peak_ts": "2024-09-26T00:00:00", + "peak_level": 8.27 + }, + { + "crossing": "2024-09-30T03:00:00", + "peak_ts": "2024-09-30T06:00:00", + "peak_level": 5.99 + }, + { + "crossing": "2024-10-03T06:00:00", + "peak_ts": "2024-10-05T07:00:00", + "peak_level": 9.93 + } + ], + "variants": { + "rise_rain": { + "mae": 0.13901842325128816, + "mae_above_2p5": 0.2887416310966684, + "brier_warn": 0.0072040857753137046, + "events": [ + { + "crossing": "2024-09-24T10:00:00", + "lead_h": 19.0, + "peak_level": 8.27, + "peak_pred_24h_before": 8.212734363860193 + }, + { + "crossing": "2024-09-30T03:00:00", + "lead_h": 9.0, + "peak_level": 5.99, + "peak_pred_24h_before": 5.47565489914091 + }, + { + "crossing": "2024-10-03T06:00:00", + "lead_h": 55.0, + "peak_level": 9.93, + "peak_pred_24h_before": 8.780278464915938 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.15101681306989473, + "mae_above_2p5": 0.3221006670360109, + "brier_warn": 0.008559364941781224, + "events": [ + { + "crossing": "2024-09-24T10:00:00", + "lead_h": 18.0, + "peak_level": 8.27, + "peak_pred_24h_before": 7.997091483400909 + }, + { + "crossing": "2024-09-30T03:00:00", + "lead_h": 9.0, + "peak_level": 5.99, + "peak_pred_24h_before": 5.502287543201456 + }, + { + "crossing": "2024-10-03T06:00:00", + "lead_h": 62.0, + "peak_level": 9.93, + "peak_pred_24h_before": 9.103449390810866 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2025, + "n_train": 54577, + "n_test": 4392, + "events": [ + { + "crossing": "2025-09-26T06:00:00", + "peak_ts": "2025-09-27T21:00:00", + "peak_level": 6.64 + }, + { + "crossing": "2025-10-03T06:00:00", + "peak_ts": "2025-10-03T12:00:00", + "peak_level": 6.14 + } + ], + "variants": { + "rise_rain": { + "mae": 0.17206685418009837, + "mae_above_2p5": 0.3840520085986099, + "brier_warn": 0.015824852891785323, + "events": [ + { + "crossing": "2025-09-26T06:00:00", + "lead_h": 6.0, + "peak_level": 6.64, + "peak_pred_24h_before": 5.73043665401637 + }, + { + "crossing": "2025-10-03T06:00:00", + "lead_h": 8.0, + "peak_level": 6.14, + "peak_pred_24h_before": 5.609919760117381 + } + ], + "false_alarm_episodes": 2 + }, + "rise_rain_fc24_missing": { + "mae": 0.17368174393168262, + "mae_above_2p5": 0.3756076666801506, + "brier_warn": 0.014794520234147523, + "events": [ + { + "crossing": "2025-09-26T06:00:00", + "lead_h": 5.0, + "peak_level": 6.64, + "peak_pred_24h_before": 5.8809952440025 + }, + { + "crossing": "2025-10-03T06:00:00", + "lead_h": 35.0, + "peak_level": 6.14, + "peak_pred_24h_before": 5.6215400644356235 + } + ], + "false_alarm_episodes": 2 + } + } + } + ] + }, + { + "station": "P.21", + "warn_thr": 3.2, + "folds": [ + { + "year": 2021, + "n_train": 19099, + "n_test": 4392, + "events": [], + "variants": { + "rise_rain": { + "mae": 0.31976019104736464, + "mae_above_2p5": 0.49280187683998045, + "brier_warn": 4.599376306869306e-06, + "events": [], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.3238133154815793, + "mae_above_2p5": 0.49187453544981, + "brier_warn": 4.609640062881583e-06, + "events": [], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2022, + "n_train": 27860, + "n_test": 4392, + "events": [ + { + "crossing": "2022-09-12T05:00:00", + "peak_ts": "2022-09-12T05:00:00", + "peak_level": 3.21 + }, + { + "crossing": "2022-10-02T04:00:00", + "peak_ts": "2022-10-03T21:00:00", + "peak_level": 3.92 + } + ], + "variants": { + "rise_rain": { + "mae": 0.188862265489962, + "mae_above_2p5": 0.22738565836625405, + "brier_warn": 0.009674093199871036, + "events": [ + { + "crossing": "2022-09-12T05:00:00", + "lead_h": 0.0, + "peak_level": 3.21, + "peak_pred_24h_before": 2.49615362708868 + }, + { + "crossing": "2022-10-02T04:00:00", + "lead_h": 1.0, + "peak_level": 3.92, + "peak_pred_24h_before": 3.907087039152045 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.19791868293297726, + "mae_above_2p5": 0.2739869896442161, + "brier_warn": 0.0103072884296978, + "events": [ + { + "crossing": "2022-09-12T05:00:00", + "lead_h": 0.0, + "peak_level": 3.21, + "peak_pred_24h_before": 2.4722245571508683 + }, + { + "crossing": "2022-10-02T04:00:00", + "lead_h": 0.0, + "peak_level": 3.92, + "peak_pred_24h_before": 3.845235629130954 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2023, + "n_train": 36614, + "n_test": 4392, + "events": [ + { + "crossing": "2023-10-10T06:00:00", + "peak_ts": "2023-10-10T06:00:00", + "peak_level": 3.2 + } + ], + "variants": { + "rise_rain": { + "mae": 0.17283257678933697, + "mae_above_2p5": 0.39237224972369483, + "brier_warn": 0.012908395601661503, + "events": [ + { + "crossing": "2023-10-10T06:00:00", + "lead_h": 15.0, + "peak_level": 3.2, + "peak_pred_24h_before": 1.8019095357375217 + } + ], + "false_alarm_episodes": 3 + }, + "rise_rain_fc24_missing": { + "mae": 0.1953027570197896, + "mae_above_2p5": 0.37345466270918226, + "brier_warn": 0.005336247430915456, + "events": [ + { + "crossing": "2023-10-10T06:00:00", + "lead_h": 14.0, + "peak_level": 3.2, + "peak_pred_24h_before": 1.8796882667592647 + } + ], + "false_alarm_episodes": 3 + } + } + }, + { + "year": 2024, + "n_train": 45398, + "n_test": 4385, + "events": [ + { + "crossing": "2024-09-01T03:00:00", + "peak_ts": "2024-09-01T08:00:00", + "peak_level": 4.3 + }, + { + "crossing": "2024-09-23T06:00:00", + "peak_ts": "2024-09-23T15:00:00", + "peak_level": 4.68 + }, + { + "crossing": "2024-10-03T05:00:00", + "peak_ts": "2024-10-04T20:00:00", + "peak_level": 5.21 + } + ], + "variants": { + "rise_rain": { + "mae": 0.20733600140331027, + "mae_above_2p5": 0.4173432645880738, + "brier_warn": 0.013730524899775963, + "events": [ + { + "crossing": "2024-09-01T03:00:00", + "lead_h": 0.0, + "peak_level": 4.3, + "peak_pred_24h_before": 1.7973215452596953 + }, + { + "crossing": "2024-09-23T06:00:00", + "lead_h": 2.0, + "peak_level": 4.68, + "peak_pred_24h_before": 2.670881094688805 + }, + { + "crossing": "2024-10-03T05:00:00", + "lead_h": 33.0, + "peak_level": 5.21, + "peak_pred_24h_before": 4.999142187693036 + } + ], + "false_alarm_episodes": 0 + }, + "rise_rain_fc24_missing": { + "mae": 0.2222046566812092, + "mae_above_2p5": 0.43423994130380844, + "brier_warn": 0.01611358822165686, + "events": [ + { + "crossing": "2024-09-01T03:00:00", + "lead_h": 0.0, + "peak_level": 4.3, + "peak_pred_24h_before": 1.7847018175539242 + }, + { + "crossing": "2024-09-23T06:00:00", + "lead_h": 1.0, + "peak_level": 4.68, + "peak_pred_24h_before": 2.640820102168781 + }, + { + "crossing": "2024-10-03T05:00:00", + "lead_h": 32.0, + "peak_level": 5.21, + "peak_pred_24h_before": 4.749354278189956 + } + ], + "false_alarm_episodes": 0 + } + } + }, + { + "year": 2025, + "n_train": 54150, + "n_test": 4392, + "events": [ + { + "crossing": "2025-07-23T15:00:00", + "peak_ts": "2025-07-24T04:00:00", + "peak_level": 3.63 + }, + { + "crossing": "2025-08-26T19:00:00", + "peak_ts": "2025-08-26T23:00:00", + "peak_level": 3.48 + }, + { + "crossing": "2025-09-12T21:00:00", + "peak_ts": "2025-09-12T22:00:00", + "peak_level": 3.21 + }, + { + "crossing": "2025-09-22T20:00:00", + "peak_ts": "2025-09-22T21:00:00", + "peak_level": 3.27 + } + ], + "variants": { + "rise_rain": { + "mae": 0.1851599197920176, + "mae_above_2p5": 0.29411810271343447, + "brier_warn": 0.021849546400589632, + "events": [ + { + "crossing": "2025-07-23T15:00:00", + "lead_h": 3.0, + "peak_level": 3.63, + "peak_pred_24h_before": 2.4024371626847 + }, + { + "crossing": "2025-08-26T19:00:00", + "lead_h": 4.0, + "peak_level": 3.48, + "peak_pred_24h_before": 2.492383184457359 + }, + { + "crossing": "2025-09-12T21:00:00", + "lead_h": 18.0, + "peak_level": 3.21, + "peak_pred_24h_before": 2.5690976144688706 + }, + { + "crossing": "2025-09-22T20:00:00", + "lead_h": 1.0, + "peak_level": 3.27, + "peak_pred_24h_before": 2.382265050508778 + } + ], + "false_alarm_episodes": 2 + }, + "rise_rain_fc24_missing": { + "mae": 0.21936157981899435, + "mae_above_2p5": 0.3121582241464527, + "brier_warn": 0.025980409031233827, + "events": [ + { + "crossing": "2025-07-23T15:00:00", + "lead_h": 2.0, + "peak_level": 3.63, + "peak_pred_24h_before": 2.096100137296138 + }, + { + "crossing": "2025-08-26T19:00:00", + "lead_h": 69.0, + "peak_level": 3.48, + "peak_pred_24h_before": 2.3678295230657067 + }, + { + "crossing": "2025-09-12T21:00:00", + "lead_h": 19.0, + "peak_level": 3.21, + "peak_pred_24h_before": 2.783241642058761 + }, + { + "crossing": "2025-09-22T20:00:00", + "lead_h": 2.0, + "peak_level": 3.27, + "peak_pred_24h_before": 2.279238629032059 + } + ], + "false_alarm_episodes": 7 + } + } + } + ] + } +] \ No newline at end of file diff --git a/src/ml/data.py b/src/ml/data.py index 906ba66..1ce3b77 100644 --- a/src/ml/data.py +++ b/src/ml/data.py @@ -146,15 +146,22 @@ def _fetch_hii_levels( with engine.connect() as connection: df = pd.read_sql(text(query), connection, params=params) + return _prepare_hii_levels(df) + + +def _prepare_hii_levels(df: pd.DataFrame) -> pd.DataFrame: + """Hourly mirror rows. Several readings in one hour: the latest wins + (stable sort on the raw time, so the choice never depends on row order).""" df = df.dropna(subset=["wl_msl"]) if df.empty: return df - df["timestamp"] = pd.to_datetime(df["timestamp"]).dt.floor("h") + df = df.copy() + df["timestamp"] = pd.to_datetime(df["timestamp"]) + df = df.sort_values("timestamp", kind="stable") + df["timestamp"] = df["timestamp"].dt.floor("h") df["wl_msl"] = pd.to_numeric(df["wl_msl"], errors="coerce") df["discharge"] = pd.to_numeric(df["discharge"], errors="coerce") - df = df.sort_values("timestamp").drop_duplicates( - subset=["station_code", "timestamp"], keep="last" - ) + df = df.drop_duplicates(subset=["station_code", "timestamp"], keep="last") return df diff --git a/src/ml/evaluate.py b/src/ml/evaluate.py index 0ae8e0c..bc5dfdb 100644 --- a/src/ml/evaluate.py +++ b/src/ml/evaluate.py @@ -72,6 +72,7 @@ class Variant: use_dam: bool = False, use_fc48: bool = False, qsigma: bool = False, + fc24_missing_at_serve: bool = False, ): self.name = name self.target = target # 'abs' or 'rise' @@ -84,6 +85,10 @@ class Variant: # behaviour of the deployed model exactly, since p>=0.5 alerts are # sigma-independent) and quantile heads ONLY for a per-row sigma. self.qsigma = qsigma + # Test rows get rain_fc24 = NaN: what production served until the + # 2026-10 fix (the rain series was cut to the gauge grid before the + # forward sum, so the newest row never had forecast rain). + self.fc24_missing_at_serve = fc24_missing_at_serve def fit_predict(self, X_tr, y_abs_tr, X_te) -> Tuple[np.ndarray, np.ndarray]: if not self.use_rain: @@ -108,6 +113,9 @@ class Variant: X_te = X_te.drop(columns=drop) elif "rain_fc48" not in X_tr.columns: raise ValueError(f"{self.name} requires the rain series") + if self.fc24_missing_at_serve and "rain_fc24" in X_te.columns: + X_te = X_te.copy() + X_te["rain_fc24"] = np.nan level_tr = X_tr["level"] level_te = X_te["level"].to_numpy() y_tr = (y_abs_tr - level_tr) if self.target == "rise" else y_abs_tr @@ -164,6 +172,14 @@ VARIANTS: Dict[str, Variant] = { "rise_rain_qsigma": Variant( "rise_rain_qsigma", target="rise", use_rain=True, qsigma=True ), + # 2026-10-18: the deployed rise_rain model as it was actually served + # before the rain_fc24 serving fix (forecast rain missing at test time). + "rise_rain_fc24_missing": Variant( + "rise_rain_fc24_missing", + target="rise", + use_rain=True, + fc24_missing_at_serve=True, + ), } # Dam variants are opt-in by name: they require dam columns that only exist @@ -176,6 +192,7 @@ DEFAULT_VARIANTS = [ if not v.use_dam and not v.use_fc48 and not v.qsigma + and not v.fc24_missing_at_serve and not (v.quantile and v.use_rain) ] diff --git a/src/ml/features.py b/src/ml/features.py index 3b9dd3b..603a300 100644 --- a/src/ml/features.py +++ b/src/ml/features.py @@ -209,6 +209,7 @@ def _hours_since_observed(mask_col: pd.Series) -> pd.Series: RAIN_FEATURES = ("rain_6h", "rain_24h", "rain_72h", "rain_fc24") +FORWARD_RAIN_H = 24 # rain_fc24 window DAM_FEATURES = ("dam_storage_pct", "dam_storage_pct_d3", "dam_inflow", "dam_outflow") # Stations hydrologically downstream of the Mae Ngat confluence (Ping mainstem @@ -297,9 +298,25 @@ def build_features( cols["rain_6h"] = r.rolling(6, min_periods=1).sum() cols["rain_24h"] = r.rolling(24, min_periods=1).sum() cols["rain_72h"] = r.rolling(72, min_periods=1).sum() - # forward sum over (t, t+24]: shift(-1) starts the window at t+1 + # forward sum over (t, t+24]: shift(-1) starts the window at t+1. + # Computed on the grid extended by the 24 h after its last row: the + # last row is the one serving predicts from, and cutting the rain + # series to the grid made its forecast rain always NaN (unlike every + # training row since 2021). Earlier rows keep their values. + future = ( + pd.date_range(idx[-1], periods=FORWARD_RAIN_H + 1, freq="h")[1:] + if len(idx) + else idx[:0] + ) + r_ext = rain.reindex(idx.append(future)) cols["rain_fc24"] = ( - r.shift(-1).iloc[::-1].rolling(24, min_periods=1).sum().iloc[::-1] + r_ext.shift(-1) + .iloc[::-1] + .rolling(FORWARD_RAIN_H, min_periods=1) + .sum() + .iloc[::-1] + .iloc[: len(idx)] + .set_axis(idx) ) if dam is not None and station in DAM_STATIONS: diff --git a/tests/test_rain_serving.py b/tests/test_rain_serving.py new file mode 100644 index 0000000..55c0a36 --- /dev/null +++ b/tests/test_rain_serving.py @@ -0,0 +1,64 @@ +"""rain_fc24 at the newest grid row (the one serving predicts from).""" + +import numpy as np +import pandas as pd + +from src.ml import features +from src.ml.data import _prepare_hii_levels + + +def _grid(hours=48): + idx = pd.date_range("2026-09-01", periods=hours, freq="h") + df = pd.DataFrame( + { + "timestamp": idx, + "station_code": "P.1", + "water_level": np.linspace(1.0, 2.0, hours), + "discharge": np.nan, + } + ) + return features.make_hourly_grid(df) + + +def test_forward_rain_reaches_past_the_last_reading(): + grid = _grid() + idx = grid.observed.index + rain = pd.Series( + np.arange(1, len(idx) + 49, dtype=float), + index=pd.date_range(idx[0], periods=len(idx) + 48, freq="h"), + ) + X = features.build_features(grid, "P.1", rain=rain) + last = idx[-1] + expected = rain.loc[last + pd.Timedelta(hours=1) : last + pd.Timedelta(hours=24)] + assert X.loc[last, "rain_fc24"] == expected.sum() + # rows whose window lies inside the grid are what they always were + old = rain.reindex(idx).shift(-1).iloc[::-1].rolling(24, min_periods=1).sum() + inner = idx[:-24] + np.testing.assert_allclose(X.loc[inner, "rain_fc24"], old.iloc[::-1].loc[inner]) + + +def test_forward_rain_nan_without_forecast_hours(): + grid = _grid() + idx = grid.observed.index + rain = pd.Series(1.0, index=idx) # nothing after the last reading + X = features.build_features(grid, "P.1", rain=rain) + assert np.isnan(X.loc[idx[-1], "rain_fc24"]) + assert X.loc[idx[-2], "rain_fc24"] == 1.0 + + +def test_hii_rows_latest_reading_in_an_hour_wins(): + raw = pd.DataFrame( + { + "timestamp": pd.to_datetime( + ["2026-10-01 10:40", "2026-10-01 10:10", "2026-10-01 11:00"] + ), + "station_code": "P.1", + "wl_msl": [301.5, 301.2, 301.7], + "discharge": [10.0, 9.0, 11.0], + } + ) + out = _prepare_hii_levels(raw) + assert list(out["wl_msl"]) == [301.5, 301.7] + assert list(out["timestamp"]) == list( + pd.to_datetime(["2026-10-01 10:00", "2026-10-01 11:00"]) + )