#!/usr/bin/env python3 """Backtest DEB baseline vs METAR/observation-calibrated intraday path. This script intentionally mirrors the frontend chart logic at a data-science level: DEB baseline path = hourly forecast curve + (DEB daily high - OM daily high) calibrated path = DEB path + recent observation bias * fade-to-evening It uses only local data. The best dataset is SQLite runtime state because it can contain: - open_meteo_cache_store: hourly forecast curves - official_intraday_observations_store: intraday anchor observations - daily_records_store / truth_records_store: final actual high If a city/date lacks any of those pieces, it is skipped and reported. """ from __future__ import annotations import argparse import csv import json import math import sqlite3 import statistics import sys from dataclasses import dataclass from datetime import datetime from pathlib import Path from typing import Any, Iterable ROOT = Path(__file__).resolve().parents[1] if str(ROOT) not in sys.path: sys.path.insert(0, str(ROOT)) try: from src.data_collection.city_registry import CITY_REGISTRY from src.analysis.settlement_rounding import apply_city_settlement except Exception: # pragma: no cover - script fallback for partial envs CITY_REGISTRY = {} def apply_city_settlement(_city: str, value: float | None) -> int | None: return None if value is None else round(value) def sf(value: Any) -> float | None: try: if value is None or value == "": return None num = float(value) return num if math.isfinite(num) else None except Exception: return None def hm_to_minutes(value: str | None) -> int | None: if not value: return None text = str(value).strip() if "T" in text: text = text.split("T", 1)[1] text = text[:5] try: hh, mm = text.split(":")[:2] h = int(hh) m = int(mm) if not (0 <= h <= 23 and 0 <= m <= 59): return None return h * 60 + m except Exception: return None def interp(times: list[str], values: list[float | None], minute: int) -> float | None: pts: list[tuple[int, float]] = [] for t, v in zip(times, values): m = hm_to_minutes(t) y = sf(v) if m is not None and y is not None: pts.append((m, y)) if not pts: return None pts.sort() if minute <= pts[0][0]: return pts[0][1] if minute >= pts[-1][0]: return pts[-1][1] for (lm, lv), (rm, rv) in zip(pts, pts[1:]): if lm <= minute <= rm: if rm == lm: return rv ratio = (minute - lm) / (rm - lm) return lv + (rv - lv) * ratio return pts[-1][1] def clamp_delta(value: float, lo: float = -4.0, hi: float = 4.0) -> float: return min(max(value, lo), hi) @dataclass class Observation: time: str temp: float @dataclass class SampleResult: city: str date: str current_time: str obs_count: int actual_high: float deb_high: float calibrated_high: float deb_abs_error: float calibrated_abs_error: float delta_vs_deb: float bucket_deb_hit: bool | None bucket_calibrated_hit: bool | None def dedupe_observations(rows: Iterable[Observation]) -> list[Observation]: by_time: dict[str, Observation] = {} for row in rows: minute = hm_to_minutes(row.time) if minute is None: continue key = f"{minute // 60:02d}:{minute % 60:02d}" existing = by_time.get(key) if existing is None or row.temp >= existing.temp: by_time[key] = Observation(key, row.temp) return sorted(by_time.values(), key=lambda r: hm_to_minutes(r.time) or 0) def calibrated_future_path( *, times: list[str], deb_path: list[float | None], observations: list[Observation], current_minute: int, reversion_minute: int, ) -> tuple[list[float | None], float | None]: usable: list[tuple[int, float]] = [] for obs in dedupe_observations(observations): minute = hm_to_minutes(obs.time) if minute is None or minute > current_minute + 30: continue expected = interp(times, deb_path, minute) if expected is None: continue usable.append((minute, clamp_delta(obs.temp - expected))) usable = usable[-3:] if not usable: return [None for _ in times], None total = 0.0 weight_total = 0.0 for idx, (_minute, delta) in enumerate(usable): weight = idx + 1 total += delta * weight weight_total += weight adjustment = round(clamp_delta(total / max(weight_total, 1.0)), 1) last_minute = next((m for m in reversed([hm_to_minutes(t) for t in times]) if m is not None), current_minute + 360) return_to = reversion_minute if reversion_minute > current_minute else last_minute if return_to <= current_minute: return_to = current_minute + 360 out: list[float | None] = [] for t, base in zip(times, deb_path): minute = hm_to_minutes(t) if minute is None or minute < current_minute or base is None: out.append(None) continue progress = min(max((minute - current_minute) / max(return_to - current_minute, 1), 0.0), 1.0) decay = (1 - progress) ** 1.35 out.append(round(base + adjustment * decay, 1)) return out, adjustment def connect(db_path: Path) -> sqlite3.Connection: con = sqlite3.connect(str(db_path)) con.row_factory = sqlite3.Row return con def get_actual_high(con: sqlite3.Connection, city: str, date: str) -> float | None: row = con.execute( "select actual_high from daily_records_store where city=? and target_date=?", (city, date), ).fetchone() if row and sf(row["actual_high"]) is not None: return sf(row["actual_high"]) row = con.execute( "select actual_high from truth_records_store where city=? and target_date=? and is_final=1 order by updated_at desc limit 1", (city, date), ).fetchone() return sf(row["actual_high"]) if row else None def get_daily_record(con: sqlite3.Connection, city: str, date: str) -> dict[str, Any] | None: row = con.execute( "select deb_prediction, payload_json from daily_records_store where city=? and target_date=?", (city, date), ).fetchone() if not row: return None payload = {} try: payload = json.loads(row["payload_json"] or "{}") except Exception: payload = {} payload.setdefault("deb_prediction", row["deb_prediction"]) return payload def cache_key_for_city(city: str) -> str | None: meta = CITY_REGISTRY.get(city) or {} lat = sf(meta.get("lat")) lon = sf(meta.get("lon")) if lat is None or lon is None: return None unit = "f" if meta.get("use_fahrenheit") else "c" return f"{lat:.4f}:{lon:.4f}:14:{unit}" def load_hourly_forecast(con: sqlite3.Connection, city: str, date: str) -> tuple[list[str], list[float | None], float | None, str | None]: key = cache_key_for_city(city) if not key: return [], [], None, None row = con.execute( "select payload_json, updated_at from open_meteo_cache_store where source_kind='forecast' and cache_key=? order by updated_at desc limit 1", (key,), ).fetchone() if not row: return [], [], None, None try: payload = json.loads(row["payload_json"] or "{}") except Exception: return [], [], None, None data = payload.get("data") or payload hourly = data.get("hourly") or {} raw_times = hourly.get("time") or hourly.get("times") or [] raw_temps = hourly.get("temperature_2m") or hourly.get("temps") or [] times: list[str] = [] temps: list[float | None] = [] for ts, temp in zip(raw_times, raw_temps): text = str(ts) if not text.startswith(date): continue times.append(text.split("T", 1)[1][:5] if "T" in text else text[:5]) temps.append(sf(temp)) om_high = max([v for v in temps if v is not None], default=None) updated_at = None try: updated_at = datetime.fromtimestamp(float(row["updated_at"])).isoformat() except Exception: pass return times, temps, om_high, updated_at def station_codes_for_city(city: str) -> set[str]: meta = CITY_REGISTRY.get(city) or {} values = [ meta.get("settlement_station_code"), meta.get("icao"), *(meta.get("settlement_station_candidates") or []), ] return {str(v).strip().upper() for v in values if str(v or "").strip()} def load_observations(con: sqlite3.Connection, city: str, date: str) -> list[Observation]: codes = station_codes_for_city(city) if not codes: return [] placeholders = ",".join("?" for _ in codes) rows = con.execute( f""" select observation_time, value from official_intraday_observations_store where target_date=? and upper(station_code) in ({placeholders}) order by observation_time asc """, (date, *sorted(codes)), ).fetchall() obs = [] for row in rows: temp = sf(row["value"]) time = str(row["observation_time"] or "")[:5] if temp is not None and hm_to_minutes(time) is not None: obs.append(Observation(time, temp)) return dedupe_observations(obs) def bucket_hit(city: str, predicted: float, actual: float) -> bool | None: try: return apply_city_settlement(city, predicted) == apply_city_settlement(city, actual) except Exception: return None def evaluate_city_date(con: sqlite3.Connection, city: str, date: str, min_obs: int) -> list[SampleResult]: daily = get_daily_record(con, city, date) if not daily: return [] actual_high = get_actual_high(con, city, date) deb_high = sf(daily.get("deb_prediction")) if actual_high is None or deb_high is None: return [] times, temps, om_high, _updated_at = load_hourly_forecast(con, city, date) if not times or om_high is None: # Fallback: use Open-Meteo daily forecast from the daily record only. # This cannot evaluate path shape, so skip rather than pretend. return [] offset = deb_high - om_high deb_path = [round(t + offset, 1) if t is not None else None for t in temps] observations = load_observations(con, city, date) if len(observations) < min_obs: return [] sunset = "18:00" reversion_minute = hm_to_minutes(sunset) or 18 * 60 results: list[SampleResult] = [] for idx in range(min_obs - 1, len(observations)): current_obs = observations[idx] current_minute = hm_to_minutes(current_obs.time) if current_minute is None: continue used_obs = observations[: idx + 1] calibrated_path, adjustment = calibrated_future_path( times=times, deb_path=deb_path, observations=used_obs, current_minute=current_minute, reversion_minute=reversion_minute, ) future_values = [v for v in calibrated_path if v is not None] observed_so_far = max(o.temp for o in used_obs) calibrated_high = max([observed_so_far, *future_values], default=observed_so_far) results.append( SampleResult( city=city, date=date, current_time=current_obs.time, obs_count=len(used_obs), actual_high=actual_high, deb_high=deb_high, calibrated_high=calibrated_high, deb_abs_error=abs(deb_high - actual_high), calibrated_abs_error=abs(calibrated_high - actual_high), delta_vs_deb=adjustment if adjustment is not None else 0.0, bucket_deb_hit=bucket_hit(city, deb_high, actual_high), bucket_calibrated_hit=bucket_hit(city, calibrated_high, actual_high), ) ) return results def load_path_snapshot_rows( con: sqlite3.Connection, *, cities: list[str], dates_filter: set[str], ) -> list[dict[str, Any]]: try: con.execute("select 1 from intraday_path_snapshots_store limit 1").fetchone() except sqlite3.Error: return [] params: list[Any] = [] clauses: list[str] = [] if cities: clauses.append("city in (" + ",".join("?" for _ in cities) + ")") params.extend(cities) if dates_filter: clauses.append("target_date in (" + ",".join("?" for _ in dates_filter) + ")") params.extend(sorted(dates_filter)) where = " where " + " and ".join(clauses) if clauses else "" rows = con.execute( f"select payload_json from intraday_path_snapshots_store{where} order by id asc", params, ).fetchall() out: list[dict[str, Any]] = [] for row in rows: try: payload = json.loads(row["payload_json"]) except Exception: continue if isinstance(payload, dict): out.append(payload) return out def observation_rows_from_snapshot(snapshot: dict[str, Any]) -> list[Observation]: rows: list[Observation] = [] for key in ("metar_today_obs", "settlement_today_obs"): raw_rows = snapshot.get(key) if not isinstance(raw_rows, list): continue for item in raw_rows: if not isinstance(item, dict): continue temp = sf(item.get("temp")) time_text = str(item.get("time") or "").strip()[:5] if temp is not None and hm_to_minutes(time_text) is not None: rows.append(Observation(time_text, temp)) return dedupe_observations(rows) def evaluate_path_snapshot( con: sqlite3.Connection, snapshot: dict[str, Any], min_obs: int, ) -> SampleResult | None: city = str(snapshot.get("city") or "").strip().lower() date = str(snapshot.get("target_date") or snapshot.get("date") or "").strip() if not city or not date: return None actual_high = get_actual_high(con, city, date) deb_high = sf(snapshot.get("deb_prediction")) if actual_high is None or deb_high is None: return None path = snapshot.get("deb_base_path") or {} times = path.get("times") if isinstance(path, dict) else [] deb_path = path.get("temps") if isinstance(path, dict) else [] if not isinstance(times, list) or not isinstance(deb_path, list) or not times: return None deb_values = [sf(v) for v in deb_path] observations = observation_rows_from_snapshot(snapshot) if len(observations) < min_obs: return None local_time = str(snapshot.get("local_time") or "").strip() current_minute = hm_to_minutes(local_time) if current_minute is None: current_minute = hm_to_minutes(observations[-1].time) if current_minute is None: return None forecast = snapshot.get("forecast") if isinstance(snapshot.get("forecast"), dict) else {} reversion_minute = hm_to_minutes(forecast.get("sunset")) or hm_to_minutes("18:00") or 18 * 60 calibrated_path, adjustment = calibrated_future_path( times=[str(t) for t in times], deb_path=deb_values, observations=observations, current_minute=current_minute, reversion_minute=reversion_minute, ) future_values = [v for v in calibrated_path if v is not None] current = snapshot.get("current") if isinstance(snapshot.get("current"), dict) else {} max_so_far = sf(current.get("max_so_far")) observed_so_far = max([o.temp for o in observations] + ([max_so_far] if max_so_far is not None else [])) calibrated_high = max([observed_so_far, *future_values], default=observed_so_far) return SampleResult( city=city, date=date, current_time=local_time or observations[-1].time, obs_count=len(observations), actual_high=actual_high, deb_high=deb_high, calibrated_high=calibrated_high, deb_abs_error=abs(deb_high - actual_high), calibrated_abs_error=abs(calibrated_high - actual_high), delta_vs_deb=adjustment if adjustment is not None else 0.0, bucket_deb_hit=bucket_hit(city, deb_high, actual_high), bucket_calibrated_hit=bucket_hit(city, calibrated_high, actual_high), ) def summarize(samples: list[SampleResult]) -> dict[str, Any]: if not samples: return {"samples": 0} deb_errors = [s.deb_abs_error for s in samples] cal_errors = [s.calibrated_abs_error for s in samples] improved = [s for s in samples if s.calibrated_abs_error < s.deb_abs_error] worsened = [s for s in samples if s.calibrated_abs_error > s.deb_abs_error] deb_hits = [s.bucket_deb_hit for s in samples if s.bucket_deb_hit is not None] cal_hits = [s.bucket_calibrated_hit for s in samples if s.bucket_calibrated_hit is not None] return { "samples": len(samples), "city_dates": len({(s.city, s.date) for s in samples}), "deb_mae": round(statistics.mean(deb_errors), 3), "calibrated_mae": round(statistics.mean(cal_errors), 3), "mae_delta_cal_minus_deb": round(statistics.mean(cal_errors) - statistics.mean(deb_errors), 3), "improved_samples": len(improved), "worsened_samples": len(worsened), "unchanged_samples": len(samples) - len(improved) - len(worsened), "deb_bucket_hit_rate": round(sum(1 for x in deb_hits if x) / len(deb_hits), 3) if deb_hits else None, "calibrated_bucket_hit_rate": round(sum(1 for x in cal_hits if x) / len(cal_hits), 3) if cal_hits else None, } def write_csv(path: Path, samples: list[SampleResult]) -> None: path.parent.mkdir(parents=True, exist_ok=True) with path.open("w", newline="", encoding="utf-8") as f: writer = csv.DictWriter(f, fieldnames=list(SampleResult.__dataclass_fields__.keys())) writer.writeheader() for s in samples: writer.writerow(s.__dict__) def main() -> int: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--db", default=str(ROOT / "data" / "polyweather.db")) parser.add_argument("--city", action="append", help="City key; can be repeated. Defaults to cities found in daily_records_store.") parser.add_argument("--date", action="append", help="YYYY-MM-DD; can be repeated.") parser.add_argument("--min-obs", type=int, default=2) parser.add_argument( "--source", choices=("strict", "snapshots", "both"), default="both", help="strict uses reconstructed legacy stores; snapshots uses intraday_path_snapshots_store.", ) parser.add_argument("--output", default=str(ROOT / "tmp_metar_calibration_backtest.csv")) args = parser.parse_args() con = connect(Path(args.db)) cities = args.city if not cities: cities = [r[0] for r in con.execute("select distinct city from daily_records_store order by city").fetchall()] dates_filter = set(args.date or []) all_samples: list[SampleResult] = [] skipped = {"no_records_or_inputs": 0, "snapshots_unusable": 0} if args.source in {"strict", "both"}: for city in cities: rows = con.execute( "select distinct target_date from daily_records_store where city=? order by target_date", (city,), ).fetchall() for row in rows: date = row[0] if dates_filter and date not in dates_filter: continue samples = evaluate_city_date(con, city, date, args.min_obs) if samples: all_samples.extend(samples) else: skipped["no_records_or_inputs"] += 1 if args.source in {"snapshots", "both"}: snapshot_rows = load_path_snapshot_rows( con, cities=cities or [], dates_filter=dates_filter, ) for snapshot in snapshot_rows: sample = evaluate_path_snapshot(con, snapshot, args.min_obs) if sample: all_samples.append(sample) else: skipped["snapshots_unusable"] += 1 summary = summarize(all_samples) write_csv(Path(args.output), all_samples) print(json.dumps({"summary": summary, "skipped": skipped, "output": args.output}, ensure_ascii=False, indent=2)) if not all_samples: print( "No usable samples. Need strict store matches or rows in intraday_path_snapshots_store with later actual_high.", file=sys.stderr, ) return 2 return 0 if __name__ == "__main__": raise SystemExit(main())