#!/usr/bin/env python3
"""
Forecast Skill Evaluation — measures how well the propagation scoring algorithm
predicts contact distances as the weather profile ages.

For each contact, we look up HRRR profiles at 0h (analysis/ground truth),
1h, 3h, 6h, 12h, and 24h before the contact time. These lagged profiles
represent the weather data that would have been available N hours earlier.
We then compute how much the score↔distance correlation degrades with lead time.

Output:
  - Markdown table to stdout
  - docs/algo-reports/forecast-2026-08-01.json
  - docs/algo-reports/forecast-2026-08-01.md

Usage:
    python3 scripts/validate_forecast.py

Env vars:
    PROP_PROD_DB_URL — PostgreSQL connection string
"""

from __future__ import annotations

import datetime as dt
import json
import math
import os
import sys
from collections import defaultdict
from pathlib import Path

try:
    import psycopg
    from psycopg.rows import dict_row
except ImportError:
    sys.exit("psycopg is required: pip install 'psycopg[binary]>=3.1'")

REPO_ROOT = Path(__file__).resolve().parent.parent
WEIGHTS_PATH = REPO_ROOT / "priv" / "algo" / "band_weights.json"
REPORT_DIR = REPO_ROOT / "docs" / "algo-reports"

# ── Hardcoded band configuration (mirrors validate_algo.py) ────────────────────

SUNRISE_TABLE = [7.4, 7.3, 7.0, 6.7, 6.35, 6.25, 6.35, 6.65, 6.9, 7.1, 7.35, 7.45]

REFRACTIVITY_THRESHOLDS = [
    (-200, 98, 85),
    (-150, 92, 80),
    (-100, 82, 72),
    (-75, 68, 62),
    (-55, 55, 55),
    (-40, 48, 48),
]
REFRACTIVITY_DEFAULT = (42, 42)

HUMIDITY_BENEFICIAL_THRESHOLDS = [(4, 55), (7, 70), (10, 82), (14, 90), (18, 95), (22, 88)]
HUMIDITY_BENEFICIAL_DEFAULT = 75

BAND_CONFIGS: dict[int, dict] = {
    50: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.0, "rain_alpha": 1.0,
         "seasonal_base": {1: 30, 2: 32, 3: 25, 4: 50, 5: 70, 6: 95, 7: 95, 8: 70, 9: 65, 10: 70, 11: 60, 12: 25},
         "seasonal_adj": {}},
    144: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.0, "rain_alpha": 1.0,
          "seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
          "seasonal_adj": {}},
    222: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.0, "rain_alpha": 1.0,
          "seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
          "seasonal_adj": {}},
    432: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.0, "rain_alpha": 1.0,
          "seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
          "seasonal_adj": {}},
    902: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.0, "rain_alpha": 1.0,
          "seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
          "seasonal_adj": {}},
    1296: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.0, "rain_alpha": 1.0,
           "seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
           "seasonal_adj": {}},
    2304: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.001, "rain_alpha": 1.15,
           "seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
           "seasonal_adj": {}},
    3400: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.002, "rain_alpha": 1.20,
           "seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
           "seasonal_adj": {}},
    5760: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.005, "rain_alpha": 1.25,
           "seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
           "seasonal_adj": {}},
    10000: {"humidity_effect": "beneficial", "humidity_penalty": 0.0, "rain_k": 0.010, "rain_alpha": 1.28,
            "seasonal_base": {1: 38, 2: 40, 3: 22, 4: 55, 5: 68, 6: 90, 7: 95, 8: 75, 9: 78, 10: 88, 11: 78, 12: 25},
            "seasonal_adj": {}},
    24000: {"humidity_effect": "harmful", "humidity_penalty": 1.6, "rain_k": 0.070, "rain_alpha": 1.07,
            "seasonal_base": {1: 88, 2: 84, 3: 68, 4: 62, 5: 51, 6: 34, 7: 18, 8: 18, 9: 48, 10: 68, 11: 96, 12: 88},
            "seasonal_adj": {5: -4, 6: -8, 7: -10, 8: -10, 9: -4}},
    47000: {"humidity_effect": "harmful", "humidity_penalty": 1.0, "rain_k": 0.187, "rain_alpha": 0.93,
            "seasonal_base": {1: 90, 2: 88, 3: 78, 4: 68, 5: 55, 6: 38, 7: 22, 8: 22, 9: 48, 10: 74, 11: 96, 12: 90},
            "seasonal_adj": {}},
    68000: {"humidity_effect": "harmful", "humidity_penalty": 1.4, "rain_k": 0.310, "rain_alpha": 0.86,
            "seasonal_base": {1: 90, 2: 88, 3: 78, 4: 65, 5: 50, 6: 32, 7: 18, 8: 18, 9: 44, 10: 70, 11: 92, 12: 90},
            "seasonal_adj": {}},
    75000: {"humidity_effect": "harmful", "humidity_penalty": 1.2, "rain_k": 0.345, "rain_alpha": 0.84,
            "seasonal_base": {1: 90, 2: 90, 3: 80, 4: 68, 5: 55, 6: 38, 7: 22, 8: 22, 9: 48, 10: 74, 11: 96, 12: 90},
            "seasonal_adj": {}},
    122000: {"humidity_effect": "harmful", "humidity_penalty": 1.0, "rain_k": 0.498, "rain_alpha": 0.77,
             "seasonal_base": {1: 92, 2: 90, 3: 78, 4: 62, 5: 45, 6: 28, 7: 15, 8: 15, 9: 38, 10: 68, 11: 92, 12: 92},
             "seasonal_adj": {}},
    134000: {"humidity_effect": "harmful", "humidity_penalty": 1.3, "rain_k": 0.520, "rain_alpha": 0.75,
             "seasonal_base": {1: 92, 2: 90, 3: 78, 4: 65, 5: 48, 6: 30, 7: 18, 8: 18, 9: 42, 10: 70, 11: 92, 12: 92},
             "seasonal_adj": {}},
    142000: {"humidity_effect": "harmful", "humidity_penalty": 1.4, "rain_k": 0.530, "rain_alpha": 0.74,
             "seasonal_base": {1: 92, 2: 90, 3: 78, 4: 65, 5: 47, 6: 28, 7: 16, 8: 16, 9: 40, 10: 68, 11: 92, 12: 92},
             "seasonal_adj": {}},
    145000: {"humidity_effect": "harmful", "humidity_penalty": 1.5, "rain_k": 0.535, "rain_alpha": 0.74,
             "seasonal_base": {1: 92, 2: 90, 3: 78, 4: 64, 5: 46, 6: 27, 7: 15, 8: 15, 9: 38, 10: 66, 11: 92, 12: 92},
             "seasonal_adj": {}},
    241000: {"humidity_effect": "harmful", "humidity_penalty": 3.0, "rain_k": 0.550, "rain_alpha": 0.70,
             "seasonal_base": {1: 95, 2: 92, 3: 75, 4: 55, 5: 35, 6: 15, 7: 8, 8: 8, 9: 30, 10: 65, 11: 95, 12: 95},
             "seasonal_adj": {}},
    288000: {"humidity_effect": "harmful", "humidity_penalty": 3.5, "rain_k": 0.560, "rain_alpha": 0.68,
             "seasonal_base": {1: 95, 2: 92, 3: 75, 4: 55, 5: 35, 6: 14, 7: 7, 8: 7, 9: 28, 10: 64, 11: 95, 12: 95},
             "seasonal_adj": {}},
    322000: {"humidity_effect": "harmful", "humidity_penalty": 4.0, "rain_k": 0.570, "rain_alpha": 0.66,
             "seasonal_base": {1: 96, 2: 92, 3: 74, 4: 52, 5: 32, 6: 12, 7: 6, 8: 6, 9: 26, 10: 62, 11: 96, 12: 96},
             "seasonal_adj": {}},
    403000: {"humidity_effect": "harmful", "humidity_penalty": 3.0, "rain_k": 0.580, "rain_alpha": 0.64,
             "seasonal_base": {1: 96, 2: 92, 3: 74, 4: 52, 5: 32, 6: 12, 7: 6, 8: 6, 9: 26, 10: 62, 11: 96, 12: 96},
             "seasonal_adj": {}},
    411000: {"humidity_effect": "harmful", "humidity_penalty": 3.2, "rain_k": 0.580, "rain_alpha": 0.64,
             "seasonal_base": {1: 96, 2: 92, 3: 74, 4: 52, 5: 32, 6: 12, 7: 6, 8: 6, 9: 26, 10: 62, 11: 96, 12: 96},
             "seasonal_adj": {}},
}

FACTOR_ORDER = [
    "humidity", "time_of_day", "td_depression", "refractivity",
    "sky", "season", "wind", "rain", "pwat", "pressure",
]

# Lead times to evaluate (hours before contact)
LEAD_TIMES = [0, 1, 3, 6, 12, 24]

# ── Scoring functions (identical to validate_algo.py) ──────────────────────────

def absolute_humidity(temp_c: float, dewpoint_c: float) -> float:
    e_sat = 6.112 * math.exp(17.67 * dewpoint_c / (dewpoint_c + 243.5))
    return 217.0 * e_sat / (temp_c + 273.15)


def c_to_f(c: float) -> float:
    return c * 9.0 / 5.0 + 32.0


def score_humidity(abs_hum: float | None, band_cfg: dict) -> int:
    if abs_hum is None:
        return HUMIDITY_BENEFICIAL_DEFAULT
    effect = band_cfg["humidity_effect"]
    if effect == "beneficial":
        for max_h, s in HUMIDITY_BENEFICIAL_THRESHOLDS:
            if abs_hum < max_h:
                return s
        return HUMIDITY_BENEFICIAL_DEFAULT
    else:
        penalty = band_cfg["humidity_penalty"]
        r = abs_hum * penalty
        if r <= 6: return 100
        if r <= 9: return round(95 - (r - 6) / 3 * 20)
        if r <= 13: return round(75 - (r - 9) / 4 * 30)
        if r <= 18: return round(45 - (r - 13) / 5 * 35)
        return max(0, round(10 - (r - 18) * 2))


def score_time_of_day(utc_hour: int, utc_minute: int, month: int, lon: float) -> int:
    offset = lon / 15.0
    local = (utc_hour + utc_minute / 60.0 + offset + 24) % 24
    sunrise = SUNRISE_TABLE[month - 1]
    d = local - sunrise
    if -1.5 <= d <= 1.5: return 100
    if 1.5 < d <= 3.0: return 78
    if -3.0 <= d < -1.5: return 82
    if 3.0 < d <= 6.0: return 38
    if local >= 20 or local <= 1: return 72
    if d > 6: return 18
    return 55


def score_td_depression(temp_c: float | None, dewpoint_c: float | None, band_cfg: dict) -> int:
    if temp_c is None or dewpoint_c is None:
        return 50
    dep_f = c_to_f(temp_c) - c_to_f(dewpoint_c)
    effect = band_cfg["humidity_effect"]
    if effect == "beneficial":
        if dep_f < 3: return 40
        if dep_f < 8: return 75
        if dep_f < 14: return 85
        if dep_f < 22: return 70
        return 55
    else:
        if dep_f > 22: return 96
        if dep_f > 14: return 80
        if dep_f > 8: return 60
        if dep_f > 4: return 38
        return 18


def score_refractivity(gradient: float | None, band_cfg: dict) -> int:
    if gradient is None:
        return 50
    effect = band_cfg["humidity_effect"]
    for max_grad, ben, harm in REFRACTIVITY_THRESHOLDS:
        if gradient < max_grad:
            return ben if effect == "beneficial" else harm
    b_def, h_def = REFRACTIVITY_DEFAULT
    return b_def if effect == "beneficial" else h_def


def score_sky(cloud_pct: float | None) -> int:
    if cloud_pct is None:
        return 50
    if cloud_pct <= 6: return 100
    if cloud_pct <= 25: return 88
    if cloud_pct <= 50: return 60
    if cloud_pct <= 87: return 25
    return 5


def score_season(month: int, band_cfg: dict) -> int:
    base = band_cfg["seasonal_base"].get(month, 50)
    adj = band_cfg["seasonal_adj"].get(month, 0)
    return round(min(100, max(0, base + adj)))


def score_wind(speed_kts: float | None) -> int:
    if speed_kts is None:
        return 50
    if speed_kts < 5: return 100
    if speed_kts < 10: return 90
    if speed_kts < 15: return 75
    if speed_kts < 20: return 55
    if speed_kts < 25: return 35
    return 15


def score_rain(rate_mmhr: float | None, band_cfg: dict) -> int:
    if rate_mmhr is None or rate_mmhr == 0:
        return 100
    k = band_cfg["rain_k"]
    alpha = band_cfg["rain_alpha"]
    if k == 0:
        return 100
    gamma = k * (rate_mmhr ** alpha)
    if gamma < 0.1: return 95
    if gamma < 0.5: return 75
    if gamma < 1.0: return 50
    if gamma < 2.0: return 25
    if gamma < 5.0: return 10
    return 0


def score_pwat(pwat_mm: float | None, band_cfg: dict) -> int:
    if pwat_mm is None:
        return 60
    effect = band_cfg["humidity_effect"]
    if effect == "beneficial":
        if pwat_mm < 10: return 55
        if pwat_mm < 20: return 75
        if pwat_mm < 30: return 90
        if pwat_mm < 40: return 70
        return 50
    else:
        if pwat_mm < 10: return 95
        if pwat_mm < 20: return 80
        if pwat_mm < 30: return 60
        if pwat_mm < 40: return 35
        return 15


def score_pressure(pressure_mb: float | None) -> int:
    if pressure_mb is None:
        return 50
    if pressure_mb < 980: return 88
    if pressure_mb < 990: return 82
    if pressure_mb < 1000: return 70
    if pressure_mb < 1010: return 55
    if pressure_mb < 1020: return 40
    return 30


def compute_composite_score(profile: dict | None, contact: dict, band_cfg: dict,
                             weights: dict[str, float]) -> int | None:
    """Compute the composite score for a contact using a specific weather profile.

    If profile is None (missing), returns None.
    The contact's own metadata (month, hour, lon) is used for time-of-day and
    season factors; weather-driven factors use the profile.
    """
    if profile is None:
        return None

    month = contact["month"]
    abs_hum = None
    tc = profile.get("surface_temp_c")
    dc = profile.get("surface_dewpoint_c")
    if tc is not None and dc is not None:
        abs_hum = absolute_humidity(tc, dc)

    f_humidity = score_humidity(abs_hum, band_cfg)
    f_tod = score_time_of_day(contact["utc_hour"], contact["utc_minute"], month, contact.get("lon", -97.0))
    f_td = score_td_depression(tc, dc, band_cfg)
    f_ref = score_refractivity(profile.get("min_refractivity_gradient"), band_cfg)
    f_sky = score_sky(None)
    f_season = score_season(month, band_cfg)
    f_wind = score_wind(None)
    f_rain = score_rain(None, band_cfg)
    f_pwat = score_pwat(profile.get("pwat_mm"), band_cfg)
    f_pressure = score_pressure(profile.get("surface_pressure_mb"))

    raw_factors = {
        "humidity": f_humidity,
        "time_of_day": f_tod,
        "td_depression": f_td,
        "refractivity": f_ref,
        "sky": f_sky,
        "season": f_season,
        "wind": f_wind,
        "rain": f_rain,
        "pwat": f_pwat,
        "pressure": f_pressure,
    }
    ws = sum(raw_factors[k] * weights.get(k, 0.1) for k in FACTOR_ORDER)
    return round(ws)


# ── Statistics ─────────────────────────────────────────────────────────────────

def spearman_rho(xs: list[float], ys: list[float]) -> float:
    """Compute Spearman rank correlation coefficient."""
    n = len(xs)
    if n < 3:
        return 0.0

    def rank(vals: list[float]) -> list[float]:
        indexed = sorted(enumerate(vals), key=lambda x: x[1])
        ranks = [0.0] * len(vals)
        i = 0
        while i < len(indexed):
            j = i
            while j < len(indexed) and indexed[j][1] == indexed[i][1]:
                j += 1
            avg_rank = (i + j - 1) / 2.0 + 1.0
            for k in range(i, j):
                ranks[indexed[k][0]] = avg_rank
            i = j
        return ranks

    rx = rank(xs)
    ry = rank(ys)
    mean_rx = sum(rx) / n
    mean_ry = sum(ry) / n
    num = sum((rx[i] - mean_rx) * (ry[i] - mean_ry) for i in range(n))
    den_x = math.sqrt(sum((rx[i] - mean_rx) ** 2 for i in range(n)))
    den_y = math.sqrt(sum((ry[i] - mean_ry) ** 2 for i in range(n)))
    if den_x == 0 or den_y == 0:
        return 0.0
    return num / (den_x * den_y)


# ── Data loading ───────────────────────────────────────────────────────────────

# Snap constants: HRRR grid is 0.25° resolution in practice, but let's round
# to 0.125° (1/8) for snapping — matches the "nearest grid point" approach.
SNAP_RESOLUTION = 0.125

CONTACTS_SQL = """
WITH joined AS (
  SELECT DISTINCT ON (c.id)
    c.id,
    c.band::int AS band,
    c.distance_km::float AS dist,
    c.qso_timestamp,
    EXTRACT(HOUR FROM c.qso_timestamp)::int AS utc_hour,
    EXTRACT(MINUTE FROM c.qso_timestamp)::int AS utc_minute,
    EXTRACT(MONTH FROM c.qso_timestamp)::int AS month,
    EXTRACT(YEAR FROM c.qso_timestamp)::int AS year,
    (c.pos1->>'lon')::float AS lon,
    (c.pos1->>'lat')::float AS lat
  FROM contacts c
  JOIN hrrr_profiles h
    ON h.lat BETWEEN (c.pos1->>'lat')::float - 0.07
                 AND (c.pos1->>'lat')::float + 0.07
   AND h.lon BETWEEN (c.pos1->>'lon')::float - 0.07
                 AND (c.pos1->>'lon')::float + 0.07
   AND h.valid_time BETWEEN c.qso_timestamp - INTERVAL '1 hour'
                       AND c.qso_timestamp + INTERVAL '1 hour'
  WHERE c.pos1 IS NOT NULL
    AND c.distance_km BETWEEN 0 AND 3000
    AND c.flagged_invalid IS NOT TRUE
    AND c.qso_timestamp >= '2020-01-01'
  ORDER BY c.id, ABS(EXTRACT(EPOCH FROM h.valid_time - c.qso_timestamp))
)
SELECT * FROM joined
WHERE band >= 50
ORDER BY id
"""


def snap_coord(val: float) -> float:
    """Round a coordinate to the nearest HRRR grid multiple."""
    return round(val / SNAP_RESOLUTION) * SNAP_RESOLUTION


def load_contacts(conn) -> list[dict]:
    """Load contact → HRRR joined rows (same pattern as validate_algo.py)."""
    print("  executing contact↔HRRR join…", file=sys.stderr)
    with conn.cursor(row_factory=dict_row) as cur:
        cur.execute(CONTACTS_SQL)
        rows = cur.fetchall()
    print(f"  loaded {len(rows)} joined rows", file=sys.stderr)
    return rows


def load_lagged_profiles(conn, contacts: list[dict]) -> dict[tuple[int, int], dict]:
    """For each lag, load one HRRR profile per contact at the snapped location.

    Returns a dict keyed by (contact_id, lead_hours) → profile row.
    Contacts are snapped to the nearest 0.125° grid point.
    """
    # First, build the set of unique (lat, lon, valid_time) tuples we need
    # Group by contact: we need one profile per (contact, lead) pair
    by_snapped = defaultdict(list)
    for c in contacts:
        slat = snap_coord(c["lat"])
        slon = snap_coord(c["lon"])
        ts = c["qso_timestamp"]
        for lead in LEAD_TIMES:
            by_snapped[(slat, slon, lead)].append((c["id"], ts, lead))

    print(f"  fetching {len(by_snapped)} unique (lat, lon, lead) combinations…", file=sys.stderr)

    # Batch fetch: one query per lead time for efficiency
    result: dict[tuple[int, int], dict] = {}  # (contact_id, lead_hours) → profile

    for lead in LEAD_TIMES:
        # Get all (lat, lon) pairs for this lead
        lats = sorted({k[0] for k in by_snapped if k[2] == lead})
        lons = sorted({k[1] for k in by_snapped if k[2] == lead})
        if not lats:
            continue

        # Use ANY for efficient batch lookup when many points
        sql = """
        SELECT lat, lon, valid_time,
               surface_temp_c, surface_dewpoint_c,
               surface_pressure_mb, pwat_mm, hpbl_m,
               min_refractivity_gradient
        FROM hrrr_profiles
        WHERE lat = ANY(%s)
          AND lon = ANY(%s)
          AND valid_time = ANY(
            SELECT date_trunc('hour', qso_timestamp) - make_interval(hours => %s)
            FROM contacts
            WHERE id = ANY(%s)
              AND qso_timestamp >= '2020-01-01'
          )
        """

        # Actually, the ANY subquery approach is complex. Let's do it differently:
        # For each contact, compute valid_time and do an IN clause.
        # Better: precompute all valid_times and batch-query by (lat, lon, valid_time) tuples.
        # But we can't pass tuples to ANY easily. Let's use a VALUES join.

        # Simpler approach: one query that joins contacts to hrrr_profiles for each lag
        contact_ids = [c["id"] for c in contacts]

        # Build a VALUES table approach for efficiency
        lag_sql = f"""
        SELECT p.surface_temp_c, p.surface_dewpoint_c,
               p.surface_pressure_mb, p.pwat_mm, p.hpbl_m,
               p.min_refractivity_gradient,
               c.id AS contact_id, {lead} AS lead_hours
        FROM contacts c
        CROSS JOIN LATERAL (
          SELECT *
          FROM hrrr_profiles h
          WHERE h.lat = {snap_coord_fn_sql('(c.pos1->>''lat'')::float')}
            AND h.lon = {snap_coord_fn_sql('(c.pos1->>''lon'')::float')}
            AND h.valid_time = date_trunc('hour', c.qso_timestamp) - make_interval(hours => {lead})
          LIMIT 1
        ) p
        WHERE c.id = ANY(%s)
          AND c.qso_timestamp >= '2020-01-01'
        """

        # Actually the CROSS JOIN LATERAL with the snap function in SQL is messy.
        # Let me take a simpler but effective approach: use a VALUES clause to batch.

        # Build the set of (lat, lon, valid_times) to query
        # Since the number of contacts might be large, use a temp table or VALUES approach

        # Alternative: For each unique (snapped_lat, snapped_lon), query all lags at once
        pass  # placeholder — see actual implementation below

    return result


def snap_coord_fn_sql(formula: str) -> str:
    """Generate SQL for rounding a formula to the nearest SNAP_RESOLUTION multiple."""
    r = SNAP_RESOLUTION
    return f"(round(({formula}) / {r}) * {r})"


def fetch_all_lagged_profiles(conn, contacts: list[dict]) -> dict[tuple[str, int], dict]:
    """Fetch all lagged HRRR profiles efficiently.

    Strategy: Uses a temp table with all (snapped_lat, snapped_lon, lead, base_ts)
    tuples, then joins against hrrr_profiles in a single query.

    Returns: {(contact_id_str, lead_hours): profile_row} where profile_row is a dict
    with the weather fields.
    """
    if not contacts:
        return {}

    print(f"  fetching lagged profiles for {len(contacts)} contacts at leads {LEAD_TIMES}…", file=sys.stderr)

    # Build lookup rows: (contact_id, snapped_lat, snapped_lon, lead_hours, base_ts)
    # We generate this in Python and pass to PostgreSQL via a VALUES clause

    rows = []
    for c in contacts:
        slat = snap_coord(c["lat"])
        slon = snap_coord(c["lon"])
        base_ts = c["qso_timestamp"]
        for lead in LEAD_TIMES:
            rows.append((c["id"], slat, slon, lead, base_ts))

    if not rows:
        return {}

    # Use a temp table to batch the profile lookups efficiently.
    # Without ON COMMIT DROP since autocommit=True.
    with conn.cursor() as cur:
        cur.execute("DROP TABLE IF EXISTS _forecast_lookups")
        cur.execute("""
        CREATE TEMP TABLE _forecast_lookups (
            contact_id text,
            snapped_lat double precision,
            snapped_lon double precision,
            lead_hours int,
            base_ts timestamptz
        )
        """)
        # Insert in batches of 5000 rows to stay under param limits
        batch_size = 5000
        for i in range(0, len(rows), batch_size):
            batch = rows[i:i + batch_size]
            values_clauses = []
            params = []
            for contact_id, slat, slon, lead, base_ts in batch:
                cid_str = str(contact_id)
                values_clauses.append("(%s::text, %s::double precision, %s::double precision, %s::int, %s::timestamptz)")
                params.extend([cid_str, slat, slon, lead, base_ts])
            sql = "INSERT INTO _forecast_lookups (contact_id, snapped_lat, snapped_lon, lead_hours, base_ts) VALUES " + ", ".join(values_clauses)
            cur.execute(sql, params)

    # Now query hrrr_profiles joining against the lookup table
    print(f"  querying hrrr_profiles with {len(rows)} lookup tuples…", file=sys.stderr)
    with conn.cursor(row_factory=dict_row) as cur:
        cur.execute("""
        SELECT
            fl.contact_id,
            fl.lead_hours,
            h.surface_temp_c,
            h.surface_dewpoint_c,
            h.surface_pressure_mb,
            h.pwat_mm,
            h.hpbl_m,
            h.min_refractivity_gradient
        FROM _forecast_lookups fl
        JOIN hrrr_profiles h
          ON h.lat = fl.snapped_lat
         AND h.lon = fl.snapped_lon
         AND h.valid_time = date_trunc('hour', fl.base_ts) - make_interval(hours => fl.lead_hours)
        """)
        result_rows = cur.fetchall()

    # Cleanup
    with conn.cursor() as cur:
        cur.execute("DROP TABLE IF EXISTS _forecast_lookups")

    result: dict[tuple[str, int], dict] = {}
    for r in result_rows:
        key = (r["contact_id"], r["lead_hours"])
        result[key] = {
            "surface_temp_c": r["surface_temp_c"],
            "surface_dewpoint_c": r["surface_dewpoint_c"],
            "surface_pressure_mb": r["surface_pressure_mb"],
            "pwat_mm": r["pwat_mm"],
            "hpbl_m": r["hpbl_m"],
            "min_refractivity_gradient": r["min_refractivity_gradient"],
        }

    print(f"  fetched {len(result)} profile rows", file=sys.stderr)

    # Report coverage per lead
    for lead in LEAD_TIMES:
        found = sum(1 for (cid, lh) in result if lh == lead)
        print(f"    lead-{lead}h: {found}/{len(contacts)} profiles found", file=sys.stderr)

    return result


# ── Main ───────────────────────────────────────────────────────────────────────

def main() -> int:
    dsn = os.environ.get("PROP_PROD_DB_URL")
    if not dsn:
        print("error: PROP_PROD_DB_URL not set", file=sys.stderr)
        return 2

    today = dt.date.today().isoformat()
    report_dir = REPORT_DIR
    report_dir.mkdir(parents=True, exist_ok=True)

    # ── Load weights ────────────────────────────────────────────────────────
    print("• loading band weights", file=sys.stderr)
    if not WEIGHTS_PATH.exists():
        print(f"  warning: {WEIGHTS_PATH} not found — using global defaults", file=sys.stderr)
        weights_json = {"global_weights": {
            "humidity": 0.1262, "time_of_day": 0.038, "td_depression": 0.101,
            "refractivity": 0.0986, "sky": 0.0841, "season": 0.1134,
            "wind": 0.0841, "rain": 0.1431, "pwat": 0.1147, "pressure": 0.0967,
        }, "band_overrides": {}}
    else:
        with open(WEIGHTS_PATH) as f:
            weights_json = json.load(f)

    global_weights: dict[str, float] = weights_json["global_weights"]
    band_overrides: dict[str, dict] = weights_json.get("band_overrides", {})

    def get_weights(band_mhz: int) -> dict[str, float]:
        key = str(band_mhz)
        if key in band_overrides and "weights" in band_overrides[key]:
            return band_overrides[key]["weights"]
        return dict(global_weights)

    # ── Load data ───────────────────────────────────────────────────────────
    print("• connecting to database", file=sys.stderr)
    with psycopg.connect(dsn, autocommit=True) as conn:
        with conn.cursor() as cur:
            cur.execute("SET statement_timeout = '30min'")
        contacts = load_contacts(conn)

        if not contacts:
            print("error: no contacts loaded — check DB connection and data", file=sys.stderr)
            return 3

        print(f"• loaded {len(contacts)} contacts (≥2020-01-01)", file=sys.stderr)

        # ── Fetch lagged profiles ──────────────────────────────────────────
        print("• fetching lagged HRRR profiles", file=sys.stderr)
        all_profiles = fetch_all_lagged_profiles(conn, contacts)

    # ── Score each contact at each lag ──────────────────────────────────────
    print("• scoring contacts at each lead time", file=sys.stderr)
    # Track which contacts have ALL required profiles (lag-0 through lag-24)
    scored_by_band: dict[int, dict[int, list]] = defaultdict(lambda: defaultdict(list))
    # {band_mhz: {lead_hours: [(score, distance), ...]}}

    contacts_with_all = 0
    for c in contacts:
        band_mhz = c["band"]
        band_cfg = BAND_CONFIGS.get(band_mhz)
        if band_cfg is None:
            continue

        w = get_weights(band_mhz)
        cid_str = str(c["id"])

        # Check if all profiles exist
        has_all = True
        scores_at_lag = {}
        for lead in LEAD_TIMES:
            key = (cid_str, lead)
            profile = all_profiles.get(key)
            if profile is None:
                has_all = False
                break
            score = compute_composite_score(profile, c, band_cfg, w)
            if score is None:
                has_all = False
                break
            scores_at_lag[lead] = (score, c["dist"])

        if has_all:
            contacts_with_all += 1
            for lead, (score, dist) in scores_at_lag.items():
                scored_by_band[band_mhz][lead].append((score, dist))

    print(f"• contacts with all {len(LEAD_TIMES)} profiles: {contacts_with_all}/{len(contacts)}", file=sys.stderr)

    # ── Per-band per-lag Spearman ρ ────────────────────────────────────────
    print("• computing per-band per-lag Spearman ρ", file=sys.stderr)
    per_band_results = {}
    for band_mhz in sorted(scored_by_band.keys()):
        by_lead = scored_by_band[band_mhz]
        # All leads must have the same N (contacts with all profiles)
        ns = [len(pairs) for pairs in by_lead.values()]
        n = ns[0] if ns else 0

        if n < 30:
            print(f"  band {band_mhz} MHz: {n} contacts — skipping (<30 threshold)", file=sys.stderr)
            continue
        if not all(x == n for x in ns):
            print(f"  band {band_mhz} MHz: inconsistent counts — skipping", file=sys.stderr)
            continue

        rhos = {}
        for lead in LEAD_TIMES:
            pairs = by_lead.get(lead, [])
            scores = [s for s, d in pairs]
            dists = [d for s, d in pairs]
            rhos[lead] = spearman_rho(scores, dists)

        rho_0 = rhos.get(0, 0.0)
        per_band_results[band_mhz] = {
            "n": n,
            "rhos": rhos,
            "delta_6h": round(rho_0 - rhos.get(6, rho_0), 4),
            "delta_24h": round(rho_0 - rhos.get(24, rho_0), 4),
        }
        print(f"  band {band_mhz} MHz: n={n}, ρ(0h)={rhos.get(0, 0):+.4f}, "
              f"ρ(6h)={rhos.get(6, 0):+.4f}, ρ(24h)={rhos.get(24, 0):+.4f}", file=sys.stderr)

    # ── Build report ────────────────────────────────────────────────────────
    generated_at = dt.datetime.now(dt.timezone.utc).replace(microsecond=0)

    json_payload = {
        "generated_at": generated_at.isoformat(timespec="seconds"),
        "generated_by": "scripts/validate_forecast.py",
        "schema_version": 1,
        "summary": {
            "total_contacts": len(contacts),
            "contacts_with_all_profiles": contacts_with_all,
            "lead_times_hours": LEAD_TIMES,
            "bands_evaluated": len(per_band_results),
        },
        "per_band": {
            str(b): {
                "n": info["n"],
                "rhos": {str(k): round(v, 4) for k, v in info["rhos"].items()},
                "delta_6h": info["delta_6h"],
                "delta_24h": info["delta_24h"],
            }
            for b, info in sorted(per_band_results.items())
        },
    }

    json_path = report_dir / f"forecast-{today}.json"
    json_path.write_text(json.dumps(json_payload, indent=2, sort_keys=False) + "\n")
    print(f"• wrote {json_path}", file=sys.stderr)

    # ── Render Markdown table ───────────────────────────────────────────────
    md = []
    md.append(f"# Forecast Skill Degradation — {today}\n")
    md.append(
        "> Auto-generated by `scripts/validate_forecast.py`. Measures how well "
        "the propagation scoring algorithm predicts contact distances as the "
        "weather profile ages. The lag-0h profile uses actual conditions at "
        "contact time (ground truth). Lag-Nh uses the weather profile from "
        "N hours earlier — representing the forecast that would have been "
        "available at that lead time.\n"
    )
    md.append(
        f"- **Total contacts** (≥2020): {len(contacts):,}"
    )
    md.append(
        f"- **Contacts with all {len(LEAD_TIMES)} profiles**: {contacts_with_all:,} " +
        f"({contacts_with_all * 100 // max(1, len(contacts))}%)"
    )
    md.append(f"- **Generated**: {generated_at.isoformat(timespec='seconds')}\n")

    md.append(
        "Higher ρ = better distance prediction. Δρ(0h→Nh) = ρ(0h) − ρ(Nh) "
        "measures skill degradation — positive values mean the forecast skill "
        "is worse than using current conditions.\n"
    )

    md.append("## Per-Band Per-Lag Spearman ρ\n")

    # Build header
    col_headers = "| Band | N |"
    col_sep = "|------|--:|"
    for lead in LEAD_TIMES:
        col_headers += f" ρ({lead}h) |"
        col_sep += "-------:|" if lead == 0 else "------:|"
    col_headers += " Δρ(0→6h) | Δρ(0→24h) |"
    col_sep += "----------|------------|"

    md.append(col_headers)
    md.append(col_sep)

    # Sort bands by contact count descending for display
    sorted_bands = sorted(per_band_results.keys(), key=lambda b: per_band_results[b]["n"], reverse=True)
    for band_mhz in sorted_bands:
        info = per_band_results[band_mhz]
        row = f"| {band_mhz} MHz | {info['n']} |"
        for lead in LEAD_TIMES:
            rho = info["rhos"].get(lead, 0.0)
            row += f" {rho:+.4f} |"
        row += f" {info['delta_6h']:+.4f} | {info['delta_24h']:+.4f} |"
        md.append(row)

    # ── Key findings ────────────────────────────────────────────────────────
    md.append("")
    md.append("## Key Findings\n")

    # Check monotonic degradation
    monotonic_bands = 0
    non_monotonic_bands = []
    for band_mhz in sorted_bands:
        info = per_band_results[band_mhz]
        rhos = info["rhos"]
        prev = None
        monotonic = True
        for lead in LEAD_TIMES:
            rho = rhos.get(lead, 0.0)
            if prev is not None and rho > prev + 0.001:  # slight tolerance
                monotonic = False
            prev = rho
        if monotonic:
            monotonic_bands += 1
        else:
            non_monotonic_bands.append(band_mhz)

    md.append(f"- **Monotonic degradation** (ρ decreases with lead time): {monotonic_bands}/{len(sorted_bands)} bands")

    if non_monotonic_bands:
        md.append(f"- **Non-monotonic bands** (ρ increases at some lead): {', '.join(f'{b} MHz' for b in non_monotonic_bands)}")
        md.append(
            "  - This may indicate the algorithm puts disproportionate weight "
            "on noisy short-term features that vary within hours. If ρ(6h) > ρ(0h), "
            "the current-conditions score is noisier than the 6h-lagged score — "
            "weather features at lag-0 may introduce variance that degrades the "
            "correlation."
        )

    # @6h vs @0h comparison
    rho_6h_better = 0
    for band_mhz in sorted_bands:
        info = per_band_results[band_mhz]
        r0 = info["rhos"].get(0, 0.0)
        r6 = info["rhos"].get(6, 0.0)
        if r6 > r0:
            rho_6h_better += 1
    if rho_6h_better > 0:
        md.append(f"- **ρ(6h) > ρ(0h)**: {rho_6h_better} bands — the 6-hour lagged score outperforms current-conditions")

    # Average degradation
    if per_band_results:
        avg_delta_6 = sum(info["delta_6h"] for info in per_band_results.values()) / len(per_band_results)
        avg_delta_24 = sum(info["delta_24h"] for info in per_band_results.values()) / len(per_band_results)
        md.append(f"- **Mean Δρ(0→6h)**: {avg_delta_6:+.4f}")
        md.append(f"- **Mean Δρ(0→24h)**: {avg_delta_24:+.4f}")
    else:
        md.append("- **No bands passed the ≥30 contact threshold.**")

    # Caveats
    md.append("")
    md.append("## Caveats\n")
    md.append(
        "- **Lagged profiles use analysis data**: The HRRR profiles at T−N "
        "represent actual weather conditions at that earlier time, not a true "
        "forecast run initialized at T−N. Real HRRR forecasts would contain "
        "model error growth; these measurements capture the degradation from "
        "weather *evolution* alone (the lower bound of forecast skill loss)."
    )
    md.append(
        "- **Snap to nearest grid point**: Profiles are looked up at the HRRR "
        "grid point nearest to the contact location (rounded to 0.125°). "
        "A single point approximates the full path's conditions."
    )
    md.append(
        "- **Filter threshold**: Only bands with ≥30 contacts having ALL "
        "6 profiles (lag-0 through lag-24) are included in the table."
    )
    md.append(
        "- **Censored data**: Contacts are confirmed success events only. "
        "Correlation ρ measures discrimination among successful contacts."
    )
    md.append("")

    md_path = report_dir / f"forecast-{today}.md"
    md_text = "\n".join(md)
    md_path.write_text(md_text + "\n")
    print(f"• wrote {md_path}", file=sys.stderr)

    # ── Print to stdout ─────────────────────────────────────────────────────
    print(md_text)

    return 0


if __name__ == "__main__":
    raise SystemExit(main())