diff --git a/src/ml/data.py b/src/ml/data.py index 6684e79..eb0834f 100644 --- a/src/ml/data.py +++ b/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}") diff --git a/src/ml/features.py b/src/ml/features.py index 45bef6d..28b0180 100644 --- a/src/ml/features.py +++ b/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.103": [("P.67", 6), ("P.21", 8), ("P.75", 11), ("P.4A", 11), ("P.92", 14), ("P.20", 16)], + "P.1": [ + ("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, } diff --git a/src/ml/predict.py b/src/ml/predict.py index a0d5602..9e6fe5e 100644 --- a/src/ml/predict.py +++ b/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_danger = _clip_probability(_sigmoid_probability(predicted_max, danger_thr, HEURISTIC_SIGMA)) + p_warning = _clip_probability( + _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) diff --git a/src/ml/train.py b/src/ml/train.py index 49331c6..b3b95ac 100644 --- a/src/ml/train.py +++ b/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.metrics import average_precision_score, brier_score_loss, mean_absolute_error, mean_squared_error +from sklearn.ensemble import ( + 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"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) + horizon_metrics[f"pr_auc_{label_name}"] = float( + average_precision_score(y_true, y_score) + ) + 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)] - per_horizon[h]["events"] = _events_with_lead_time(observed_test_level, warn_thr, p_warning_test) + observed_test_level = observed_test_level.loc[ + 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.add_argument("--stations", default="all", help="'all' or a comma-separated list of station codes") + parser = argparse.ArgumentParser( + 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("--start", default=None, help="ISO date; earliest measurement to load") - parser.add_argument("--end", default=None, help="ISO date; latest measurement to load") + parser.add_argument( + "--skip-eval", + 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") - logger.info(f"Done: {trained}/{len(stations)} stations trained. metrics.json written to {args.models_dir}") + trained = sum( + 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__":