删除 LGBM 全部代码和模型文件,EMOS 简化为纯 legacy 高斯分桶模式

This commit is contained in:
2569718930@qq.com
2026-05-18 22:05:55 +08:00
parent aec47adda1
commit 0e0aad3171
26 changed files with 11 additions and 149413 deletions
-834
View File
@@ -1,834 +0,0 @@
from __future__ import annotations
import json
import math
import os
from datetime import datetime, timezone
from typing import Any, Dict, Iterable, List, Optional, Tuple
import numpy as np
from src.analysis.settlement_rounding import apply_city_settlement, is_exact_settlement_city
ENGINE_MODE_LEGACY = "legacy"
ENGINE_MODE_EMOS_SHADOW = "emos_shadow"
ENGINE_MODE_EMOS_PRIMARY = "emos_primary"
DEFAULT_ENGINE_MODE = ENGINE_MODE_EMOS_PRIMARY
VALID_ENGINE_MODES = {
ENGINE_MODE_LEGACY,
ENGINE_MODE_EMOS_SHADOW,
ENGINE_MODE_EMOS_PRIMARY,
}
DEFAULT_CALIBRATION_FILE = os.path.join(
os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))),
"artifacts",
"probability_calibration",
"default.json",
)
_CALIBRATION_CACHE: Dict[str, Dict[str, Any]] = {}
_CALIBRATION_MTIME: Dict[str, float] = {}
def _sf(value: Any) -> Optional[float]:
if value is None:
return None
try:
return float(value)
except Exception:
return None
def _coalesce_float(value: Any, default: float) -> float:
parsed = _sf(value)
return default if parsed is None else parsed
def _mean(values: Iterable[float]) -> Optional[float]:
values = list(values)
return (sum(values) / len(values)) if values else None
def resolve_probability_engine_mode(explicit_mode: Optional[str] = None) -> str:
mode = str(
explicit_mode
or os.getenv("POLYWEATHER_PROBABILITY_ENGINE")
or DEFAULT_ENGINE_MODE
).strip().lower()
if mode not in VALID_ENGINE_MODES:
return DEFAULT_ENGINE_MODE
return mode
def load_calibration(calibration_path: Optional[str] = None) -> Dict[str, Any]:
path = str(
calibration_path
or os.getenv("POLYWEATHER_PROBABILITY_CALIBRATION_FILE")
or DEFAULT_CALIBRATION_FILE
).strip()
if not path:
return {}
if not os.path.exists(path):
return {}
try:
mtime = os.path.getmtime(path)
cached = _CALIBRATION_CACHE.get(path)
if cached and _CALIBRATION_MTIME.get(path) == mtime:
return cached
with open(path, "r", encoding="utf-8") as fh:
data = json.load(fh)
if not isinstance(data, dict):
return {}
_CALIBRATION_CACHE[path] = data
_CALIBRATION_MTIME[path] = mtime
return data
except Exception:
return {}
def build_probability_features(
city_name: str,
raw_mu: Optional[float],
raw_sigma: Optional[float],
deb_prediction: Optional[float],
ens_data: Optional[Dict[str, Any]],
current_forecasts: Optional[Dict[str, Any]],
max_so_far: Optional[float],
peak_status: str,
local_hour_frac: Optional[float],
) -> Dict[str, Any]:
ens_data = ens_data or {}
current_forecasts = current_forecasts or {}
forecast_values = [
v for v in (_sf(val) for val in current_forecasts.values()) if v is not None
]
forecast_values.sort()
forecast_median = None
if forecast_values:
forecast_median = forecast_values[len(forecast_values) // 2]
ens_median = _sf(ens_data.get("median"))
ens_p10 = _sf(ens_data.get("p10"))
ens_p90 = _sf(ens_data.get("p90"))
ensemble_spread = None
if ens_p10 is not None and ens_p90 is not None and ens_p90 >= ens_p10:
ensemble_spread = max(0.1, (ens_p90 - ens_p10) / 2.56)
elif len(forecast_values) >= 2:
ensemble_spread = max(0.6, (forecast_values[-1] - forecast_values[0]) / 2.0)
elif raw_sigma is not None:
ensemble_spread = max(0.1, raw_sigma)
baseline = deb_prediction if deb_prediction is not None else raw_mu
max_so_far_gap = None
if baseline is not None and max_so_far is not None:
max_so_far_gap = baseline - max_so_far
peak_flag = 0.0
if peak_status == "in_window":
peak_flag = 0.5
elif peak_status == "past":
peak_flag = 1.0
return {
"city": str(city_name or "").strip().lower(),
"raw_mu": raw_mu,
"raw_sigma": raw_sigma,
"deb_prediction": deb_prediction,
"ens_median": ens_median,
"ens_p10": ens_p10,
"ens_p90": ens_p90,
"forecast_median": forecast_median,
"forecast_spread": forecast_values[-1] - forecast_values[0]
if len(forecast_values) >= 2
else None,
"ensemble_spread": ensemble_spread,
"max_so_far": max_so_far,
"max_so_far_gap": max_so_far_gap,
"peak_status": peak_status,
"peak_flag": peak_flag,
"local_hour_frac": local_hour_frac,
"model_count": len(forecast_values),
}
def _normal_cdf(x: float, mean: float, sigma: float) -> float:
return 0.5 * (1.0 + math.erf((x - mean) / (sigma * math.sqrt(2.0))))
def _normal_pdf(x: float) -> float:
return math.exp(-(x ** 2) / 2.0) / math.sqrt(2.0 * math.pi)
def _bucket_probabilities(
mu: float,
sigma: float,
max_so_far: Optional[float],
city_name: str,
) -> Tuple[List[Dict[str, Any]], List[Tuple[int, float]], List[Dict[str, Any]]]:
sigma = max(0.1, float(sigma))
min_possible = (
apply_city_settlement(city_name, max_so_far) if max_so_far is not None else -999
)
probs: Dict[int, float] = {}
search_range = max(2, int(sigma * 2.5))
is_exact = is_exact_settlement_city(city_name)
target_mu = apply_city_settlement(city_name, mu)
if is_exact:
target_mu = int(math.floor(mu))
for value in range(target_mu - search_range, target_mu + search_range + 1):
if value < min_possible:
continue
if is_exact:
prob = _normal_cdf(value + 1.0, mu, sigma) - _normal_cdf(value, mu, sigma)
else:
prob = _normal_cdf(value + 0.5, mu, sigma) - _normal_cdf(value - 0.5, mu, sigma)
if prob > 0.01:
probs[value] = prob
total = sum(probs.values())
if total <= 0:
return [], [], []
normalized = {key: val / total for key, val in probs.items()}
sorted_probs = sorted(normalized.items(), key=lambda item: item[1], reverse=True)
full_distribution = []
for value, prob in sorted_probs:
if is_exact:
bucket_range = "[{0}.0~{1}.0)".format(value, value + 1)
else:
bucket_range = "[{0}~{1})".format(value - 0.5, value + 0.5)
full_distribution.append(
{
"value": int(value),
"range": bucket_range,
"probability": round(prob, 3),
}
)
return full_distribution[:4], sorted_probs, full_distribution
def _top_bucket_value(distribution: Optional[List[Dict[str, Any]]]) -> Optional[int]:
if not distribution:
return None
top = max(
(row for row in distribution if isinstance(row, dict)),
key=lambda row: float(row.get("probability") or 0.0),
default=None,
)
if not top:
return None
value = top.get("value")
return int(value) if value is not None else None
def _bucket_brier_score(
distribution: Optional[List[Dict[str, Any]]],
city_name: str,
actual_high: float,
) -> float:
actual_bucket = apply_city_settlement(city_name, actual_high)
hit_prob = 0.0
total = 0.0
for row in distribution or []:
if not isinstance(row, dict):
continue
value = row.get("value")
try:
prob = float(row.get("probability") or 0.0)
except Exception:
prob = 0.0
if value == actual_bucket:
hit_prob = prob
else:
total += prob * prob
total += (1.0 - hit_prob) ** 2
return total
def _composite_score(
mean_crps: float,
mean_mae: float,
bucket_hit_rate: float,
bucket_brier: float,
) -> float:
return mean_crps + 0.1 * mean_mae + 1.5 * (1.0 - bucket_hit_rate) + 0.75 * bucket_brier
def _blend_value(raw_value: float, calibrated_value: float, alpha: float) -> float:
return (1.0 - alpha) * raw_value + alpha * calibrated_value
def _clamp_sigma(
raw_sigma: float,
calibrated_sigma: float,
constraints: Optional[Dict[str, Any]],
) -> float:
constraints = constraints or {}
min_ratio = max(0.25, _coalesce_float(constraints.get("min_ratio"), 0.85))
max_ratio = max(min_ratio, _coalesce_float(constraints.get("max_ratio"), 1.35))
absolute_min = max(0.1, _coalesce_float(constraints.get("absolute_min"), 0.25))
absolute_max = max(
absolute_min,
_coalesce_float(constraints.get("absolute_max"), raw_sigma * max_ratio),
)
sigma_floor = max(absolute_min, raw_sigma * min_ratio)
sigma_cap = min(absolute_max, raw_sigma * max_ratio)
return min(max(calibrated_sigma, sigma_floor), sigma_cap)
def apply_probability_calibration(
city_name: str,
temp_symbol: str,
raw_mu: Optional[float],
raw_sigma: Optional[float],
max_so_far: Optional[float],
legacy_distribution: Optional[List[Dict[str, Any]]],
features: Optional[Dict[str, Any]] = None,
calibration_path: Optional[str] = None,
mode: Optional[str] = None,
) -> Dict[str, Any]:
selected_mode = resolve_probability_engine_mode(mode)
if raw_mu is None or raw_sigma is None:
return {
"mode": selected_mode,
"engine": ENGINE_MODE_LEGACY,
"distribution": legacy_distribution or [],
"distribution_all": legacy_distribution or [],
"shadow_distribution": [],
"shadow_distribution_all": [],
"raw_mu": raw_mu,
"raw_sigma": raw_sigma,
"calibrated_mu": None,
"calibrated_sigma": None,
"calibration_version": None,
"calibration_source": None,
}
calibration = load_calibration(calibration_path)
if not calibration:
return {
"mode": selected_mode,
"engine": ENGINE_MODE_LEGACY,
"distribution": legacy_distribution or [],
"distribution_all": legacy_distribution or [],
"shadow_distribution": [],
"shadow_distribution_all": [],
"raw_mu": raw_mu,
"raw_sigma": raw_sigma,
"calibrated_mu": None,
"calibrated_sigma": None,
"calibration_version": None,
"calibration_source": None,
}
features = features or {}
city_key = str(city_name or "").strip().lower()
global_params = calibration.get("global", {}) or {}
city_params = (calibration.get("cities", {}) or {}).get(city_key, {}) or {}
blending_cfg = calibration.get("blending", {}) or {}
sigma_constraints = calibration.get("sigma_constraints", {}) or {}
mu_cfg = global_params.get("mu", {}) or {}
sigma_cfg = global_params.get("sigma", {}) or {}
city_confidence = max(0.0, min(1.0, _coalesce_float(city_params.get("confidence"), 1.0)))
city_mu_bias = _coalesce_float(city_params.get("mu_bias"), 0.0) * city_confidence
city_sigma_scale = 1.0 + (
(_coalesce_float(city_params.get("sigma_scale"), 1.0) - 1.0) * city_confidence
)
deb_prediction = _sf(features.get("deb_prediction"))
ens_median = _sf(features.get("ens_median"))
max_so_far_gap = _sf(features.get("max_so_far_gap"))
peak_flag = _sf(features.get("peak_flag")) or 0.0
ensemble_spread = _sf(features.get("ensemble_spread"))
mu_intercept = _coalesce_float(mu_cfg.get("intercept"), 0.0)
mu_raw_coef = _coalesce_float(mu_cfg.get("raw_mu_coef"), 1.0)
mu_deb_coef = _coalesce_float(mu_cfg.get("deb_coef"), 0.0)
mu_ens_coef = _coalesce_float(mu_cfg.get("ens_median_coef"), 0.0)
mu_gap_coef = _coalesce_float(mu_cfg.get("max_so_far_gap_coef"), 0.0)
sigma_intercept = _coalesce_float(
sigma_cfg.get("intercept"),
math.log(max(raw_sigma, 0.1)),
)
sigma_raw_coef = _coalesce_float(sigma_cfg.get("raw_sigma_coef"), 1.0)
sigma_spread_coef = _coalesce_float(sigma_cfg.get("spread_coef"), 0.0)
sigma_peak_coef = _coalesce_float(sigma_cfg.get("peak_flag_coef"), 0.0)
sigma_gap_coef = _coalesce_float(sigma_cfg.get("max_so_far_gap_coef"), 0.0)
calibrated_mu = (
mu_intercept
+ mu_raw_coef * raw_mu
+ mu_deb_coef * (deb_prediction if deb_prediction is not None else raw_mu)
+ mu_ens_coef * (ens_median if ens_median is not None else raw_mu)
+ mu_gap_coef * (max_so_far_gap if max_so_far_gap is not None else 0.0)
+ city_mu_bias
)
sigma_log = (
sigma_intercept
+ sigma_raw_coef * math.log(max(raw_sigma, 0.1))
+ sigma_spread_coef * math.log(max(ensemble_spread or raw_sigma, 0.1))
+ sigma_peak_coef * peak_flag
+ sigma_gap_coef * (max_so_far_gap if max_so_far_gap is not None else 0.0)
)
calibrated_sigma = max(0.1, math.exp(sigma_log) * city_sigma_scale)
blend_alpha_mu = max(0.0, min(1.0, _coalesce_float(blending_cfg.get("alpha_mu"), 1.0)))
blend_alpha_sigma = max(0.0, min(1.0, _coalesce_float(blending_cfg.get("alpha_sigma"), 1.0)))
calibrated_mu = _blend_value(raw_mu, calibrated_mu, blend_alpha_mu)
if max_so_far is not None:
observed_floor = _sf(max_so_far)
if observed_floor is not None and calibrated_mu < observed_floor:
calibrated_mu = observed_floor
calibrated_sigma = max(0.1, _blend_value(raw_sigma, calibrated_sigma, blend_alpha_sigma))
calibrated_sigma = _clamp_sigma(raw_sigma, calibrated_sigma, sigma_constraints)
calibrated_distribution, calibrated_sorted, calibrated_distribution_all = _bucket_probabilities(
calibrated_mu,
calibrated_sigma,
max_so_far=max_so_far,
city_name=city_key,
)
engine = ENGINE_MODE_LEGACY
selected_distribution = legacy_distribution or []
selected_distribution_all = legacy_distribution or []
selected_sorted: List[Tuple[int, float]] = []
shadow_distribution: List[Dict[str, Any]] = []
shadow_distribution_all: List[Dict[str, Any]] = []
shadow_sorted: List[Tuple[int, float]] = []
if selected_mode == ENGINE_MODE_EMOS_PRIMARY:
engine = "emos"
selected_distribution = calibrated_distribution
selected_distribution_all = calibrated_distribution_all
selected_sorted = calibrated_sorted
elif selected_mode == ENGINE_MODE_EMOS_SHADOW:
shadow_distribution = calibrated_distribution
shadow_distribution_all = calibrated_distribution_all
shadow_sorted = calibrated_sorted
return {
"mode": selected_mode,
"engine": engine,
"distribution": selected_distribution,
"distribution_all": selected_distribution_all,
"selected_sorted_probs": selected_sorted,
"shadow_distribution": shadow_distribution,
"shadow_distribution_all": shadow_distribution_all,
"shadow_sorted_probs": shadow_sorted,
"raw_mu": raw_mu,
"raw_sigma": raw_sigma,
"calibrated_mu": calibrated_mu,
"calibrated_sigma": calibrated_sigma,
"blend_alpha_mu": blend_alpha_mu,
"blend_alpha_sigma": blend_alpha_sigma,
"calibration_version": calibration.get("version"),
"calibration_source": calibration.get("source")
or os.path.relpath(
calibration_path or DEFAULT_CALIBRATION_FILE,
os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))),
),
}
def _gaussian_crps(observation: float, mean: float, sigma: float) -> float:
sigma = max(0.1, float(sigma))
z = (observation - mean) / sigma
return sigma * (
z * (2.0 * _normal_cdf(z, 0.0, 1.0) - 1.0)
+ 2.0 * _normal_pdf(z)
- 1.0 / math.sqrt(math.pi)
)
def _fit_linear(xs: np.ndarray, ys: np.ndarray) -> np.ndarray:
if len(xs) == 0:
return np.zeros(xs.shape[1], dtype=float)
coeffs, _, _, _ = np.linalg.lstsq(xs, ys, rcond=None)
return coeffs
def fit_calibration(
samples: Iterable[Dict[str, Any]],
version: Optional[str] = None,
) -> Dict[str, Any]:
normalized_samples: List[Dict[str, Any]] = []
city_residuals: Dict[str, List[float]] = {}
city_sigma_ratios: Dict[str, List[float]] = {}
for raw_sample in samples:
actual = _sf(raw_sample.get("actual_high"))
raw_mu = _sf(raw_sample.get("raw_mu"))
raw_sigma = _sf(raw_sample.get("raw_sigma"))
if actual is None or raw_mu is None or raw_sigma is None:
continue
feature_row = {
"city": str(raw_sample.get("city") or "").strip().lower(),
"actual_high": actual,
"raw_mu": raw_mu,
"raw_sigma": max(0.1, raw_sigma),
"deb_prediction": _sf(raw_sample.get("deb_prediction")),
"ens_median": _sf(raw_sample.get("ens_median")),
"ensemble_spread": _sf(raw_sample.get("ensemble_spread")),
"max_so_far_gap": _sf(raw_sample.get("max_so_far_gap")),
"peak_flag": _sf(raw_sample.get("peak_flag")) or 0.0,
}
normalized_samples.append(feature_row)
if len(normalized_samples) < 3:
return default_calibration_payload(version=version, reason="insufficient_samples")
mu_rows = []
mu_targets = []
for sample in normalized_samples:
deb = sample["deb_prediction"] if sample["deb_prediction"] is not None else sample["raw_mu"]
ens_median = sample["ens_median"] if sample["ens_median"] is not None else sample["raw_mu"]
gap = sample["max_so_far_gap"] if sample["max_so_far_gap"] is not None else 0.0
mu_rows.append([1.0, sample["raw_mu"], deb, ens_median, gap])
mu_targets.append(sample["actual_high"])
mu_coeffs = _fit_linear(np.array(mu_rows, dtype=float), np.array(mu_targets, dtype=float))
sigma_rows = []
sigma_targets = []
mu_predictions = []
for idx, sample in enumerate(normalized_samples):
predicted_mu = float(np.dot(mu_coeffs, np.array(mu_rows[idx], dtype=float)))
mu_predictions.append(predicted_mu)
residual = max(abs(sample["actual_high"] - predicted_mu), 0.1)
spread = max(sample["ensemble_spread"] or sample["raw_sigma"], 0.1)
gap = sample["max_so_far_gap"] if sample["max_so_far_gap"] is not None else 0.0
sigma_rows.append([1.0, math.log(sample["raw_sigma"]), math.log(spread), sample["peak_flag"], gap])
sigma_targets.append(math.log(residual))
city_residuals.setdefault(sample["city"], []).append(sample["actual_high"] - predicted_mu)
city_sigma_ratios.setdefault(sample["city"], []).append(residual / max(sample["raw_sigma"], 0.1))
sigma_coeffs = _fit_linear(np.array(sigma_rows, dtype=float), np.array(sigma_targets, dtype=float))
crps_values = []
for idx, sample in enumerate(normalized_samples):
predicted_mu = mu_predictions[idx]
sigma_log = float(np.dot(sigma_coeffs, np.array(sigma_rows[idx], dtype=float)))
predicted_sigma = max(0.1, math.exp(sigma_log))
crps_values.append(_gaussian_crps(sample["actual_high"], predicted_mu, predicted_sigma))
city_params: Dict[str, Dict[str, Any]] = {}
for city, residuals in city_residuals.items():
if len(residuals) < 3:
continue
sigma_ratios = city_sigma_ratios.get(city) or [1.0]
confidence = max(0.25, min(1.0, len(residuals) / 8.0))
city_params[city] = {
"samples": len(residuals),
"mu_bias": round(sum(residuals) / len(residuals), 6),
"sigma_scale": round(
max(0.5, min(2.0, sum(sigma_ratios) / len(sigma_ratios))),
6,
),
"confidence": round(confidence, 6),
}
legacy_crps_values = []
legacy_mae_values = []
legacy_bucket_hits = []
legacy_bucket_briers = []
candidate_predictions = []
for idx, sample in enumerate(normalized_samples):
city = sample["city"]
city_meta = city_params.get(city, {})
city_confidence = max(
0.0,
min(1.0, _coalesce_float(city_meta.get("confidence"), 1.0)),
)
city_mu_bias = _coalesce_float(city_meta.get("mu_bias"), 0.0) * city_confidence
city_sigma_scale = 1.0 + (
(_coalesce_float(city_meta.get("sigma_scale"), 1.0) - 1.0) * city_confidence
)
legacy_mu = sample["raw_mu"]
legacy_sigma = sample["raw_sigma"]
actual_high = sample["actual_high"]
legacy_crps_values.append(_gaussian_crps(actual_high, legacy_mu, legacy_sigma))
legacy_mae_values.append(abs(legacy_mu - actual_high))
legacy_bucket_hits.append(
1.0
if apply_city_settlement(city, legacy_mu)
== apply_city_settlement(city, actual_high)
else 0.0
)
legacy_distribution, _, _ = _bucket_probabilities(
legacy_mu,
legacy_sigma,
max_so_far=None,
city_name=city,
)
legacy_bucket_briers.append(
_bucket_brier_score(legacy_distribution, city, actual_high)
)
calibrated_mu = mu_predictions[idx] + city_mu_bias
sigma_log = float(np.dot(sigma_coeffs, np.array(sigma_rows[idx], dtype=float)))
calibrated_sigma = max(0.1, math.exp(sigma_log) * city_sigma_scale)
candidate_predictions.append(
{
"city": city,
"actual_high": actual_high,
"raw_mu": legacy_mu,
"raw_sigma": legacy_sigma,
"calibrated_mu": calibrated_mu,
"calibrated_sigma": calibrated_sigma,
}
)
legacy_mean_crps = _mean(legacy_crps_values) or 0.0
legacy_mean_mae = _mean(legacy_mae_values) or 0.0
legacy_bucket_hit_rate = _mean(legacy_bucket_hits) or 0.0
legacy_bucket_brier = _mean(legacy_bucket_briers) or 0.0
sigma_constraints = {
"min_ratio": 0.85,
"max_ratio": 1.35,
"absolute_min": 0.25,
"absolute_max": 3.0,
}
legacy_score = _composite_score(
legacy_mean_crps,
legacy_mean_mae,
legacy_bucket_hit_rate,
legacy_bucket_brier,
)
guardrails = {
"max_mae_increase": 0.02,
"max_bucket_hit_drop": 0.01,
"max_bucket_brier_increase": 0.05,
}
best_alpha_mu = 0.0
best_alpha_sigma = 0.0
best_score = legacy_score
best_metrics = {
"mean_crps": legacy_mean_crps,
"mean_mae": legacy_mean_mae,
"bucket_hit_rate": legacy_bucket_hit_rate,
"bucket_brier": legacy_bucket_brier,
}
alpha_grid = [step / 20.0 for step in range(21)]
for alpha_mu in alpha_grid:
for alpha_sigma in alpha_grid:
crps_values = []
mae_values = []
bucket_hits = []
bucket_briers = []
for row in candidate_predictions:
mu_hat = _blend_value(row["raw_mu"], row["calibrated_mu"], alpha_mu)
sigma_hat = _clamp_sigma(
row["raw_sigma"],
max(
0.1,
_blend_value(row["raw_sigma"], row["calibrated_sigma"], alpha_sigma),
),
sigma_constraints,
)
actual_high = row["actual_high"]
city = row["city"]
crps_values.append(_gaussian_crps(actual_high, mu_hat, sigma_hat))
mae_values.append(abs(mu_hat - actual_high))
distribution, _, _ = _bucket_probabilities(
mu_hat,
sigma_hat,
max_so_far=None,
city_name=city,
)
predicted_bucket = _top_bucket_value(distribution)
actual_bucket = apply_city_settlement(city, actual_high)
bucket_hits.append(1.0 if predicted_bucket == actual_bucket else 0.0)
bucket_briers.append(
_bucket_brier_score(distribution, city, actual_high)
)
mean_crps = _mean(crps_values) or 0.0
mean_mae = _mean(mae_values) or 0.0
bucket_hit_rate = _mean(bucket_hits) or 0.0
bucket_brier = _mean(bucket_briers) or 0.0
if mean_mae > legacy_mean_mae + guardrails["max_mae_increase"]:
continue
if bucket_hit_rate + guardrails["max_bucket_hit_drop"] < legacy_bucket_hit_rate:
continue
if bucket_brier > legacy_bucket_brier + guardrails["max_bucket_brier_increase"]:
continue
score = _composite_score(mean_crps, mean_mae, bucket_hit_rate, bucket_brier)
if score + 1e-9 < best_score:
best_score = score
best_alpha_mu = alpha_mu
best_alpha_sigma = alpha_sigma
best_metrics = {
"mean_crps": mean_crps,
"mean_mae": mean_mae,
"bucket_hit_rate": bucket_hit_rate,
"bucket_brier": bucket_brier,
}
return {
"version": version or datetime.now(timezone.utc).strftime("emos-%Y%m%d%H%M%S"),
"trained_at": datetime.now(timezone.utc).isoformat(),
"global": {
"mu": {
"intercept": round(float(mu_coeffs[0]), 8),
"raw_mu_coef": round(float(mu_coeffs[1]), 8),
"deb_coef": round(float(mu_coeffs[2]), 8),
"ens_median_coef": round(float(mu_coeffs[3]), 8),
"max_so_far_gap_coef": round(float(mu_coeffs[4]), 8),
},
"sigma": {
"intercept": round(float(sigma_coeffs[0]), 8),
"raw_sigma_coef": round(float(sigma_coeffs[1]), 8),
"spread_coef": round(float(sigma_coeffs[2]), 8),
"peak_flag_coef": round(float(sigma_coeffs[3]), 8),
"max_so_far_gap_coef": round(float(sigma_coeffs[4]), 8),
},
},
"sigma_constraints": sigma_constraints,
"selection_guardrails": guardrails,
"blending": {
"alpha_mu": round(best_alpha_mu, 6),
"alpha_sigma": round(best_alpha_sigma, 6),
},
"cities": city_params,
"metrics": {
"sample_count": len(normalized_samples),
"mean_crps": round(sum(crps_values) / len(crps_values), 6),
"legacy_mean_crps": round(legacy_mean_crps, 6),
"legacy_mean_mae": round(legacy_mean_mae, 6),
"legacy_bucket_hit_rate": round(legacy_bucket_hit_rate, 6),
"legacy_bucket_brier": round(legacy_bucket_brier, 6),
"selected_mean_crps": round(best_metrics["mean_crps"], 6),
"selected_mean_mae": round(best_metrics["mean_mae"], 6),
"selected_bucket_hit_rate": round(best_metrics["bucket_hit_rate"], 6),
"selected_bucket_brier": round(best_metrics["bucket_brier"], 6),
"selected_score": round(best_score, 6),
"legacy_score": round(legacy_score, 6),
},
}
def default_calibration_payload(
version: Optional[str] = None,
reason: str = "bootstrap",
) -> Dict[str, Any]:
return {
"version": version or "emos-bootstrap-v1",
"trained_at": datetime.now(timezone.utc).isoformat(),
"global": {
"mu": {
"intercept": 0.0,
"raw_mu_coef": 1.0,
"deb_coef": 0.0,
"ens_median_coef": 0.0,
"max_so_far_gap_coef": 0.0,
},
"sigma": {
"intercept": 0.0,
"raw_sigma_coef": 1.0,
"spread_coef": 0.0,
"peak_flag_coef": 0.0,
"max_so_far_gap_coef": 0.0,
},
},
"sigma_constraints": {
"min_ratio": 0.85,
"max_ratio": 1.35,
"absolute_min": 0.25,
"absolute_max": 3.0,
},
"selection_guardrails": {
"max_mae_increase": 0.02,
"max_bucket_hit_drop": 0.01,
"max_bucket_brier_increase": 0.05,
},
"blending": {
"alpha_mu": 1.0,
"alpha_sigma": 1.0,
},
"cities": {},
"metrics": {
"sample_count": 0,
"mean_crps": None,
"reason": reason,
},
}
def check_calibration_drift(records: list[dict[str, Any]], calibration_path: str | None = None) -> dict[str, Any]:
"""Compare recent CRPS against calibration baseline to detect drift.
Returns a dict with keys: drifted (bool), current_crps (float),
baseline_crps (float), delta_pct (float), sample_count (int), warning (str|None).
A positive delta_pct means the model is performing worse than baseline.
"""
import json
import os
if len(records) < 5:
return {"drifted": False, "sample_count": len(records), "warning": "Insufficient samples"}
path = calibration_path or DEFAULT_CALIBRATION_FILE
baseline_crps = None
try:
with open(path, "r", encoding="utf-8") as fh:
json.load(fh) # calibration file
eval_path = os.path.join(
os.path.dirname(path) if os.path.dirname(path) else os.path.join(os.path.dirname(DEFAULT_CALIBRATION_FILE)),
"evaluation_report.json",
)
with open(eval_path, "r", encoding="utf-8") as fh:
report = json.load(fh)
baseline_crps = float(report.get("metrics", {}).get("selected_mean_crps") or 0) or None
except Exception:
pass
current_crps_values: list[float] = []
for r in records:
actual = float(r.get("actual_high") or r.get("observed") or 0)
mu = float(r.get("deb_prediction") or r.get("mu") or 0)
sigma = float(r.get("sigma") or r.get("ensemble_std") or 2.0)
if actual and mu and sigma > 0:
current_crps_values.append(_gaussian_crps(actual, mu, sigma))
if not current_crps_values:
return {"drifted": False, "sample_count": 0, "warning": "No valid records for CRPS"}
current_crps = round(sum(current_crps_values) / len(current_crps_values), 6)
if baseline_crps is None or baseline_crps <= 0:
return {
"drifted": False,
"current_crps": current_crps,
"baseline_crps": baseline_crps,
"sample_count": len(records),
"warning": "No baseline available",
}
delta_pct = round((current_crps - baseline_crps) / baseline_crps * 100, 1)
DRIFT_THRESHOLD_PCT = 15.0 # warn if CRPS degraded by >15%
if delta_pct > DRIFT_THRESHOLD_PCT:
return {
"drifted": True,
"current_crps": current_crps,
"baseline_crps": baseline_crps,
"delta_pct": delta_pct,
"sample_count": len(records),
"warning": f"CRPS degraded {delta_pct}% vs baseline; consider re-running fit_calibration()",
}
return {
"drifted": False,
"current_crps": current_crps,
"baseline_crps": baseline_crps,
"delta_pct": delta_pct,
"sample_count": len(records),
}
-217
View File
@@ -1,217 +0,0 @@
from __future__ import annotations
import json
import os
from typing import Any, Dict, List, Optional
def _load_json_file(path: str) -> Dict[str, Any]:
try:
with open(path, "r", encoding="utf-8") as fh:
data = json.load(fh)
return data if isinstance(data, dict) else {}
except Exception:
return {}
def _sf(value: Any) -> Optional[float]:
if value is None:
return None
try:
return float(value)
except Exception:
return None
def _append_reason(reasons: List[str], condition: bool, message: str) -> None:
if condition:
reasons.append(message)
def _top_shadow_regressions(by_city: Dict[str, Any], limit: int = 5) -> List[Dict[str, Any]]:
rows: List[Dict[str, Any]] = []
for city, metrics in (by_city or {}).items():
if not isinstance(metrics, dict):
continue
rows.append(
{
"city": city,
"samples": int(metrics.get("samples") or 0),
"delta_mae": _sf(metrics.get("delta_mae")),
"delta_bucket_hit_rate": _sf(metrics.get("delta_bucket_hit_rate")),
"delta_bucket_brier": _sf(metrics.get("delta_bucket_brier")),
}
)
rows.sort(
key=lambda row: (
-(row["delta_bucket_brier"] or 0.0),
row["delta_bucket_hit_rate"] or 0.0,
-(row["delta_mae"] or 0.0),
)
)
return rows[:limit]
def judge_probability_rollout(
evaluation_report: Dict[str, Any],
shadow_report: Dict[str, Any],
) -> Dict[str, Any]:
thresholds = {
"evaluation_min_samples": 80,
"shadow_min_samples": 50,
"max_delta_mae": 0.05,
"min_delta_crps": -0.02,
"min_delta_bucket_hit_rate": 0.0,
"max_delta_bucket_brier_promote": 0.02,
"max_delta_bucket_brier_observe": 0.15,
}
eval_summary = (evaluation_report or {}).get("summary") or {}
eval_delta = eval_summary.get("delta") or {}
shadow_summary = (shadow_report or {}).get("summary") or {}
eval_samples = int(eval_summary.get("sample_count") or 0)
shadow_samples = int(shadow_summary.get("samples") or 0)
delta_crps = _sf(eval_delta.get("crps"))
delta_mae = _sf(eval_delta.get("mae"))
delta_hit = _sf(eval_delta.get("bucket_hit_rate"))
shadow_delta_mae = _sf(shadow_summary.get("delta_mae"))
shadow_delta_hit = _sf(shadow_summary.get("delta_bucket_hit_rate"))
shadow_delta_brier = _sf(shadow_summary.get("delta_bucket_brier"))
promote_reasons: List[str] = []
_append_reason(
promote_reasons,
eval_samples < thresholds["evaluation_min_samples"],
f"离线评估样本不足:{eval_samples} < {thresholds['evaluation_min_samples']}",
)
_append_reason(
promote_reasons,
shadow_samples < thresholds["shadow_min_samples"],
f"shadow 样本不足:{shadow_samples} < {thresholds['shadow_min_samples']}",
)
_append_reason(
promote_reasons,
delta_crps is None or delta_crps > thresholds["min_delta_crps"],
f"离线 CRPS 改善不足:delta={delta_crps}",
)
_append_reason(
promote_reasons,
delta_mae is None or delta_mae > thresholds["max_delta_mae"],
f"离线 MAE 退化超限:delta={delta_mae}",
)
_append_reason(
promote_reasons,
delta_hit is None or delta_hit < thresholds["min_delta_bucket_hit_rate"],
f"离线 bucket 命中率下降:delta={delta_hit}",
)
_append_reason(
promote_reasons,
shadow_delta_mae is None or shadow_delta_mae > thresholds["max_delta_mae"],
f"shadow MAE 退化超限:delta={shadow_delta_mae}",
)
_append_reason(
promote_reasons,
shadow_delta_hit is None or shadow_delta_hit < thresholds["min_delta_bucket_hit_rate"],
f"shadow bucket 命中率下降:delta={shadow_delta_hit}",
)
_append_reason(
promote_reasons,
shadow_delta_brier is None
or shadow_delta_brier > thresholds["max_delta_bucket_brier_promote"],
f"shadow bucket brier 退化超限:delta={shadow_delta_brier}",
)
if not promote_reasons:
decision = "promote"
summary = "离线与 shadow 指标均达标,可以考虑切换 emos_primary。"
else:
observe_reasons: List[str] = []
_append_reason(
observe_reasons,
eval_samples < thresholds["evaluation_min_samples"],
f"离线评估样本不足:{eval_samples}",
)
_append_reason(
observe_reasons,
delta_crps is None or delta_crps > thresholds["min_delta_crps"],
f"离线 CRPS 改善不足:delta={delta_crps}",
)
_append_reason(
observe_reasons,
delta_mae is None or delta_mae > thresholds["max_delta_mae"],
f"离线 MAE 退化超限:delta={delta_mae}",
)
_append_reason(
observe_reasons,
delta_hit is None or delta_hit < thresholds["min_delta_bucket_hit_rate"],
f"离线 bucket 命中率下降:delta={delta_hit}",
)
_append_reason(
observe_reasons,
shadow_samples < thresholds["shadow_min_samples"],
f"shadow 样本不足:{shadow_samples}",
)
_append_reason(
observe_reasons,
shadow_delta_mae is None or shadow_delta_mae > thresholds["max_delta_mae"],
f"shadow MAE 退化超限:delta={shadow_delta_mae}",
)
_append_reason(
observe_reasons,
shadow_delta_hit is None or shadow_delta_hit < thresholds["min_delta_bucket_hit_rate"],
f"shadow bucket 命中率下降:delta={shadow_delta_hit}",
)
_append_reason(
observe_reasons,
shadow_delta_brier is None
or shadow_delta_brier > thresholds["max_delta_bucket_brier_observe"],
f"shadow bucket brier 退化偏大:delta={shadow_delta_brier}",
)
if not observe_reasons:
decision = "observe"
summary = "离线评估达标,但 shadow 仍需继续观察,暂不切主路径。"
else:
decision = "hold"
summary = "当前指标不足以切换 emos_primary,应继续保持 shadow。"
return {
"decision": decision,
"ready_for_primary": decision == "promote",
"summary": summary,
"thresholds": thresholds,
"evaluation": {
"sample_count": eval_samples,
"delta_crps": delta_crps,
"delta_mae": delta_mae,
"delta_bucket_hit_rate": delta_hit,
},
"shadow": {
"sample_count": shadow_samples,
"delta_mae": shadow_delta_mae,
"delta_bucket_hit_rate": shadow_delta_hit,
"delta_bucket_brier": shadow_delta_brier,
},
"blocking_reasons": promote_reasons,
"worst_shadow_regressions": _top_shadow_regressions(
(shadow_report or {}).get("by_city") or {}
),
}
def build_rollout_report(
evaluation_report_path: str,
shadow_report_path: str,
) -> Dict[str, Any]:
evaluation_report = _load_json_file(evaluation_report_path)
shadow_report = _load_json_file(shadow_report_path)
decision = judge_probability_rollout(evaluation_report, shadow_report)
return {
"evaluation_report_path": evaluation_report_path,
"shadow_report_path": shadow_report_path,
"evaluation_report_exists": os.path.exists(evaluation_report_path),
"shadow_report_exists": os.path.exists(shadow_report_path),
"decision": decision,
}
@@ -1,297 +0,0 @@
from __future__ import annotations
import json
import os
from datetime import datetime
from typing import Any, Dict, List, Optional
from src.database.runtime_state import (
ProbabilitySnapshotRepository,
STATE_STORAGE_SQLITE,
TrainingFeatureRecordRepository,
get_state_storage_mode,
)
DEDUP_SCAN_LINES = 200
MU_THRESHOLD = 0.2
SIGMA_THRESHOLD = 0.15
MAX_SO_FAR_THRESHOLD = 0.2
_snapshot_repo = ProbabilitySnapshotRepository()
_training_feature_repo = TrainingFeatureRecordRepository()
def _sf(value: Any) -> Optional[float]:
if value is None:
return None
try:
return float(value)
except Exception:
return None
def _compact_snapshot(distribution: Optional[List[Dict[str, Any]]]) -> List[Dict[str, Any]]:
compact: List[Dict[str, Any]] = []
for row in distribution or []:
if not isinstance(row, dict):
continue
value = row.get("value")
probability = row.get("probability")
if value is None or probability is None:
continue
try:
compact.append(
{
"v": int(value),
"p": round(float(probability), 3),
}
)
except Exception:
continue
if len(compact) >= 4:
break
return compact
def _top_bucket(snapshot: Optional[List[Dict[str, Any]]]) -> Optional[int]:
best_value = None
best_prob = -1.0
for row in snapshot or []:
if not isinstance(row, dict):
continue
value = row.get("v")
prob = _sf(row.get("p"))
if value is None or prob is None:
continue
if prob > best_prob:
best_value = int(value)
best_prob = prob
return best_value
def _load_recent_rows(path: str, max_lines: int = DEDUP_SCAN_LINES) -> List[Dict[str, Any]]:
if not os.path.exists(path):
return []
with open(path, "r", encoding="utf-8") as fh:
lines = fh.readlines()[-max_lines:]
rows = []
for line in lines:
line = line.strip()
if not line:
continue
try:
row = json.loads(line)
except Exception:
continue
if isinstance(row, dict):
rows.append(row)
return rows
def load_snapshot_rows_for_day(
city_name: str,
target_date: str,
archive_path: Optional[str] = None,
) -> List[Dict[str, Any]]:
city_key = str(city_name or "").strip().lower()
date_key = str(target_date or "").strip()
if not city_key or not date_key:
return []
mode = get_state_storage_mode()
if mode == STATE_STORAGE_SQLITE:
return _snapshot_repo.load_rows_by_city_date(city_key, date_key)
root_dir = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
path = archive_path or os.path.join(root_dir, "data", "probability_training_snapshots.jsonl")
if not os.path.exists(path):
if mode == STATE_STORAGE_SQLITE:
return _snapshot_repo.load_rows_by_city_date(city_key, date_key)
return []
rows: List[Dict[str, Any]] = []
try:
with open(path, "r", encoding="utf-8") as fh:
for line in fh:
line = line.strip()
if not line:
continue
try:
row = json.loads(line)
except Exception:
continue
if not isinstance(row, dict):
continue
if str(row.get("city") or "").strip().lower() != city_key:
continue
if str(row.get("date") or "").strip() != date_key:
continue
rows.append(row)
except Exception:
return []
rows.sort(key=lambda row: str(row.get("timestamp") or ""))
return rows
def _should_skip_append(path: str, payload: Dict[str, Any]) -> bool:
mode = get_state_storage_mode()
if mode == STATE_STORAGE_SQLITE:
recent_rows = _snapshot_repo.load_recent_rows(
str(payload.get("city") or ""),
str(payload.get("date") or ""),
DEDUP_SCAN_LINES,
)
else:
recent_rows = _load_recent_rows(path)
city = payload.get("city")
date_str = payload.get("date")
if not city or not date_str:
return False
for row in reversed(recent_rows):
if row.get("city") != city or row.get("date") != date_str:
continue
if row.get("peak_status") != payload.get("peak_status"):
return False
if row.get("probability_mode") != payload.get("probability_mode"):
return False
current_top = _top_bucket(payload.get("prob_snapshot"))
previous_top = _top_bucket(row.get("prob_snapshot"))
current_shadow_top = _top_bucket(payload.get("shadow_prob_snapshot"))
previous_shadow_top = _top_bucket(row.get("shadow_prob_snapshot"))
if current_top != previous_top or current_shadow_top != previous_shadow_top:
return False
if abs((_sf(payload.get("raw_mu")) or 0.0) - (_sf(row.get("raw_mu")) or 0.0)) > MU_THRESHOLD:
return False
if abs((_sf(payload.get("raw_sigma")) or 0.0) - (_sf(row.get("raw_sigma")) or 0.0)) > SIGMA_THRESHOLD:
return False
if abs((_sf(payload.get("max_so_far")) or 0.0) - (_sf(row.get("max_so_far")) or 0.0)) > MAX_SO_FAR_THRESHOLD:
return False
return True
return False
def append_probability_snapshot(
city_name: str,
*,
local_date: str,
observation_time: Optional[str],
temp_symbol: str,
raw_mu: Optional[float],
raw_sigma: Optional[float],
deb_prediction: Optional[float],
ens_data: Optional[Dict[str, Any]],
current_forecasts: Optional[Dict[str, Any]],
max_so_far: Optional[float],
current_temp: Optional[float] = None,
humidity: Optional[float] = None,
wind_speed_kt: Optional[float] = None,
visibility_mi: Optional[float] = None,
local_hour: Optional[float] = None,
peak_status: Optional[str],
probabilities: Optional[List[Dict[str, Any]]],
shadow_probabilities: Optional[List[Dict[str, Any]]],
calibration_summary: Optional[Dict[str, Any]],
archive_path: Optional[str] = None,
) -> None:
city_key = str(city_name or "").strip().lower()
if not city_key:
return
root_dir = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
path = archive_path or os.path.join(
root_dir,
"data",
"probability_training_snapshots.jsonl",
)
calibration_summary = calibration_summary or {}
ens_data = ens_data or {}
current_forecasts = current_forecasts or {}
timestamp = str(observation_time or datetime.utcnow().isoformat() + "Z").strip()
payload = {
"city": city_key,
"timestamp": timestamp,
"date": local_date,
"temp_symbol": temp_symbol,
"raw_mu": _sf(raw_mu),
"raw_sigma": _sf(raw_sigma),
"deb_prediction": _sf(deb_prediction),
"ensemble": {
"p10": _sf(ens_data.get("p10")),
"median": _sf(ens_data.get("median")),
"p90": _sf(ens_data.get("p90")),
},
"multi_model": {
key: _sf(value)
for key, value in current_forecasts.items()
if _sf(value) is not None
},
"max_so_far": _sf(max_so_far),
"observation": {
"current_temp": _sf(current_temp),
"humidity": _sf(humidity),
"wind_speed_kt": _sf(wind_speed_kt),
"visibility_mi": _sf(visibility_mi),
"local_hour": _sf(local_hour),
},
"peak_status": peak_status,
"prob_snapshot": _compact_snapshot(probabilities),
"shadow_prob_snapshot": _compact_snapshot(shadow_probabilities),
"probability_engine": calibration_summary.get("engine"),
"probability_mode": calibration_summary.get("mode"),
"calibration_version": calibration_summary.get("calibration_version"),
"calibration_source": calibration_summary.get("calibration_source"),
"calibrated_mu": _sf(calibration_summary.get("calibrated_mu")),
"calibrated_sigma": _sf(calibration_summary.get("calibrated_sigma")),
}
parent = os.path.dirname(os.path.abspath(path))
if parent:
os.makedirs(parent, exist_ok=True)
if _should_skip_append(path, payload):
return
mode = get_state_storage_mode()
if mode == STATE_STORAGE_SQLITE:
_snapshot_repo.append_snapshot(payload)
_training_feature_repo.upsert_record(
city_key,
local_date,
{
"forecasts": payload.get("multi_model") or {},
"deb_prediction": payload.get("deb_prediction"),
"mu": payload.get("raw_mu"),
"probability_features": {
"raw_mu": payload.get("raw_mu"),
"raw_sigma": payload.get("raw_sigma"),
"deb_prediction": payload.get("deb_prediction"),
"ens_median": (payload.get("ensemble") or {}).get("median"),
"ensemble_spread": None,
"max_so_far": payload.get("max_so_far"),
"peak_status": payload.get("peak_status"),
},
"prob_snapshot": payload.get("prob_snapshot") or [],
"shadow_prob_snapshot": payload.get("shadow_prob_snapshot") or [],
"probability_calibration": {
"engine": payload.get("probability_engine"),
"mode": payload.get("probability_mode"),
"calibration_version": payload.get("calibration_version"),
"calibration_source": payload.get("calibration_source"),
"calibrated_mu": payload.get("calibrated_mu"),
"calibrated_sigma": payload.get("calibrated_sigma"),
},
"observation": payload.get("observation") or {},
"snapshot_timestamp": payload.get("timestamp"),
},
)
if mode != STATE_STORAGE_SQLITE:
with open(path, "a", encoding="utf-8") as fh:
fh.write(json.dumps(payload, ensure_ascii=False) + "\n")
+2 -103
View File
@@ -15,11 +15,6 @@ from src.analysis.deb_algorithm import (
update_daily_record,
_is_excluded_model_name,
)
from src.analysis.probability_calibration import (
apply_probability_calibration,
build_probability_features,
)
from src.analysis.probability_snapshot_archive import append_probability_snapshot
from src.analysis.settlement_rounding import apply_city_settlement, is_exact_settlement_city
from src.data_collection.city_registry import CITY_REGISTRY
from src.data_collection.city_risk_profiles import get_city_risk_profile
@@ -595,20 +590,7 @@ def analyze_weather_trend(
# === Probability Engine ===
probabilities: List[Dict[str, Any]] = []
probabilities_all: List[Dict[str, Any]] = []
shadow_probabilities: List[Dict[str, Any]] = []
shadow_probabilities_all: List[Dict[str, Any]] = []
forecast_miss_deg = 0.0
probability_features = None
calibration_summary = {
"mode": "legacy",
"engine": "legacy",
"raw_mu": None,
"raw_sigma": sigma,
"calibrated_mu": None,
"calibrated_sigma": None,
"calibration_version": None,
"calibration_source": None,
}
if is_dead_market:
settled_wu = apply_city_settlement(city_name, max_so_far) if max_so_far is not None else 0
@@ -662,7 +644,7 @@ def analyze_weather_trend(
f"实测最高 {max_so_far}{temp_symbol},偏差 {forecast_miss_deg}°。当前趋势: {_trend_dir}"
)
# Probability Engine
# Probability (legacy Gaussian buckets)
probs_result = calculate_prob_distribution(
mu, sigma, max_so_far, temp_symbol, city_name
)
@@ -671,45 +653,6 @@ def analyze_weather_trend(
probabilities_all = probs_result.get("probabilities_all", probabilities)
sorted_probs = probs_result.get("sorted_probs", [])
probability_features = build_probability_features(
city_name=city_name or "",
raw_mu=mu,
raw_sigma=sigma,
deb_prediction=deb_prediction,
ens_data=ens_data,
current_forecasts=current_forecasts,
max_so_far=max_so_far,
peak_status=peak_status,
local_hour_frac=local_hour_frac,
)
calibration_result = apply_probability_calibration(
city_name=city_name or "",
temp_symbol=temp_symbol,
raw_mu=mu,
raw_sigma=sigma,
max_so_far=max_so_far,
legacy_distribution=probabilities,
features=probability_features,
)
calibration_summary = {
"mode": calibration_result.get("mode", "legacy"),
"engine": calibration_result.get("engine", "legacy"),
"raw_mu": calibration_result.get("raw_mu"),
"raw_sigma": calibration_result.get("raw_sigma"),
"calibrated_mu": calibration_result.get("calibrated_mu"),
"calibrated_sigma": calibration_result.get("calibrated_sigma"),
"calibration_version": calibration_result.get("calibration_version"),
"calibration_source": calibration_result.get("calibration_source"),
}
shadow_probabilities = calibration_result.get("shadow_distribution") or []
shadow_probabilities_all = calibration_result.get("shadow_distribution_all") or shadow_probabilities
if calibration_result.get("engine") == "emos":
mu = calibration_result.get("calibrated_mu", mu)
sigma = calibration_result.get("calibrated_sigma", sigma)
probabilities = calibration_result.get("distribution") or probabilities
probabilities_all = calibration_result.get("distribution_all") or probabilities_all or probabilities
sorted_probs = calibration_result.get("selected_sorted_probs") or sorted_probs
if sorted_probs:
prob_parts = [
f"{int(t)}{temp_symbol} [{t - 0.5}~{t + 0.5}) {p * 100:.0f}%"
@@ -868,7 +811,6 @@ def analyze_weather_trend(
# === Save daily record (with μ + prob snapshot) ===
try:
_prob_list = None
_shadow_prob_list = None
if sorted_probs:
_prob_list = [
{"value": int(t), "probability": round(p, 3)}
@@ -876,12 +818,6 @@ def analyze_weather_trend(
]
elif is_dead_market and max_so_far is not None:
_prob_list = [{"value": apply_city_settlement(city_name, max_so_far), "probability": 1.0}]
if shadow_probabilities:
_shadow_prob_list = [
{"value": int(row.get("value")), "probability": round(float(row.get("probability") or 0.0), 3)}
for row in shadow_probabilities[:4]
if row.get("value") is not None
]
update_daily_record(
city_name,
@@ -891,34 +827,6 @@ def analyze_weather_trend(
deb_prediction=_deb_to_save,
mu=mu,
probabilities=_prob_list,
probability_features=probability_features,
shadow_probabilities=_shadow_prob_list,
calibration_summary=calibration_summary,
)
except Exception:
pass
try:
append_probability_snapshot(
city_name=city_name or "",
local_date=local_date_str,
observation_time=obs_time_raw or local_time_full or None,
temp_symbol=temp_symbol,
raw_mu=calibration_summary.get("raw_mu"),
raw_sigma=calibration_summary.get("raw_sigma"),
deb_prediction=_deb_to_save,
ens_data=ens_data,
current_forecasts=current_forecasts,
max_so_far=max_so_far,
current_temp=cur_temp,
humidity=_sf(primary_current.get("humidity")),
wind_speed_kt=_sf(primary_current.get("wind_speed_kt")),
visibility_mi=_sf(primary_current.get("visibility_mi")),
local_hour=local_hour_frac,
peak_status=peak_status,
probabilities=_prob_list,
shadow_probabilities=_shadow_prob_list,
calibration_summary=calibration_summary,
)
except Exception:
pass
@@ -933,16 +841,7 @@ def analyze_weather_trend(
"mu": mu,
"probabilities": probabilities,
"probabilities_all": probabilities_all or probabilities,
"shadow_probabilities": shadow_probabilities,
"shadow_probabilities_all": shadow_probabilities_all or shadow_probabilities,
"probability_engine": calibration_summary["engine"],
"probability_calibration_mode": calibration_summary["mode"],
"probability_calibration_version": calibration_summary["calibration_version"],
"probability_calibration_source": calibration_summary["calibration_source"],
"probability_raw_mu": calibration_summary["raw_mu"],
"probability_raw_sigma": calibration_summary["raw_sigma"],
"probability_calibrated_mu": calibration_summary["calibrated_mu"],
"probability_calibrated_sigma": calibration_summary["calibrated_sigma"],
"probability_engine": "legacy",
"trend_info": {
"direction": trend_direction if 'trend_direction' in dir() else "unknown",
"recent": recent_list,
-171
View File
@@ -1,171 +0,0 @@
from __future__ import annotations
import json
import os
from typing import Any, Dict, List, Optional, Tuple
from loguru import logger
from src.models.lgbm_features import (
FEATURE_NAMES,
build_runtime_feature_map,
)
_MODEL_CACHE: Dict[str, Any] = {"path": None, "mtime": None, "booster": None}
_SCHEMA_CACHE: Dict[str, Any] = {"path": None, "mtime": None, "schema": None}
def _sf(value: Any) -> Optional[float]:
if value is None:
return None
try:
return float(value)
except Exception:
return None
def _truthy_env(name: str, default: str = "false") -> bool:
return str(os.getenv(name, default)).strip().lower() in {"1", "true", "yes", "on"}
def lgbm_model_path() -> str:
root = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
return str(
os.getenv(
"POLYWEATHER_LGBM_MODEL_PATH",
os.path.join(root, "artifacts", "models", "lgbm_daily_high.txt"),
)
).strip()
def lgbm_schema_path() -> str:
root = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
return str(
os.getenv(
"POLYWEATHER_LGBM_SCHEMA_PATH",
os.path.join(root, "artifacts", "models", "lgbm_daily_high_schema.json"),
)
).strip()
def lgbm_min_history_points() -> int:
try:
return max(1, int(os.getenv("POLYWEATHER_LGBM_MIN_HISTORY_POINTS", "3")))
except Exception:
return 3
def is_lgbm_enabled() -> bool:
return _truthy_env("POLYWEATHER_LGBM_ENABLED", "false")
def _load_schema(schema_path: str) -> Optional[Dict[str, Any]]:
if not schema_path or not os.path.exists(schema_path):
return None
mtime = os.path.getmtime(schema_path)
if (
_SCHEMA_CACHE["schema"] is not None
and _SCHEMA_CACHE["path"] == schema_path
and _SCHEMA_CACHE["mtime"] == mtime
):
return _SCHEMA_CACHE["schema"]
with open(schema_path, "r", encoding="utf-8") as fh:
data = json.load(fh)
if not isinstance(data, dict):
return None
_SCHEMA_CACHE.update({"path": schema_path, "mtime": mtime, "schema": data})
return data
def _load_booster(model_path: str):
if not model_path or not os.path.exists(model_path):
return None
mtime = os.path.getmtime(model_path)
if (
_MODEL_CACHE["booster"] is not None
and _MODEL_CACHE["path"] == model_path
and _MODEL_CACHE["mtime"] == mtime
):
return _MODEL_CACHE["booster"]
try:
import lightgbm as lgb
except Exception as exc:
logger.warning(f"LGBM runtime dependency missing: {exc}")
return None
booster = lgb.Booster(model_file=model_path)
_MODEL_CACHE.update({"path": model_path, "mtime": mtime, "booster": booster})
return booster
def _vector_from_features(
feature_map: Dict[str, Optional[float]],
schema: Optional[Dict[str, Any]],
) -> List[float]:
feature_names = schema.get("feature_names") if isinstance(schema, dict) else None
ordered_names = feature_names if isinstance(feature_names, list) and feature_names else FEATURE_NAMES
vector: List[float] = []
for name in ordered_names:
value = feature_map.get(str(name))
vector.append(float("nan") if value is None else float(value))
return vector
def predict_lgbm_daily_high(
*,
city_name: str,
current_forecasts: Dict[str, Any],
deb_prediction: Optional[float],
current_temp: Optional[float],
max_so_far: Optional[float],
humidity: Optional[float],
wind_speed_kt: Optional[float],
visibility_mi: Optional[float],
local_hour: int,
local_date: str,
peak_status: str,
history_data: Optional[Dict[str, Any]] = None,
) -> Tuple[Optional[float], Dict[str, Any]]:
if not is_lgbm_enabled():
return None, {"reason": "disabled"}
schema = _load_schema(lgbm_schema_path())
booster = _load_booster(lgbm_model_path())
if schema is None or booster is None:
return None, {"reason": "artifact_missing"}
feature_map, meta = build_runtime_feature_map(
city_name=city_name,
current_forecasts=current_forecasts,
deb_prediction=deb_prediction,
current_temp=current_temp,
max_so_far=max_so_far,
humidity=humidity,
wind_speed_kt=wind_speed_kt,
visibility_mi=visibility_mi,
local_hour=local_hour,
local_date=local_date,
peak_status=peak_status,
history_data=history_data,
)
if not feature_map:
return None, meta
if int(meta.get("history_count") or 0) < lgbm_min_history_points():
return None, {
"reason": "insufficient_history",
"history_count": int(meta.get("history_count") or 0),
}
try:
vector = _vector_from_features(feature_map, schema)
prediction = booster.predict([vector], num_iteration=booster.best_iteration)
value = _sf(prediction[0] if prediction is not None else None)
if value is None:
return None, {"reason": "empty_prediction"}
return round(float(value), 1), {
"reason": "ok",
"history_count": int(meta.get("history_count") or 0),
}
except Exception as exc:
logger.warning(f"LGBM prediction failed for {city_name}: {exc}")
return None, {"reason": "predict_failed", "error": str(exc)}
-472
View File
@@ -1,472 +0,0 @@
from __future__ import annotations
import json
import os
from datetime import datetime
from statistics import mean
from typing import Any, Dict, List, Optional, Tuple
from src.analysis.deb_algorithm import load_history
from src.data_collection.city_registry import ALIASES
from src.database.runtime_state import (
DailyRecordRepository,
ProbabilitySnapshotRepository,
STATE_STORAGE_FILE,
STATE_STORAGE_SQLITE,
TrainingFeatureRecordRepository,
TruthRecordRepository,
get_state_storage_mode,
)
BASE_MODEL_COLUMNS: List[Tuple[str, str]] = [
("Open-Meteo", "open_meteo"),
("ECMWF", "ecmwf"),
("GFS", "gfs"),
("GEM", "gem"),
("JMA", "jma"),
("ICON", "icon"),
("MGM", "mgm"),
("NWS", "nws"),
]
FEATURE_NAMES: List[str] = [
"actual_high_lag_1",
"actual_high_lag_2",
"actual_high_lag_3",
"actual_high_lag_7",
"actual_high_mean_7",
"actual_high_mean_14",
"actual_high_trend_3",
*[column for _, column in BASE_MODEL_COLUMNS],
"deb_prediction",
"model_median",
"model_spread",
"current_temp",
"max_so_far",
"humidity",
"wind_speed_kt",
"visibility_mi",
"local_hour",
"month",
"weekday",
"peak_status_code",
]
PEAK_STATUS_CODES = {
"before": 0.0,
"in_window": 1.0,
"past": 2.0,
}
def _sf(value: Any) -> Optional[float]:
if value is None:
return None
try:
return float(value)
except Exception:
return None
def _parse_date(value: Any) -> Optional[datetime]:
text = str(value or "").strip()
if not text:
return None
try:
return datetime.strptime(text, "%Y-%m-%d")
except Exception:
return None
def _parse_timestamp(value: Any) -> Optional[datetime]:
text = str(value or "").strip()
if not text:
return None
try:
return datetime.fromisoformat(text.replace("Z", "+00:00"))
except Exception:
return None
def _history_file_path() -> str:
root = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
if get_state_storage_mode() == STATE_STORAGE_FILE:
return os.path.join(root, "data", "daily_records.json")
return ""
def _snapshot_archive_path() -> str:
root = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
if get_state_storage_mode() == STATE_STORAGE_FILE:
return os.path.join(root, "data", "probability_training_snapshots.jsonl")
return ""
def _normalized_city_key(city_name: str) -> str:
city_key = str(city_name or "").strip().lower()
return ALIASES.get(city_key, city_key)
def _safe_mean(values: List[Optional[float]]) -> Optional[float]:
valid = [float(v) for v in values if v is not None]
if not valid:
return None
return float(mean(valid))
def _peak_status_code(value: Any) -> Optional[float]:
status = str(value or "").strip().lower()
if not status:
return None
return PEAK_STATUS_CODES.get(status, -1.0)
def _compute_model_summary(features: Dict[str, Optional[float]]) -> Tuple[Optional[float], Optional[float]]:
values = [
features.get(column)
for _, column in BASE_MODEL_COLUMNS
if features.get(column) is not None
]
values = [float(v) for v in values if v is not None]
if not values:
return None, None
ordered = sorted(values)
median_value = ordered[len(ordered) // 2]
spread_value = ordered[-1] - ordered[0] if len(ordered) >= 2 else 0.0
return float(median_value), float(spread_value)
def load_snapshot_index(archive_path: Optional[str] = None) -> Dict[Tuple[str, str], Dict[str, Any]]:
mode = get_state_storage_mode()
if mode == STATE_STORAGE_SQLITE:
latest_rows: Dict[Tuple[str, str], Dict[str, Any]] = {}
for row in ProbabilitySnapshotRepository().load_all_rows():
if not isinstance(row, dict):
continue
city = _normalized_city_key(str(row.get("city") or ""))
date_str = str(row.get("date") or "").strip()
if not city or not date_str:
continue
key = (city, date_str)
current_best = latest_rows.get(key)
if current_best is None or str(row.get("timestamp") or "") >= str(
current_best.get("timestamp") or ""
):
latest_rows[key] = row
return latest_rows
path = archive_path or _snapshot_archive_path()
if not path or not os.path.exists(path):
return {}
latest_rows: Dict[Tuple[str, str], Dict[str, Any]] = {}
with open(path, "r", encoding="utf-8") as fh:
for line in fh:
line = line.strip()
if not line:
continue
try:
row = json.loads(line)
except Exception:
continue
if not isinstance(row, dict):
continue
city = _normalized_city_key(str(row.get("city") or ""))
date_str = str(row.get("date") or "").strip()
if not city or not date_str:
continue
key = (city, date_str)
current_best = latest_rows.get(key)
if current_best is None or str(row.get("timestamp") or "") >= str(
current_best.get("timestamp") or ""
):
latest_rows[key] = row
return latest_rows
def _extract_history_rows(
history_data: Dict[str, Any],
city_name: str,
exclude_date: Optional[str] = None,
) -> List[Tuple[str, float]]:
city_key = _normalized_city_key(city_name)
city_rows = history_data.get(city_key) if isinstance(history_data, dict) else None
if not isinstance(city_rows, dict):
return []
rows: List[Tuple[str, float]] = []
for date_str, record in city_rows.items():
if exclude_date and str(date_str) >= str(exclude_date):
continue
if not isinstance(record, dict):
continue
actual = _sf(record.get("actual_high"))
if actual is None:
continue
rows.append((str(date_str), float(actual)))
rows.sort(key=lambda item: item[0])
return rows
def _lag(values: List[float], distance: int) -> Optional[float]:
if len(values) < distance:
return None
return float(values[-distance])
def build_runtime_feature_map(
*,
city_name: str,
current_forecasts: Dict[str, Any],
deb_prediction: Optional[float],
current_temp: Optional[float],
max_so_far: Optional[float],
humidity: Optional[float],
wind_speed_kt: Optional[float],
visibility_mi: Optional[float],
local_hour: int,
local_date: str,
peak_status: str,
history_data: Optional[Dict[str, Any]] = None,
) -> Tuple[Optional[Dict[str, Optional[float]]], Dict[str, Any]]:
data = history_data if isinstance(history_data, dict) else load_history(_history_file_path())
history_rows = _extract_history_rows(data, city_name, exclude_date=local_date)
history_values = [value for _, value in history_rows]
if not history_values:
return None, {"reason": "no_history", "history_count": 0}
date_obj = _parse_date(local_date)
if date_obj is None:
return None, {"reason": "invalid_date", "history_count": len(history_values)}
features: Dict[str, Optional[float]] = {
"actual_high_lag_1": _lag(history_values, 1),
"actual_high_lag_2": _lag(history_values, 2),
"actual_high_lag_3": _lag(history_values, 3),
"actual_high_lag_7": _lag(history_values, 7),
"actual_high_mean_7": _safe_mean(history_values[-7:]),
"actual_high_mean_14": _safe_mean(history_values[-14:]),
"actual_high_trend_3": (
history_values[-1] - history_values[-3] if len(history_values) >= 3 else None
),
"deb_prediction": _sf(deb_prediction),
"current_temp": _sf(current_temp),
"max_so_far": _sf(max_so_far),
"humidity": _sf(humidity),
"wind_speed_kt": _sf(wind_speed_kt),
"visibility_mi": _sf(visibility_mi),
"local_hour": float(local_hour),
"month": float(date_obj.month),
"weekday": float(date_obj.weekday()),
"peak_status_code": _peak_status_code(peak_status),
}
for model_name, column in BASE_MODEL_COLUMNS:
features[column] = _sf(current_forecasts.get(model_name))
model_median, model_spread = _compute_model_summary(features)
features["model_median"] = model_median
features["model_spread"] = model_spread
return features, {
"reason": "ok",
"history_count": len(history_values),
}
def _features_to_vector(features: Dict[str, Optional[float]], feature_names: Optional[List[str]] = None) -> List[float]:
ordered_names = feature_names or FEATURE_NAMES
vector: List[float] = []
for name in ordered_names:
value = features.get(name)
vector.append(float("nan") if value is None else float(value))
return vector
def build_training_samples(
history_data: Optional[Dict[str, Any]] = None,
snapshot_index: Optional[Dict[Tuple[str, str], Dict[str, Any]]] = None,
) -> List[Dict[str, Any]]:
mode = get_state_storage_mode()
if isinstance(history_data, dict):
runtime_history = history_data
elif mode == STATE_STORAGE_SQLITE:
runtime_history = DailyRecordRepository().load_all()
else:
runtime_history = load_history(_history_file_path())
if isinstance(history_data, dict):
truth_history = runtime_history
training_feature_history = {}
elif mode == STATE_STORAGE_SQLITE:
truth_history = TruthRecordRepository().load_all()
training_feature_history = TrainingFeatureRecordRepository().load_all()
else:
truth_history = runtime_history
training_feature_history = {}
snapshots = snapshot_index if isinstance(snapshot_index, dict) else load_snapshot_index()
samples: List[Dict[str, Any]] = []
excluded_keys: set[tuple[str, str]] = set()
for (city_name, date_str), snapshot in (snapshots or {}).items():
if not isinstance(snapshot, dict):
continue
truth_row = ((truth_history.get(city_name) or {}).get(str(date_str)) or {})
target = _sf(truth_row.get("actual_high"))
if target is None:
runtime_record = ((runtime_history.get(city_name) or {}).get(str(date_str)) or {})
target = _sf(runtime_record.get("actual_high"))
if target is None:
continue
observation = snapshot.get("observation") if isinstance(snapshot.get("observation"), dict) else {}
current_forecasts = snapshot.get("multi_model") if isinstance(snapshot.get("multi_model"), dict) else {}
local_hour = _sf(observation.get("local_hour"))
if local_hour is None:
timestamp = _parse_timestamp(snapshot.get("timestamp"))
local_hour = float(timestamp.hour) if timestamp is not None else 12.0
feature_map, meta = build_runtime_feature_map(
city_name=city_name,
current_forecasts=current_forecasts,
deb_prediction=_sf(snapshot.get("deb_prediction")) or _sf(snapshot.get("raw_mu")),
current_temp=_sf(observation.get("current_temp")),
max_so_far=_sf(snapshot.get("max_so_far")),
humidity=_sf(observation.get("humidity")),
wind_speed_kt=_sf(observation.get("wind_speed_kt")),
visibility_mi=_sf(observation.get("visibility_mi")),
local_hour=int(local_hour),
local_date=str(date_str),
peak_status=str(snapshot.get("peak_status") or "before"),
history_data=truth_history,
)
if not feature_map:
continue
samples.append(
{
"city": _normalized_city_key(city_name),
"date": str(date_str),
"target": float(target),
"features": feature_map,
"vector": _features_to_vector(feature_map),
"history_count": int(meta.get("history_count") or 0),
"deb_prediction": _sf(snapshot.get("deb_prediction")) or _sf(snapshot.get("raw_mu")),
"forecasts": {
key: _sf(value)
for key, value in current_forecasts.items()
if _sf(value) is not None
},
"sample_source": "snapshot",
"settlement_source": truth_row.get("settlement_source"),
"settlement_station_code": truth_row.get("settlement_station_code"),
"truth_version": truth_row.get("truth_version"),
"truth_updated_by": truth_row.get("updated_by"),
"truth_updated_at": truth_row.get("truth_updated_at"),
}
)
excluded_keys.add((_normalized_city_key(city_name), str(date_str)))
training_source = training_feature_history or runtime_history or {}
for city_name, city_records in training_source.items():
if not isinstance(city_records, dict):
continue
ordered_dates = sorted(city_records.keys())
for date_str in ordered_dates:
normalized_city = _normalized_city_key(city_name)
if (normalized_city, str(date_str)) in excluded_keys:
continue
record = city_records.get(date_str)
if not isinstance(record, dict):
continue
truth_row = ((truth_history.get(normalized_city) or {}).get(str(date_str)) or {})
target = _sf(truth_row.get("actual_high"))
if target is None:
target = _sf(((runtime_history.get(normalized_city) or {}).get(str(date_str)) or {}).get("actual_high"))
forecasts = record.get("forecasts") if isinstance(record.get("forecasts"), dict) else {}
if target is None or not forecasts:
continue
feature_map, meta = build_runtime_feature_map(
city_name=city_name,
current_forecasts=forecasts,
deb_prediction=_sf(record.get("deb_prediction")),
current_temp=None,
max_so_far=None,
humidity=None,
wind_speed_kt=None,
visibility_mi=None,
local_hour=12,
local_date=str(date_str),
peak_status="before",
history_data=truth_history,
)
if not feature_map:
continue
snapshot = snapshots.get((_normalized_city_key(city_name), str(date_str))) or {}
observation = snapshot.get("observation") if isinstance(snapshot.get("observation"), dict) else {}
feature_map["max_so_far"] = _sf(snapshot.get("max_so_far"))
feature_map["current_temp"] = _sf(observation.get("current_temp"))
feature_map["humidity"] = _sf(observation.get("humidity"))
feature_map["wind_speed_kt"] = _sf(observation.get("wind_speed_kt"))
feature_map["visibility_mi"] = _sf(observation.get("visibility_mi"))
local_hour = _sf(observation.get("local_hour"))
if local_hour is not None:
feature_map["local_hour"] = local_hour
else:
timestamp = _parse_timestamp(snapshot.get("timestamp"))
if timestamp is not None:
feature_map["local_hour"] = float(timestamp.hour)
peak_code = _peak_status_code(snapshot.get("peak_status"))
if peak_code is not None:
feature_map["peak_status_code"] = peak_code
model_median, model_spread = _compute_model_summary(feature_map)
feature_map["model_median"] = model_median
feature_map["model_spread"] = model_spread
samples.append(
{
"city": _normalized_city_key(city_name),
"date": str(date_str),
"target": float(target),
"features": feature_map,
"vector": _features_to_vector(feature_map),
"history_count": int(meta.get("history_count") or 0),
"deb_prediction": _sf(record.get("deb_prediction")),
"forecasts": {
key: _sf(value)
for key, value in forecasts.items()
if _sf(value) is not None
},
"sample_source": "daily_record",
"settlement_source": truth_row.get("settlement_source"),
"settlement_station_code": truth_row.get("settlement_station_code"),
"truth_version": truth_row.get("truth_version"),
"truth_updated_by": truth_row.get("updated_by"),
"truth_updated_at": truth_row.get("truth_updated_at"),
}
)
samples.sort(key=lambda row: (row["date"], row["city"]))
return samples
def schema_payload(
*,
model_path: str,
sample_count: int,
train_count: int,
validation_count: int,
metrics: Dict[str, Any],
) -> Dict[str, Any]:
return {
"model_type": "LightGBMRegressor",
"target": "actual_high",
"horizon": "D0",
"feature_names": FEATURE_NAMES,
"base_model_columns": [column for _, column in BASE_MODEL_COLUMNS],
"model_path": model_path,
"sample_count": sample_count,
"train_count": train_count,
"validation_count": validation_count,
"metrics": metrics,
"generated_at": datetime.utcnow().isoformat() + "Z",
}