fix: use multi-model peak window for gaussian probability

2026-05-28 10:03:29 +08:00
parent fb4a9d7e09
commit be549daf50
4 changed files with 148 additions and 18 deletions
@@ -48,6 +48,91 @@ def _sf(v):
        return None


+def _median(values: List[float]) -> Optional[float]:
+    if not values:
+        return None
+    sorted_values = sorted(values)
+    mid = len(sorted_values) // 2
+    if len(sorted_values) % 2:
+        return sorted_values[mid]
+    return (sorted_values[mid - 1] + sorted_values[mid]) / 2.0
+
+
+def _peak_hours_from_hourly_values(
+    hourly_values: List[Tuple[str, float]],
+    *,
+    tolerance: float = 0.3,
+) -> List[str]:
+    if not hourly_values:
+        return []
+    peak_value = max(value for _, value in hourly_values)
+    return [
+        time_part
+        for time_part, value in hourly_values
+        if abs(value - peak_value) <= tolerance
+    ]
+
+
+def _resolve_peak_hours(
+    weather_data: dict,
+    local_date_str: str,
+    open_meteo_times: Optional[List[Any]] = None,
+    open_meteo_temps: Optional[List[Any]] = None,
+    open_meteo_peak: Optional[Any] = None,
+) -> List[str]:
+    """Resolve the local high-temperature window, preferring multi-model hourly consensus."""
+    multi_model = weather_data.get("multi_model") if isinstance(weather_data, dict) else {}
+    if isinstance(multi_model, dict):
+        hourly_times = multi_model.get("hourly_times") or []
+        hourly_forecasts = multi_model.get("hourly_forecasts") or {}
+        if isinstance(hourly_forecasts, dict) and hourly_times:
+            hourly_values: List[Tuple[str, float]] = []
+            for idx, raw_time in enumerate(hourly_times):
+                t_str = str(raw_time or "")
+                if not t_str.startswith(local_date_str) or "T" not in t_str:
+                    continue
+                time_part = t_str.split("T", 1)[1][:5]
+                try:
+                    hour = int(time_part[:2])
+                except Exception:
+                    continue
+                if not 8 <= hour <= 19:
+                    continue
+                values = []
+                for model_name, series in hourly_forecasts.items():
+                    if _is_excluded_model_name(model_name):
+                        continue
+                    if not isinstance(series, (list, tuple)) or idx >= len(series):
+                        continue
+                    value = _sf(series[idx])
+                    if value is not None:
+                        values.append(value)
+                median_value = _median(values)
+                if median_value is not None:
+                    hourly_values.append((time_part, median_value))
+            peak_hours = _peak_hours_from_hourly_values(hourly_values)
+            if peak_hours:
+                return peak_hours
+
+    om_peak = _sf(open_meteo_peak)
+    if open_meteo_times and open_meteo_temps and om_peak is not None:
+        peak_hours = []
+        for t_raw, temp_raw in zip(open_meteo_times, open_meteo_temps):
+            t_str = str(t_raw or "")
+            temp = _sf(temp_raw)
+            if temp is None or not t_str.startswith(local_date_str) or "T" not in t_str:
+                continue
+            time_part = t_str.split("T", 1)[1][:5]
+            try:
+                hour = int(time_part[:2])
+            except Exception:
+                continue
+            if 8 <= hour <= 19 and abs(temp - om_peak) <= 0.2:
+                peak_hours.append(time_part)
+        return peak_hours
+    return []
+
+
 def _resolve_settlement_source_label(city_name: Optional[str]) -> str:
    if not city_name:
        return "METAR"
@@ -420,16 +505,16 @@ def analyze_weather_trend(

    is_cooling = trend_direction == "falling"

-    om_today = daily.get("temperature_2m_max", [None])[0]
+    om_today = _sf(current_forecasts.get("Open-Meteo"))

    # === Peak hours ===
-    peak_hours = []
-    if times and temps and om_today is not None:
-        for t_str, temp in zip(times, temps):
-            if t_str.startswith(local_date_str) and abs(temp - om_today) <= 0.2:
-                hour = int(t_str.split("T")[1][:2])
-                if 8 <= hour <= 19:
-                    peak_hours.append(t_str.split("T")[1][:5])
+    peak_hours = _resolve_peak_hours(
+        weather_data,
+        local_date_str,
+        times,
+        temps,
+        om_today,
+    )
    if peak_hours:
        first_peak_h = int(peak_hours[0].split(":")[0])
        last_peak_h = int(peak_hours[-1].split(":")[0])
@@ -53,3 +53,10 @@ def test_probability_engine_uses_enriched_multi_model_snapshot():
    source = (ROOT / "web" / "analysis_service.py").read_text(encoding="utf-8")

    assert 'raw["multi_model"] = mm' in source
+
+
+def test_city_detail_peak_window_uses_shared_multi_model_resolver():
+    source = (ROOT / "web" / "analysis_service.py").read_text(encoding="utf-8")
+
+    assert "from src.analysis.trend_engine import _resolve_peak_hours" in source
+    assert "peak_hours = _resolve_peak_hours(" in source
@@ -25,6 +25,7 @@ def _make_weather_data(
    local_time="2026-03-04 14:30",
    recent_temps=None,
    multi_model=None,
+    multi_model_hourly=None,
    recent_obs=None,
 ):
    """Build a minimal weather_data dict for testing."""
@@ -60,7 +61,7 @@ def _make_weather_data(
            "p10": ens_p10,
            "p90": ens_p90,
        },
-        "multi_model": {"forecasts": multi_model or {}},
+        "multi_model": {"forecasts": multi_model or {}, **(multi_model_hourly or {})},
        "nws": {},
    }
    return data
@@ -153,16 +154,58 @@ class TestMuCalculation:
            recent_temps=[("09:00", 25.0), ("08:00", 24.0), ("07:00", 23.0)],
            multi_model={},
        )
-        data["open-meteo"]["utc_offset"] = 8 * 60 * 60
        data["open-meteo"]["daily"]["time"] = ["2026-05-27", "2026-05-28"]
        data["open-meteo"]["daily"]["temperature_2m_max"] = [24.0, 31.0]
        data["open-meteo"]["hourly"]["time"] = [f"2026-05-28T{h:02d}:00" for h in range(24)]

-        _, _, sd = analyze_weather_trend(data, "°C", "wuhan")
+        _, _, sd = analyze_weather_trend(data, "°C", "test_city")

        assert sd["current_forecasts"]["Open-Meteo"] == 31.0
        assert sd["mu"] is not None and sd["mu"] >= 30.0

+    @patch("src.analysis.trend_engine.calculate_dynamic_weights", return_value=(None, ""))
+    @patch("src.analysis.trend_engine.get_deb_accuracy", return_value=None)
+    @patch("src.analysis.trend_engine.update_daily_record")
+    def test_multi_model_peak_window_prevents_early_open_meteo_bust(
+        self, _udr, _deb_acc, _dw
+    ):
+        """If Open-Meteo peaks early but multi-models peak later, μ must not anchor to morning actuals."""
+        hourly_times = [f"2026-03-04T{h:02d}:00" for h in range(24)]
+        data = _make_weather_data(
+            cur_temp=25.0,
+            max_so_far=25.0,
+            om_today_high=31.5,
+            ens_median=None,
+            ens_p10=None,
+            ens_p90=None,
+            local_time="2026-03-04 09:50",
+            recent_temps=[("09:00", 25.0), ("08:00", 24.5), ("07:00", 24.0)],
+            multi_model={
+                "ECMWF": 30.8,
+                "GFS": 32.0,
+                "ICON": 28.8,
+                "GEM": 29.8,
+            },
+            multi_model_hourly={
+                "hourly_times": hourly_times,
+                "hourly_forecasts": {
+                    "ECMWF": [24, 24, 24, 24, 24, 24, 24, 25, 26, 27, 28, 29, 30, 30.5, 30.8, 30.8, 30.4, 29, 28, 27, 26, 25, 24, 24],
+                    "GFS": [24, 24, 24, 24, 24, 24, 24, 25, 26, 27, 29, 30, 31, 31.5, 32.0, 32.0, 31.2, 30, 28, 27, 26, 25, 24, 24],
+                    "ICON": [23, 23, 23, 23, 23, 23, 23, 24, 25, 26, 27, 28, 28.5, 28.8, 28.8, 28.8, 28.4, 28, 27, 26, 25, 24, 23, 23],
+                    "GEM": [24, 24, 24, 24, 24, 24, 24, 25, 26, 27, 28, 29, 29.5, 29.8, 29.8, 29.8, 29.4, 29, 28, 27, 26, 25, 24, 24],
+                },
+            },
+        )
+        data["open-meteo"]["hourly"]["time"] = hourly_times
+        data["open-meteo"]["hourly"]["temperature_2m"] = [
+            23, 23, 23, 23, 23, 24, 25, 27, 30, 31.5, 29, 28, 27, 26, 25, 24, 24, 23, 23, 22, 22, 22, 22, 22
+        ]
+
+        _, _, sd = analyze_weather_trend(data, "°C", "test_city")
+
+        assert sd["peak_status"] == "before"
+        assert sd["mu"] is not None and sd["mu"] >= 29.0
+

 # ─── Tests: Dead Market ───

@@ -30,6 +30,7 @@ from src.analysis.deb_hourly_correction import (
    get_cached_hourly_peak_corrector,
 )
 from src.analysis.settlement_rounding import apply_city_settlement
+from src.analysis.trend_engine import _resolve_peak_hours
 from src.data_collection.country_networks import build_country_network_snapshot
 from src.data_collection.city_registry import ALIASES, CITY_REGISTRY
 from src.data_collection.city_time import get_city_utc_offset_seconds
@@ -1187,13 +1188,7 @@ def _analyze(
            h_lifted_index = [None for _ in parsed_obs]
            h_boundary_layer_height = [None for _ in parsed_obs]

-    peak_hours = []
-    if h_times and h_temps and om_today is not None:
-        for ts, tmp in zip(h_times, h_temps):
-            if ts.startswith(local_date_str) and abs(tmp - om_today) <= 0.2:
-                hr = int(ts.split("T")[1][:2])
-                if 8 <= hr <= 19:
-                    peak_hours.append(ts.split("T")[1][:5])
+    peak_hours = _resolve_peak_hours(raw, local_date_str, h_times, h_temps, om_today)

    first_peak_h = int(peak_hours[0].split(":")[0]) if peak_hours else 13
    last_peak_h = int(peak_hours[-1].split(":")[0]) if peak_hours else 15