From fe6f8b43b02810c785a9da441d5de2f26f290ed7 Mon Sep 17 00:00:00 2001
From: "2569718930@qq.com" <2569718930@qq.com>
Date: Thu, 21 May 2026 21:03:49 +0800
Subject: [PATCH] =?UTF-8?q?DEB=20=E6=8E=A5=E5=85=A5=E5=B0=8F=E6=97=B6?=
 =?UTF-8?q?=E7=BA=A7=E8=AF=AF=E5=B7=AE=E8=AE=A1=E7=AE=97=EF=BC=8C=E5=A4=9A?=
 =?UTF-8?q?=E6=A8=A1=E5=9E=8B=E6=9D=83=E9=87=8D=E5=9F=BA=E4=BA=8E=E6=AF=8F?=
 =?UTF-8?q?=E6=97=A5+=E5=B0=8F=E6=97=B6=20MAE=20=E8=9E=8D=E5=90=88?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

新增 compute_hourly_model_errors 聚合逐模型小时 MAE/RMSE

新增 _blend_mae 按样本数加权混合每日/小时误差（24样本=70%小时权重）

calculate_dynamic_weights 从 daily_record 读取 hourly_error 参与权重计算

update_daily_record 接受并持久化 hourly_error 字段

Tested: ruff check, pytest 186/186
---
 src/analysis/deb_algorithm.py  | 116 ++++++++++++++++++++++++++++++++-
 tests/test_deb_model_family.py |  57 ++++++++++++++++
 2 files changed, 172 insertions(+), 1 deletion(-)

diff --git a/src/analysis/deb_algorithm.py b/src/analysis/deb_algorithm.py
index d50014e5..482237d5 100644
--- a/src/analysis/deb_algorithm.py
+++ b/src/analysis/deb_algorithm.py
@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import os
 import json
 from datetime import datetime, timedelta
@@ -140,6 +142,83 @@ def _collapse_forecasts_for_deb(current_forecasts):
     return collapsed
 
 
+def compute_hourly_model_errors(
+    hourly_forecasts: dict[str, list[float]],
+    hourly_actuals: list[float],
+) -> dict[str, dict[str, float]]:
+    """
+    Compute per-model hourly error aggregation for a single day.
+
+    hourly_forecasts: {model_name: [t0, t1, ..., tN]}  — N-hour forecast per model
+    hourly_actuals:  [t0, t1, ..., tM]                 — actual hourly temps
+
+    Returns: {model_name: {"mae": float, "rmse": float, "samples": int}}
+    """
+    if not hourly_actuals or not hourly_forecasts:
+        return {}
+
+    n = min(len(hourly_actuals), min(len(v) for v in hourly_forecasts.values()))
+    if n < 6:  # require at least 6 valid hours
+        return {}
+
+    actuals = hourly_actuals[:n]
+    result: dict[str, dict[str, float]] = {}
+
+    for model, preds in hourly_forecasts.items():
+        if not isinstance(preds, (list, tuple)) or len(preds) < n:
+            continue
+        try:
+            valid = preds[:n]
+        except (TypeError, IndexError):
+            continue
+
+        abs_errors = []
+        sq_errors = []
+        for p, a in zip(valid, actuals):
+            try:
+                pv = float(p)
+                av = float(a)
+            except (TypeError, ValueError):
+                continue
+            abs_errors.append(abs(pv - av))
+            sq_errors.append((pv - av) ** 2)
+
+        samples = len(abs_errors)
+        if samples < 6:
+            continue
+
+        mae = sum(abs_errors) / samples
+        rmse = (sum(sq_errors) / samples) ** 0.5
+        result[model] = {"mae": round(mae, 2), "rmse": round(rmse, 2), "samples": samples}
+
+    return result
+
+
+def _blend_mae(daily_mae: float, hourly_error: dict[str, float] | None) -> float:
+    """
+    Blend daily MAE and hourly MAE into a single error metric.
+
+    hourly_error  = {"mae": float, "rmse": float, "samples": int}
+
+    Weight split: daily=0.3, hourly=0.7 when hourly samples >= 24;
+    scales linearly from daily-only (0 samples) to full blend (≥24 samples).
+    """
+    if not hourly_error or not isinstance(hourly_error, dict):
+        return daily_mae
+
+    h_mae = float(hourly_error.get("mae", daily_mae))
+    h_samples = int(hourly_error.get("samples", 0))
+
+    if h_samples < 6:
+        return daily_mae
+
+    # Blend ratio: the more hourly samples, the more we trust hourly MAE
+    h_weight = min(0.7, h_samples / 24 * 0.7)
+    d_weight = 1.0 - h_weight
+
+    return round(daily_mae * d_weight + h_mae * h_weight, 2)
+
+
 def load_history(filepath):
     global _history_cache, _history_mtime
     mode = get_state_storage_mode()
@@ -800,6 +879,7 @@ def update_daily_record(
     probability_features=None,
     shadow_probabilities=None,
     calibration_summary=None,
+    hourly_error=None,
 ):
     """
     保存/更新某城市某天的各个模型预报与最终实测值
@@ -809,6 +889,7 @@ def update_daily_record(
     mu: float, 概率引擎中心值（用于 μ MAE 追踪）
     probabilities: list[dict], 概率分布快照（用于 Brier Score 校准）
         例如 [{"value": 25, "probability": 0.8}, {"value": 26, "probability": 0.2}]
+    hourly_error: dict, 逐模型小时级误差聚合 {"ECMWF": {"mae": 1.1, "rmse": 1.4, "samples": 68}, ...}
     """
     project_root = os.path.dirname(
         os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
@@ -882,6 +963,27 @@ def update_daily_record(
             merged_forecasts[model_name] = model_value
         elif model_name not in merged_forecasts:
             merged_forecasts[model_name] = model_value
+    old_hourly_error = existing.get("hourly_error")
+    # Merge hourly_error: keep existing per-model data and overlay new values
+    merged_hourly_error = dict(old_hourly_error) if isinstance(old_hourly_error, dict) else {}
+    if isinstance(hourly_error, dict):
+        for model, err in hourly_error.items():
+            if isinstance(err, dict) and all(
+                k in err for k in ("mae", "samples")
+            ):
+                # Prefer the entry with more samples
+                old_entry = merged_hourly_error.get(model)
+                if (
+                    not isinstance(old_entry, dict)
+                    or int(err.get("samples", 0)) >= int(old_entry.get("samples", 0))
+                ):
+                    merged_hourly_error[model] = {
+                        "mae": round(float(err["mae"]), 2),
+                        "rmse": round(float(err.get("rmse", err["mae"])), 2),
+                        "samples": int(err["samples"]),
+                    }
+    next_hourly_error = merged_hourly_error if merged_hourly_error else None
+
     next_mu = round(mu, 2) if mu is not None else None
     if (
         old_actual == actual_high
@@ -901,6 +1003,7 @@ def update_daily_record(
             compact_calibration is None
             or existing.get("probability_calibration") == compact_calibration
         )
+        and old_hourly_error == next_hourly_error
     ):
         return
 
@@ -925,6 +1028,8 @@ def update_daily_record(
         existing["shadow_prob_snapshot"] = compact_shadow_probs
     if compact_calibration is not None:
         existing["probability_calibration"] = compact_calibration
+    if next_hourly_error is not None:
+        existing["hourly_error"] = next_hourly_error
 
     if actual_high is not None:
         try:
@@ -1025,6 +1130,7 @@ def calculate_dynamic_weights(city_name, current_forecasts, lookback_days=7, dec
         record = city_data[date_str]
         actual = record.get("actual_high")
         past_forecasts = record.get("forecasts", {})
+        past_hourly_error = record.get("hourly_error")
 
         if actual is None:
             continue
@@ -1038,7 +1144,15 @@ def calculate_dynamic_weights(city_name, current_forecasts, lookback_days=7, dec
                     av = float(actual)
                 except (TypeError, ValueError):
                     continue
-                errors[model].append((abs(pv - av), decay_weight))
+                daily_error = abs(pv - av)
+                # Blend with hourly error when available
+                h_err = (
+                    past_hourly_error.get(model)
+                    if isinstance(past_hourly_error, dict)
+                    else None
+                )
+                blended_error = _blend_mae(daily_error, h_err)
+                errors[model].append((blended_error, decay_weight))
 
         days_used += 1
         if days_used >= lookback_days:
diff --git a/tests/test_deb_model_family.py b/tests/test_deb_model_family.py
index 082f7810..e420a465 100644
--- a/tests/test_deb_model_family.py
+++ b/tests/test_deb_model_family.py
@@ -89,3 +89,60 @@ def test_deb_weighted_path_uses_deduped_family_values(monkeypatch):
     assert blended < 30.0
     assert "ICON-D2" in info
     assert "家族去重" in info
+
+
+def test_compute_hourly_model_errors_basic():
+    from src.analysis.deb_algorithm import compute_hourly_model_errors
+
+    hourly_forecasts = {
+        "ECMWF": [15.0, 15.5, 16.0, 16.5, 17.0, 17.5],
+        "GFS":   [14.5, 15.0, 15.5, 16.0, 16.5, 17.0],
+    }
+    hourly_actuals = [15.2, 15.8, 16.3, 16.8, 17.4, 18.0]
+
+    result = compute_hourly_model_errors(hourly_forecasts, hourly_actuals)
+
+    assert "ECMWF" in result
+    assert "GFS" in result
+    assert result["ECMWF"]["samples"] == 6
+    assert result["GFS"]["samples"] == 6
+    # ECMWF predicted lower, should have smaller MAE
+    assert result["ECMWF"]["mae"] < result["GFS"]["mae"]
+    assert result["ECMWF"]["rmse"] > 0
+    assert result["GFS"]["rmse"] > 0
+
+
+def test_compute_hourly_model_errors_too_few_samples():
+    from src.analysis.deb_algorithm import compute_hourly_model_errors
+
+    hourly_forecasts = {"ECMWF": [15.0, 15.5]}
+    hourly_actuals = [15.0, 15.5]
+
+    result = compute_hourly_model_errors(hourly_forecasts, hourly_actuals)
+    assert result == {}  # < 6 samples, returns empty
+
+
+def test_blend_mae_with_hourly():
+    from src.analysis.deb_algorithm import _blend_mae
+
+    # Full 24 hourly samples → 0.7 hourly weight
+    h_err = {"mae": 0.8, "rmse": 1.0, "samples": 24}
+    blended = _blend_mae(1.5, h_err)
+    # 1.5 * 0.3 + 0.8 * 0.7 = 0.45 + 0.56 = 1.01
+    assert abs(blended - 1.01) < 0.01
+
+    # 12 hourly samples → ~0.35 hourly weight
+    h_err2 = {"mae": 1.0, "rmse": 1.2, "samples": 12}
+    blended2 = _blend_mae(2.0, h_err2)
+    # 12/24 * 0.7 = 0.35 hourly weight
+    expected = 2.0 * (1 - 0.35) + 1.0 * 0.35  # 1.30 + 0.35 = 1.65
+    assert abs(blended2 - expected) < 0.01
+
+    # No hourly error → returns daily MAE unchanged
+    blended3 = _blend_mae(2.5, None)
+    assert blended3 == 2.5
+
+    # Too few samples (< 6) → returns daily MAE unchanged
+    h_err4 = {"mae": 0.5, "samples": 3}
+    blended4 = _blend_mae(3.0, h_err4)
+    assert blended4 == 3.0