fix: backtest/review findings in the flood-ML package

From the adversarial review and threshold backtest (swarm verification):

- predict.py: when a bundle's trained thresholds differ from the current
  config (deploy before retrain), skip its stale classifier heads and
  derive p_warning/p_danger from the regression + sigma against the
  CURRENT thresholds - the dashboard can no longer show contradictory
  old-threshold classifier output next to new-threshold stages
- features.py: decouple the low-coverage regression-label rescue from
  the warning threshold (now the station's own p97.5 level); the old
  coupling silently dropped 34% of P.5's regression training rows and
  cost +46% MAE when its threshold rose
- features.py: P.82 danger 3.80 -> 3.75 (3.80 was above the station's
  8-year maximum of 3.78, so danger could never train or fire)
- data.py / predict.py: anchor models/cache paths to the repo root; the
  relative paths silently returned zero rows when run from another CWD
- annotate P.4A thresholds as low-confidence (11 supporting readings)

47 tests pass. Retrain required for the label-rescue and P.82 changes
to reach the classifier heads.

src/ml/data.py

@@ -24,7 +24,9 @@ from .features import UPSTREAM_LEADS
 logger = logging.getLogger(__name__)
 
 DEFAULT_API_URL = "http://100.81.167.42:8000"
-CACHE_DIR = Path("models/cache")
+# Anchored to the repo root so training/prediction work from any CWD; a relative
+# path here silently produced 0 rows when the CLI ran outside the repo root.
+CACHE_DIR = Path(__file__).resolve().parents[2] / "models" / "cache"
 _MEASUREMENT_COLUMNS = ["timestamp", "station_code", "water_level", "discharge"]
 
 
@@ -95,7 +97,9 @@ def _fetch_from_db(
     return _normalize_long(df)
 
 
-def _fetch_station_from_api(api_url: str, station_code: str, hours: int, limit: int = 100000) -> pd.DataFrame:
+def _fetch_station_from_api(
+    api_url: str, station_code: str, hours: int, limit: int = 100000
+) -> pd.DataFrame:
     import requests
 
     response = requests.get(
@@ -133,7 +137,9 @@ def _fetch_from_api(
     return _normalize_long(df)
 
 
-def _write_cache(df: pd.DataFrame, cache_dir: Path, source: str, discharge_maybe_synthetic: bool) -> None:
+def _write_cache(
+    df: pd.DataFrame, cache_dir: Path, source: str, discharge_maybe_synthetic: bool
+) -> None:
     cache_dir.mkdir(parents=True, exist_ok=True)
     for code, group in df.groupby("station_code"):
         path = cache_dir / f"{code}.csv.gz"
@@ -184,17 +190,23 @@ def load_measurements(
         try:
             df = _fetch_from_db(resolved_db_url, stations, start, end)
             if use_cache:
-                _write_cache(df, cache_dir, source="postgres", discharge_maybe_synthetic=False)
+                _write_cache(
+                    df, cache_dir, source="postgres", discharge_maybe_synthetic=False
+                )
             return df
         except Exception as error:
-            logger.warning(f"PostgreSQL fetch failed, falling back to HTTP API: {error}")
+            logger.warning(
+                f"PostgreSQL fetch failed, falling back to HTTP API: {error}"
+            )
 
     try:
         api_stations = stations or _default_stations()
         df = _fetch_from_api(api_url, api_stations, start, end)
         if not df.empty:
             if use_cache:
-                _write_cache(df, cache_dir, source="api", discharge_maybe_synthetic=True)
+                _write_cache(
+                    df, cache_dir, source="api", discharge_maybe_synthetic=True
+                )
             return df
     except Exception as error:
         logger.warning(f"HTTP API fetch failed: {error}")

src/ml/features.py

@@ -18,9 +18,8 @@ logger = logging.getLogger(__name__)
 # Static configuration
 # ---------------------------------------------------------------------------
 
-# Per-station (warning, danger) level thresholds in meters. "*" is the default
-# applied to any station without an explicit override.
 # Per-station (warning, danger) levels in metres on each gauge's own datum.
+# "*" is the default applied to any station without an explicit override.
 # Calibrated 2026-08-10 from the DB's discharge_percent (RID % of channel
 # capacity): warning = median level at 75-85% capacity, danger = median level
 # at 95-105%. P.1 instead uses the official Chiang Mai inundation map keyed to
@@ -32,6 +31,9 @@ THRESHOLDS: Dict[str, Tuple[float, float]] = {
     "P.103": (5.95, 6.75),
     "P.20": (2.35, 2.80),
     "P.21": (3.20, 3.60),
+    # P.4A: LOW CONFIDENCE — its sensor was dead 2019-2024 and only 11 readings
+    # ever reached 3.40 m, so the capacity calibration rests on very few points.
+    # It only affects the heuristic sigmoid (P.4A is NOT_TRAINABLE).
     "P.4A": (3.40, 3.90),
     "P.5": (4.55, 4.95),
     "P.67": (2.45, 2.90),
@@ -39,7 +41,9 @@ THRESHOLDS: Dict[str, Tuple[float, float]] = {
     "P.76": (5.35, 5.45),
     "P.77": (2.85, 3.35),
     "P.81": (5.15, 6.30),
-    "P.82": (3.40, 3.80),
+    # P.82 never reached 100% capacity in the record (max level 3.78, max 96.4%);
+    # danger sits just below the observed maximum so the head can actually train.
+    "P.82": (3.40, 3.75),
     "P.84": (3.45, 3.90),
     "P.85": (2.90, 3.35),
     "P.87": (3.75, 4.05),
@@ -73,8 +77,23 @@ BASIN_ANCHOR = "P.1"
 # lists, for each station, the (upstream_code, lead_hours) pairs to use as
 # routed-upstream input features when forecasting `station`.
 UPSTREAM_LEADS: Dict[str, List[Tuple[str, int]]] = {
-    "P.1": [("P.103", 1), ("P.67", 7), ("P.21", 9), ("P.75", 12), ("P.4A", 12), ("P.92", 15), ("P.20", 17)],
+    "P.1": [
-    "P.103": [("P.67", 6), ("P.21", 8), ("P.75", 11), ("P.4A", 11), ("P.92", 14), ("P.20", 16)],
+        ("P.103", 1),
+        ("P.67", 7),
+        ("P.21", 9),
+        ("P.75", 12),
+        ("P.4A", 12),
+        ("P.92", 15),
+        ("P.20", 17),
+    ],
+    "P.103": [
+        ("P.67", 6),
+        ("P.21", 8),
+        ("P.75", 11),
+        ("P.4A", 11),
+        ("P.92", 14),
+        ("P.20", 16),
+    ],
     "P.21": [("P.67", 1), ("P.75", 3), ("P.4A", 3), ("P.92", 6), ("P.20", 8)],
     "P.67": [("P.75", 5), ("P.4A", 5), ("P.92", 8), ("P.20", 10)],
     "P.75": [("P.92", 3), ("P.20", 5)],
@@ -140,9 +159,13 @@ def make_hourly_grid(df_long: pd.DataFrame) -> HourlyGrid:
     df["timestamp"] = pd.to_datetime(df["timestamp"]).dt.floor("h")
     df = df.drop_duplicates(subset=["station_code", "timestamp"], keep="last")
 
-    full_index = pd.date_range(df["timestamp"].min(), df["timestamp"].max(), freq="h", name="timestamp")
+    full_index = pd.date_range(
+        df["timestamp"].min(), df["timestamp"].max(), freq="h", name="timestamp"
+    )
 
-    wide = df.pivot(index="timestamp", columns="station_code", values=["water_level", "discharge"])
+    wide = df.pivot(
+        index="timestamp", columns="station_code", values=["water_level", "discharge"]
+    )
     wide = wide.reorder_levels([1, 0], axis=1).sort_index(axis=1)
     wide = wide.reindex(full_index)
 
@@ -154,7 +177,9 @@ def make_hourly_grid(df_long: pd.DataFrame) -> HourlyGrid:
     return HourlyGrid(observed=observed, filled=filled, mask=mask)
 
 
-def _series(grid_frame: pd.DataFrame, station: str, field: str, index: pd.Index) -> pd.Series:
+def _series(
+    grid_frame: pd.DataFrame, station: str, field: str, index: pd.Index
+) -> pd.Series:
     """Fetch a (station, field) column, or an all-NaN series if the station is absent."""
     if (station, field) in grid_frame.columns:
         return grid_frame[(station, field)]
@@ -214,7 +239,9 @@ def build_features(grid: HourlyGrid, station: str) -> pd.DataFrame:
         cols[f"{upstream_code}_level_lag_{lead_h}"] = u_level.shift(lead_h)
         cols[f"{upstream_code}_level_lag_{lead_h + 3}"] = u_level.shift(lead_h + 3)
         cols[f"{upstream_code}_rise_6_lag_{lead_h}"] = u_rise_6.shift(lead_h)
-        cols[f"{upstream_code}_rollmax_24_lag_{near_lag}"] = u_rollmax_24.shift(near_lag)
+        cols[f"{upstream_code}_rollmax_24_lag_{near_lag}"] = u_rollmax_24.shift(
+            near_lag
+        )
 
     if station != BASIN_ANCHOR:
         p1_level = _series(grid.filled, BASIN_ANCHOR, "water_level", idx)
@@ -244,11 +271,23 @@ def _future_window_stats(col: pd.Series, horizon_h: int) -> Tuple[pd.Series, pd.
     return fut_max, fut_count
 
 
-def build_labels(grid: HourlyGrid, station: str, horizons: Tuple[int, ...] = (6, 12, 24)) -> pd.DataFrame:
+def build_labels(
+    grid: HourlyGrid, station: str, horizons: Tuple[int, ...] = (6, 12, 24)
+) -> pd.DataFrame:
     """Build max-level and threshold-exceedance labels for one target station."""
     idx = grid.observed.index
     observed_level = _series(grid.observed, station, "water_level", idx)
     warn_thr, danger_thr = get_thresholds(station)
+    # Low-coverage windows are still usable regression labels when they contain a
+    # rare high reading. Anchor "rare" to the station's own distribution (p97.5),
+    # NOT to warn_thr: coupling it to the configurable threshold made raising a
+    # station's threshold silently shrink its regression training set (P.5 lost
+    # 34% of rows and +46% MAE when its warning went 3.0 -> 4.55).
+    rescue_thr = (
+        float(observed_level.quantile(0.975))
+        if observed_level.notna().any()
+        else np.inf
+    )
 
     out: Dict[str, pd.Series] = {}
     for horizon_h in horizons:
@@ -264,7 +303,7 @@ def build_labels(grid: HourlyGrid, station: str, horizons: Tuple[int, ...] = (6,
         exceed_danger[fut_max >= danger_thr] = 1.0
         exceed_danger[enough_cov & exceed_danger.isna()] = 0.0
 
-        max_level_valid = fut_max.where(enough_cov | (fut_max >= warn_thr))
+        max_level_valid = fut_max.where(enough_cov | (fut_max >= rescue_thr))
 
         out[f"max_level_{horizon_h}"] = max_level_valid
         out[f"exceed_warn_{horizon_h}"] = exceed_warn
@@ -297,13 +336,17 @@ def build_matrix(
     Y = Y.loc[keep]
 
     positive_counts = {
-        col: int(Y[col].sum()) for col in Y.columns if col.startswith("exceed_") and Y[col].notna().any()
+        col: int(Y[col].sum())
+        for col in Y.columns
+        if col.startswith("exceed_") and Y[col].notna().any()
     }
     meta = {
         "station_code": station,
         "n_rows": int(len(X)),
         "span": (
-            (X.index.min().isoformat(), X.index.max().isoformat()) if len(X) else (None, None)
+            (X.index.min().isoformat(), X.index.max().isoformat())
+            if len(X)
+            else (None, None)
         ),
         "positive_counts": positive_counts,
     }

src/ml/predict.py

@@ -20,6 +20,9 @@ from . import features
 
 logger = logging.getLogger(__name__)
 
+# Anchored to the repo root so the API finds trained bundles regardless of CWD.
+DEFAULT_MODELS_DIR = Path(__file__).resolve().parents[2] / "models"
+
 DEFAULT_HORIZONS: Tuple[int, ...] = (6, 12, 24)
 STALE_AFTER_H = 6.0
 HEURISTIC_SIGMA = 0.3
@@ -60,7 +63,9 @@ def _readings_to_long_df(readings_by_station: Dict[str, List[dict]]) -> pd.DataF
                 }
             )
     if not rows:
-        return pd.DataFrame(columns=["timestamp", "station_code", "water_level", "discharge"])
+        return pd.DataFrame(
+            columns=["timestamp", "station_code", "water_level", "discharge"]
+        )
     df = pd.DataFrame(rows)
     return df.dropna(subset=["timestamp"])
 
@@ -90,8 +95,12 @@ def _heuristic_forecast(
     results = []
     for horizon_h in horizons:
         predicted_max = max(current_level + rate * horizon_h * 0.7, current_level)
-        p_warning = _clip_probability(_sigmoid_probability(predicted_max, warn_thr, HEURISTIC_SIGMA))
+        p_warning = _clip_probability(
-        p_danger = _clip_probability(_sigmoid_probability(predicted_max, danger_thr, HEURISTIC_SIGMA))
+            _sigmoid_probability(predicted_max, warn_thr, HEURISTIC_SIGMA)
+        )
+        p_danger = _clip_probability(
+            _sigmoid_probability(predicted_max, danger_thr, HEURISTIC_SIGMA)
+        )
         p_danger = min(p_danger, p_warning)
         results.append(
             {
@@ -121,12 +130,28 @@ def _model_forecast(
 ) -> List[dict]:
     warn_thr = bundle["thresholds"]["warning"]
     danger_thr = bundle["thresholds"]["danger"]
+    # If thresholds changed since this bundle was trained, its classifier heads
+    # answer the OLD question (labels for the old levels) while stages/config use
+    # the new ones — a silent contradiction on the dashboard. Until a retrain,
+    # answer the current question consistently: use the regression + sigma against
+    # the configured thresholds and skip the stale heads.
+    cfg_warn, cfg_danger = features.get_thresholds(station_code)
+    thresholds_stale = (cfg_warn, cfg_danger) != (warn_thr, danger_thr)
+    if thresholds_stale:
+        logger.warning(
+            f"{station_code}: bundle thresholds ({warn_thr}, {danger_thr}) differ from "
+            f"configured ({cfg_warn}, {cfg_danger}); using regression-derived probabilities "
+            "until the model is retrained"
+        )
+        warn_thr, danger_thr = cfg_warn, cfg_danger
 
     feature_row = features.build_features(grid, station_code).loc[[as_of]]
     expected_columns = bundle["feature_names"]
     missing = [c for c in expected_columns if c not in feature_row.columns]
     if missing:
-        logger.error(f"Feature mismatch for {station_code} (missing {missing}); falling back to heuristic")
+        logger.error(
+            f"Feature mismatch for {station_code} (missing {missing}); falling back to heuristic"
+        )
         return None
     feature_row = feature_row[expected_columns]
 
@@ -139,13 +164,17 @@ def _model_forecast(
         predicted_max = max(float(reg.predict(feature_row)[0]), current_level)
         sigma_h = bundle["sigma"].get(horizon_h, HEURISTIC_SIGMA)
 
-        warn_head = bundle["heads"].get(f"warn_{horizon_h}")
+        warn_head = (
+            None if thresholds_stale else bundle["heads"].get(f"warn_{horizon_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)
 
-        danger_head = bundle["heads"].get(f"danger_{horizon_h}")
+        danger_head = (
+            None if thresholds_stale else bundle["heads"].get(f"danger_{horizon_h}")
+        )
         if danger_head is not None:
             p_danger = float(danger_head.predict_proba(feature_row)[0][1])
         else:
@@ -209,19 +238,35 @@ def _forecast_station(
 
     warn_thr, danger_thr = features.get_thresholds(station_code)
     t_minus_3 = as_of - pd.Timedelta(hours=3)
-    level_t_minus_3 = float(observed_level.loc[t_minus_3]) if t_minus_3 in observed_level.index else None
+    level_t_minus_3 = (
+        float(observed_level.loc[t_minus_3])
+        if t_minus_3 in observed_level.index
+        else None
+    )
 
     bundle_path = models_dir / f"flood_{station_code}.joblib"
     if not bundle_path.exists() or staleness_h > STALE_AFTER_H:
         return _heuristic_forecast(
-            station_code, as_of, current_level, level_t_minus_3, warn_thr, danger_thr, horizons
+            station_code,
+            as_of,
+            current_level,
+            level_t_minus_3,
+            warn_thr,
+            danger_thr,
+            horizons,
         )
 
     bundle = _load_bundle(bundle_path)
     model_results = _model_forecast(station_code, grid, bundle, as_of, current_level)
     if model_results is None:
         return _heuristic_forecast(
-            station_code, as_of, current_level, level_t_minus_3, warn_thr, danger_thr, horizons
+            station_code,
+            as_of,
+            current_level,
+            level_t_minus_3,
+            warn_thr,
+            danger_thr,
+            horizons,
         )
 
     # Per-horizon heads that were skipped at train time (e.g. too few positives) still
@@ -233,7 +278,13 @@ def _forecast_station(
         else:
             filled.extend(
                 _heuristic_forecast(
-                    station_code, as_of, current_level, level_t_minus_3, warn_thr, danger_thr, (horizon_h,)
+                    station_code,
+                    as_of,
+                    current_level,
+                    level_t_minus_3,
+                    warn_thr,
+                    danger_thr,
+                    (horizon_h,),
                 )
             )
     return filled
@@ -241,7 +292,7 @@ def _forecast_station(
 
 def get_forecasts(
     readings_by_station: Dict[str, List[dict]],
-    models_dir: Union[str, Path] = "models",
+    models_dir: Union[str, Path] = DEFAULT_MODELS_DIR,
     now: Optional[Union[datetime.datetime, str]] = None,
 ) -> List[dict]:
     """Produce flood forecasts for every station present in `readings_by_station`.
@@ -264,7 +315,9 @@ def get_forecasts(
     results: List[dict] = []
     for station_code in readings_by_station.keys():
         try:
-            results.extend(_forecast_station(station_code, grid, models_dir, now, DEFAULT_HORIZONS))
+            results.extend(
+                _forecast_station(station_code, grid, models_dir, now, DEFAULT_HORIZONS)
+            )
         except Exception as error:
             logger.error(f"Forecast failed for station {station_code}: {error}")
     return results
@@ -272,7 +325,7 @@ def get_forecasts(
 
 def get_latest_forecasts(
     db_url: Optional[str] = None,
-    models_dir: Union[str, Path] = "models",
+    models_dir: Union[str, Path] = DEFAULT_MODELS_DIR,
     hours: int = 336,
 ) -> List[dict]:
     """Convenience wrapper for web_api: load the latest window from the DB/API and forecast.
@@ -289,10 +342,14 @@ def get_latest_forecasts(
     readings_by_station: Dict[str, List[dict]] = {}
     if not df_long.empty:
         for station_code, group in df_long.groupby("station_code"):
-            readings_by_station[station_code] = group[["timestamp", "water_level", "discharge"]].to_dict("records")
+            readings_by_station[station_code] = group[
+                ["timestamp", "water_level", "discharge"]
+            ].to_dict("records")
 
     expected_stations = set(features.UPSTREAM_LEADS.keys())
     for missing_station in expected_stations - set(readings_by_station.keys()):
-        logger.warning(f"No recent data for station {missing_station}; omitting from forecasts")
+        logger.warning(
+            f"No recent data for station {missing_station}; omitting from forecasts"
+        )
 
     return get_forecasts(readings_by_station, models_dir=models_dir)

src/ml/train.py

@@ -20,8 +20,16 @@ import joblib
 import numpy as np
 import pandas as pd
 import sklearn
-from sklearn.ensemble import HistGradientBoostingClassifier, HistGradientBoostingRegressor
+from sklearn.ensemble import (
-from sklearn.metrics import average_precision_score, brier_score_loss, mean_absolute_error, mean_squared_error
+    HistGradientBoostingClassifier,
+    HistGradientBoostingRegressor,
+)
+from sklearn.metrics import (
+    average_precision_score,
+    brier_score_loss,
+    mean_absolute_error,
+    mean_squared_error,
+)
 
 from . import features
 from .data import DEFAULT_API_URL, load_measurements, resolve_db_url
@@ -51,7 +59,11 @@ HGB_PARAMS = {
 def _git_short_sha() -> str:
     try:
         result = subprocess.run(
-            ["git", "rev-parse", "--short", "HEAD"], capture_output=True, text=True, timeout=5, check=True
+            ["git", "rev-parse", "--short", "HEAD"],
+            capture_output=True,
+            text=True,
+            timeout=5,
+            check=True,
         )
         sha = result.stdout.strip()
         return sha or "nogit"
@@ -64,12 +76,20 @@ def _make_regressor(overrides: Optional[dict] = None) -> HistGradientBoostingReg
     return HistGradientBoostingRegressor(loss="squared_error", **params)
 
 
-def _make_classifier(overrides: Optional[dict] = None) -> HistGradientBoostingClassifier:
+def _make_classifier(
+    overrides: Optional[dict] = None,
+) -> HistGradientBoostingClassifier:
     params = {**HGB_PARAMS, **(overrides or {})}
     return HistGradientBoostingClassifier(**params)
 
 
-def _safe_fit(estimator, X: pd.DataFrame, y: pd.Series, head_key: str, skipped_heads: Dict[str, str]):
+def _safe_fit(
+    estimator,
+    X: pd.DataFrame,
+    y: pd.Series,
+    head_key: str,
+    skipped_heads: Dict[str, str],
+):
     """Fit an estimator, converting any failure (e.g. HistGradientBoosting's binning
     step rejecting an all-NaN/constant feature column) into a recorded skip rather
     than a station-killing exception."""
@@ -82,7 +102,9 @@ def _safe_fit(estimator, X: pd.DataFrame, y: pd.Series, head_key: str, skipped_h
         return None
 
 
-def _recall_at_far(y_true: np.ndarray, y_score: np.ndarray, target_far: float) -> Optional[float]:
+def _recall_at_far(
+    y_true: np.ndarray, y_score: np.ndarray, target_far: float
+) -> Optional[float]:
     """Recall at the score threshold whose false-positive rate over true negatives is <= target_far."""
     y_true = np.asarray(y_true)
     y_score = np.asarray(y_score)
@@ -98,7 +120,9 @@ def _recall_at_far(y_true: np.ndarray, y_score: np.ndarray, target_far: float) -
     return tp / n_pos
 
 
-def _p_warning_series(head, reg, X: pd.DataFrame, threshold: float, sigma: float) -> pd.Series:
+def _p_warning_series(
+    head, reg, X: pd.DataFrame, threshold: float, sigma: float
+) -> pd.Series:
     """Model score if a classifier head exists, else the sigmoid-derived fallback probability."""
     if head is not None:
         return pd.Series(head.predict_proba(X)[:, 1], index=X.index)
@@ -117,12 +141,24 @@ def _find_events(observed_level: pd.Series, warn_thr: float) -> List[dict]:
             start = t
         elif not is_above and start is not None:
             window = observed_level.loc[start:prev_t]
-            events.append({"crossed_warn_at": start, "peak_time": window.idxmax(), "peak_level": float(window.max())})
+            events.append(
+                {
+                    "crossed_warn_at": start,
+                    "peak_time": window.idxmax(),
+                    "peak_level": float(window.max()),
+                }
+            )
             start = None
         prev_t = t
     if start is not None:
         window = observed_level.loc[start:]
-        events.append({"crossed_warn_at": start, "peak_time": window.idxmax(), "peak_level": float(window.max())})
+        events.append(
+            {
+                "crossed_warn_at": start,
+                "peak_time": window.idxmax(),
+                "peak_level": float(window.max()),
+            }
+        )
     return events
 
 
@@ -142,9 +178,13 @@ def _events_with_lead_time(
     events = _find_events(observed_level_test, warn_thr)
     for event in events:
         first_alert_at = _first_alert_at(p_warning_test, event["crossed_warn_at"])
-        event["first_alert_at"] = first_alert_at.isoformat() if first_alert_at is not None else None
+        event["first_alert_at"] = (
+            first_alert_at.isoformat() if first_alert_at is not None else None
+        )
         if first_alert_at is not None:
-            lead_hours = (event["crossed_warn_at"] - first_alert_at).total_seconds() / 3600.0
+            lead_hours = (
+                event["crossed_warn_at"] - first_alert_at
+            ).total_seconds() / 3600.0
         else:
             lead_hours = None
         event["lead_hours"] = lead_hours
@@ -166,7 +206,10 @@ def train_station(
     """Train every head for one station. Returns (bundle_or_None, station_metrics)."""
     X, Y, meta = features.build_matrix(df_long, station, horizons)
     if meta["n_rows"] < MIN_ROWS_TO_TRAIN:
-        return None, {"status": "failed", "reason": f"only {meta['n_rows']} usable rows (< {MIN_ROWS_TO_TRAIN})"}
+        return None, {
+            "status": "failed",
+            "reason": f"only {meta['n_rows']} usable rows (< {MIN_ROWS_TO_TRAIN})",
+        }
 
     warn_thr, danger_thr = features.get_thresholds(station)
     feature_names = list(X.columns)
@@ -176,7 +219,9 @@ def train_station(
         test_mask = pd.Series(False, index=X.index)
     else:
         train_mask = X.index <= pd.Timestamp(split_train_end)
-        test_mask = (X.index >= pd.Timestamp(split_test_start)) & (X.index <= pd.Timestamp(split_test_end))
+        test_mask = (X.index >= pd.Timestamp(split_test_start)) & (
+            X.index <= pd.Timestamp(split_test_end)
+        )
     X_train, Y_train = X.loc[train_mask], Y.loc[train_mask]
     X_test, Y_test = X.loc[test_mask], Y.loc[test_mask]
     eval_X, eval_Y = (X, Y) if skip_eval else (X_train, Y_train)
@@ -188,7 +233,11 @@ def train_station(
     observed_grid = features.make_hourly_grid(df_long).observed
 
     for h in horizons:
-        max_col, warn_col, danger_col = f"max_level_{h}", f"exceed_warn_{h}", f"exceed_danger_{h}"
+        max_col, warn_col, danger_col = (
+            f"max_level_{h}",
+            f"exceed_warn_{h}",
+            f"exceed_danger_{h}",
+        )
         horizon_metrics: dict = {}
 
         # --- regression head (max level) ---
@@ -215,19 +264,28 @@ def train_station(
                 sigma_h = max(float(np.std(residuals)), MIN_SIGMA)
                 horizon_metrics["n_test"] = int(test_labeled.sum())
                 horizon_metrics["mae"] = float(mean_absolute_error(y_true, y_pred))
-                horizon_metrics["rmse"] = float(np.sqrt(mean_squared_error(y_true, y_pred)))
+                horizon_metrics["rmse"] = float(
+                    np.sqrt(mean_squared_error(y_true, y_pred))
+                )
                 above_2m = y_true >= 2.0
                 horizon_metrics["mae_above_2m"] = (
-                    float(mean_absolute_error(y_true[above_2m], y_pred[above_2m])) if above_2m.any() else None
+                    float(mean_absolute_error(y_true[above_2m], y_pred[above_2m]))
+                    if above_2m.any()
+                    else None
                 )
         sigma[h] = sigma_h
         horizon_metrics["sigma"] = sigma_h
 
         # --- classification heads (warn / danger) ---
         p_warning_test = None
-        for label_name, col, thr in (("warn", warn_col, warn_thr), ("danger", danger_col, danger_thr)):
+        for label_name, col, thr in (
+            ("warn", warn_col, warn_thr),
+            ("danger", danger_col, danger_thr),
+        ):
             train_labeled = eval_Y[col].notna()
-            n_pos = int(eval_Y.loc[train_labeled, col].sum()) if train_labeled.any() else 0
+            n_pos = (
+                int(eval_Y.loc[train_labeled, col].sum()) if train_labeled.any() else 0
+            )
             head_key = f"{label_name}_{h}"
             clf = None
             if n_pos >= MIN_POSITIVES_FOR_CLASSIFIER:
@@ -239,21 +297,37 @@ def train_station(
                     skipped_heads,
                 )
             else:
-                skipped_heads[head_key] = f"only {n_pos} positives in train span (< {MIN_POSITIVES_FOR_CLASSIFIER})"
+                skipped_heads[
+                    head_key
+                ] = f"only {n_pos} positives in train span (< {MIN_POSITIVES_FOR_CLASSIFIER})"
             heads[head_key] = clf
 
             if not skip_eval:
                 test_labeled = Y_test[col].notna()
                 horizon_metrics[f"base_rate_{label_name}"] = (
-                    float(Y_test.loc[test_labeled, col].mean()) if test_labeled.any() else None
+                    float(Y_test.loc[test_labeled, col].mean())
+                    if test_labeled.any()
+                    else None
                 )
-                if clf is not None and test_labeled.sum() > 0 and Y_test.loc[test_labeled, col].nunique() > 1:
+                if (
+                    clf is not None
+                    and test_labeled.sum() > 0
+                    and Y_test.loc[test_labeled, col].nunique() > 1
+                ):
                     y_true = Y_test.loc[test_labeled, col]
                     y_score = clf.predict_proba(X_test.loc[test_labeled])[:, 1]
-                    horizon_metrics[f"pr_auc_{label_name}"] = float(average_precision_score(y_true, y_score))
+                    horizon_metrics[f"pr_auc_{label_name}"] = float(
-                    horizon_metrics[f"brier_{label_name}"] = float(brier_score_loss(y_true, y_score))
+                        average_precision_score(y_true, y_score)
-                    horizon_metrics[f"recall_{label_name}_at_far1pct"] = _recall_at_far(y_true, y_score, 0.01)
+                    )
-                    horizon_metrics[f"recall_{label_name}_at_far5pct"] = _recall_at_far(y_true, y_score, 0.05)
+                    horizon_metrics[f"brier_{label_name}"] = float(
+                        brier_score_loss(y_true, y_score)
+                    )
+                    horizon_metrics[f"recall_{label_name}_at_far1pct"] = _recall_at_far(
+                        y_true, y_score, 0.01
+                    )
+                    horizon_metrics[f"recall_{label_name}_at_far5pct"] = _recall_at_far(
+                        y_true, y_score, 0.05
+                    )
                 else:
                     horizon_metrics[f"pr_auc_{label_name}"] = None
                     horizon_metrics[f"brier_{label_name}"] = None
@@ -269,8 +343,12 @@ def train_station(
         if not skip_eval and reg is not None and p_warning_test is not None:
             observed_test_level = observed_grid.get((station, "water_level"))
             if observed_test_level is not None:
-                observed_test_level = observed_test_level.loc[observed_test_level.index.isin(X_test.index)]
+                observed_test_level = observed_test_level.loc[
-                per_horizon[h]["events"] = _events_with_lead_time(observed_test_level, warn_thr, p_warning_test)
+                    observed_test_level.index.isin(X_test.index)
+                ]
+                per_horizon[h]["events"] = _events_with_lead_time(
+                    observed_test_level, warn_thr, p_warning_test
+                )
 
     # --- full refit on the ENTIRE record for the deployed artifact ---
     # This may include/exclude different heads than the eval-phase gate above (the
@@ -278,11 +356,21 @@ def train_station(
     # skipped_heads must reflect what actually ends up in the saved bundle.
     final_heads: Dict[str, object] = {}
     for h in horizons:
-        max_col, warn_col, danger_col = f"max_level_{h}", f"exceed_warn_{h}", f"exceed_danger_{h}"
+        max_col, warn_col, danger_col = (
+            f"max_level_{h}",
+            f"exceed_warn_{h}",
+            f"exceed_danger_{h}",
+        )
         head_key = f"max_{h}"
         labeled = Y[max_col].notna()
         if labeled.sum() >= MIN_ROWS_FOR_HEAD:
-            reg = _safe_fit(_make_regressor(hgb_overrides), X.loc[labeled], Y.loc[labeled, max_col], head_key, skipped_heads)
+            reg = _safe_fit(
+                _make_regressor(hgb_overrides),
+                X.loc[labeled],
+                Y.loc[labeled, max_col],
+                head_key,
+                skipped_heads,
+            )
             final_heads[head_key] = reg
             if reg is not None:
                 skipped_heads.pop(head_key, None)
@@ -296,13 +384,19 @@ def train_station(
             n_pos = int(Y.loc[train_labeled, col].sum()) if train_labeled.any() else 0
             if n_pos >= MIN_POSITIVES_FOR_CLASSIFIER:
                 clf = _safe_fit(
-                    _make_classifier(hgb_overrides), X.loc[train_labeled], Y.loc[train_labeled, col], head_key, skipped_heads
+                    _make_classifier(hgb_overrides),
+                    X.loc[train_labeled],
+                    Y.loc[train_labeled, col],
+                    head_key,
+                    skipped_heads,
                 )
                 final_heads[head_key] = clf
                 if clf is not None:
                     skipped_heads.pop(head_key, None)
             else:
-                skipped_heads[head_key] = f"only {n_pos} positives in train span (< {MIN_POSITIVES_FOR_CLASSIFIER})"
+                skipped_heads[
+                    head_key
+                ] = f"only {n_pos} positives in train span (< {MIN_POSITIVES_FOR_CLASSIFIER})"
                 final_heads[head_key] = None
 
     bundle = {
@@ -345,7 +439,11 @@ def train_all(
             continue
         try:
             bundle, station_metrics = train_station(
-                df_long, station, horizons, skip_eval=skip_eval, hgb_overrides=hgb_overrides
+                df_long,
+                station,
+                horizons,
+                skip_eval=skip_eval,
+                hgb_overrides=hgb_overrides,
             )
             if bundle is None:
                 logger.warning(f"{station}: failed ({station_metrics.get('reason')})")
@@ -377,16 +475,32 @@ def train_all(
 
 
 def main(argv: Optional[List[str]] = None) -> None:
-    logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(name)s: %(message)s")
+    logging.basicConfig(
+        level=logging.INFO, format="%(asctime)s %(levelname)s %(name)s: %(message)s"
+    )
 
-    parser = argparse.ArgumentParser(description="Train Ping River flood forecast models")
+    parser = argparse.ArgumentParser(
-    parser.add_argument("--stations", default="all", help="'all' or a comma-separated list of station codes")
+        description="Train Ping River flood forecast models"
+    )
+    parser.add_argument(
+        "--stations",
+        default="all",
+        help="'all' or a comma-separated list of station codes",
+    )
     parser.add_argument("--models-dir", default="models")
     parser.add_argument("--db-url", default=None)
     parser.add_argument("--api-url", default=DEFAULT_API_URL)
-    parser.add_argument("--skip-eval", action="store_true", help="Refit-only fast path; skip Split B evaluation")
+    parser.add_argument(
-    parser.add_argument("--start", default=None, help="ISO date; earliest measurement to load")
+        "--skip-eval",
-    parser.add_argument("--end", default=None, help="ISO date; latest measurement to load")
+        action="store_true",
+        help="Refit-only fast path; skip Split B evaluation",
+    )
+    parser.add_argument(
+        "--start", default=None, help="ISO date; earliest measurement to load"
+    )
+    parser.add_argument(
+        "--end", default=None, help="ISO date; latest measurement to load"
+    )
     args = parser.parse_args(argv)
 
     if args.stations == "all":
@@ -399,15 +513,25 @@ def main(argv: Optional[List[str]] = None) -> None:
 
     logger.info(f"Loading measurements for {len(stations)} stations...")
     df_long = load_measurements(
-        db_url=resolve_db_url(args.db_url), stations=None, start=start, end=end, api_url=args.api_url
+        db_url=resolve_db_url(args.db_url),
+        stations=None,
+        start=start,
+        end=end,
+        api_url=args.api_url,
+    )
+    logger.info(
+        f"Loaded {len(df_long)} rows spanning {df_long['timestamp'].min()} .. {df_long['timestamp'].max()}"
     )
-    logger.info(f"Loaded {len(df_long)} rows spanning {df_long['timestamp'].min()} .. {df_long['timestamp'].max()}")
 
     metrics_payload = train_all(
         df_long, stations, models_dir=Path(args.models_dir), skip_eval=args.skip_eval
     )
-    trained = sum(1 for s in metrics_payload["stations"].values() if s["status"] == "trained")
+    trained = sum(
-    logger.info(f"Done: {trained}/{len(stations)} stations trained. metrics.json written to {args.models_dir}")
+        1 for s in metrics_payload["stations"].values() if s["status"] == "trained"
+    )
+    logger.info(
+        f"Done: {trained}/{len(stations)} stations trained. metrics.json written to {args.models_dir}"
+    )
 
 
 if __name__ == "__main__":