PolyWeather/scripts/backtest_metar_calibrated_path.py

#!/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())