#!/usr/bin/env python3
"""
Algorithm Skill Validation — measures how accurately the propagation
scoring algorithm predicts contact distances on held-out data.

This is the FIRST rigorous out-of-sample validation. The existing
recalibration pipeline (`scripts/recalibrate.py`) fits per-band weights
on the full corpus with zero holdout. This script:

  1. Temporally splits contacts at 2025-01-01 (fit < → test ≥)
  2. Computes a simplified composite score using all 10 factors
  3. Reports per-band Spearman-ρ(score, distance) on the test set
  4. Produces calibration curves (binned deciles vs median/P90 distance)
  5. Compares algorithm skill vs persistence, climatology, and no-skill baselines

Output:
  - Markdown table to stdout
  - docs/algo-reports/validation-YYYY-MM-DD.json
  - docs/algo-reports/validation-YYYY-MM-DD.md

Usage:
    python3 scripts/validate_algo.py
"""

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 lib/microwaveprop/propagation/band_config.ex) ──

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

# Per-band configs: freq_mhz → {humidity_effect, humidity_penalty, rain_k, rain_alpha, seasonal_base, seasonal_adj}
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",
]

# ── Scoring functions (mirrors lib/microwaveprop/propagation/scorer.ex) ──

def absolute_humidity(temp_c: float, dewpoint_c: float) -> float:
    """Absolute humidity in g/m³ from temperature and dewpoint (both °C)."""
    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:  # harmful
        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:  # harmful
        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:  # harmful
        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 composite_score(contact: dict, band_cfg: dict, weights: dict[str, float]) -> dict:
    """Compute all 10 factor scores and the weighted composite for one contact."""
    month = contact["month"]
    f_humidity = score_humidity(contact.get("abs_humidity"), 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(contact.get("surface_temp_c"), contact.get("surface_dewpoint_c"), band_cfg)
    f_ref = score_refractivity(contact.get("min_refractivity_gradient"), band_cfg)
    f_sky = score_sky(None)  # No cloud cover data in HRRR profiles
    f_season = score_season(month, band_cfg)
    f_wind = score_wind(None)  # No wind data in HRRR profiles
    f_rain = score_rain(None, band_cfg)  # No precip data in HRRR profiles
    f_pwat = score_pwat(contact.get("pwat_mm"), band_cfg)
    f_pressure = score_pressure(contact.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 {"score": round(ws), "factors": raw_factors}


# ── 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
    # Rank the values
    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  # 1-based ranks
            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)


def percentile(sorted_vals: list[float], pct: float) -> float:
    """Compute p-th percentile from sorted values."""
    n = len(sorted_vals)
    if n == 0:
        return 0.0
    idx = (pct / 100.0) * (n - 1)
    lo = int(math.floor(idx))
    hi = int(math.ceil(idx))
    if lo == hi:
        return sorted_vals[lo]
    frac = idx - lo
    return sorted_vals[lo] * (1 - frac) + sorted_vals[hi] * frac


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

CONTACTS_HRRR_SQL = """
WITH joined AS (
  SELECT DISTINCT ON (c.id)
    c.id, c.band::int AS band, c.distance_km::float AS dist,
    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,
    h.surface_temp_c, h.surface_dewpoint_c,
    h.surface_pressure_mb, h.pwat_mm, h.hpbl_m,
    h.min_refractivity_gradient
  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
  ORDER BY c.id, ABS(EXTRACT(EPOCH FROM h.valid_time - c.qso_timestamp))
)
SELECT * FROM joined
WHERE band >= 50
ORDER BY id
"""


def load_contacts(conn, limit: int | None = None) -> list[dict]:
    """Load contact ↔ HRRR joined rows."""
    sql = CONTACTS_HRRR_SQL
    if limit:
        sql += f" LIMIT {limit}"
    print(f"  executing contact↔HRRR join...", file=sys.stderr)
    with conn.cursor(row_factory=dict_row) as cur:
        cur.execute(sql)
        rows = cur.fetchall()
    print(f"  loaded {len(rows)} joined rows", file=sys.stderr)
    return rows


# ── Baselines ─────────────────────────────────────────────────────────────────

def compute_persistence_baseline(contacts: list[dict], band_mhz: int) -> dict[tuple[int, int], float]:
    """Fit set: per (band, month) median distance."""
    by_month = defaultdict(list)
    for c in contacts:
        if c["band"] == band_mhz:
            by_month[(c["band"], c["month"])].append(c["dist"])
    return {k: sorted(v)[len(v) // 2] for k, v in by_month.items()}


def compute_climatology_baseline(contacts: list[dict]) -> dict[int, float]:
    """Fit set: per-band global median distance."""
    by_band = defaultdict(list)
    for c in contacts:
        by_band[c["band"]].append(c["dist"])
    return {b: sorted(v)[len(v) // 2] for b, v in by_band.items()}


# ── 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'")
        all_contacts = load_contacts(conn)

    # ── Split: fit (< 2025-01-01) vs test (≥ 2025-01-01) ──────────────────
    fit_contacts = [c for c in all_contacts if c["year"] < 2025]
    test_contacts = [c for c in all_contacts if c["year"] >= 2025]
    print(f"• fit set: {len(fit_contacts)} contacts (pre-2025)", file=sys.stderr)
    print(f"• test set: {len(test_contacts)} contacts (2025+)", file=sys.stderr)

    if len(fit_contacts) < 50:
        print("  WARNING: fit set < 50 contacts — baselines unreliable", file=sys.stderr)
    if len(test_contacts) < 50:
        print("  WARNING: test set < 50 contacts — validation unreliable", file=sys.stderr)

    # ── Build baselines from fit set ───────────────────────────────────────
    print("• computing baselines from fit set", file=sys.stderr)
    # Climatology (per-band global median)
    climatology = compute_climatology_baseline(fit_contacts)
    # Persistence (per-band per-month median)
    persistence = {}
    all_bands = sorted(set(c["band"] for c in all_contacts))
    for b in all_bands:
        persistence.update(compute_persistence_baseline(fit_contacts, b))
    # No-skill (overall median)
    all_dists = sorted([c["dist"] for c in fit_contacts])
    no_skill_median = all_dists[len(all_dists) // 2] if all_dists else 0.0

    # ── Score test contacts ────────────────────────────────────────────────
    print("• scoring test contacts", file=sys.stderr)
    for c in test_contacts:
        band_mhz = c["band"]
        band_cfg = BAND_CONFIGS.get(band_mhz)
        if band_cfg is None:
            c["composite_score"] = None
            continue
        # Compute derived fields
        tc = c.get("surface_temp_c")
        dc = c.get("surface_dewpoint_c")
        if tc is not None and dc is not None:
            c["abs_humidity"] = absolute_humidity(tc, dc)
        else:
            c["abs_humidity"] = None

        w = get_weights(band_mhz)
        result = composite_score(c, band_cfg, w)
        c["composite_score"] = result["score"]
        c["factor_scores"] = result["factors"]

    # ── Per-band analysis ──────────────────────────────────────────────────
    print("• computing per-band statistics", file=sys.stderr)

    per_band = {}
    for band_mhz in sorted(BAND_CONFIGS.keys()):
        band_tests = [c for c in test_contacts if c["band"] == band_mhz and c.get("composite_score") is not None]
        n = len(band_tests)
        if n == 0:
            continue

        scores = [c["composite_score"] for c in band_tests]
        dists = [c["dist"] for c in band_tests]

        # Algorithm ρ
        alg_rho = spearman_rho(scores, dists)

        # Baseline ρ
        pers_dists = []
        clim_dists = []
        ns_dists = []
        for c in band_tests:
            # Persistence: per (band, month) median from fit set
            p_med = persistence.get((band_mhz, c["month"]), no_skill_median)
            pers_dists.append(p_med)
            # Climatology: per-band global median from fit set
            c_med = climatology.get(band_mhz, no_skill_median)
            clim_dists.append(c_med)
            # No-skill
            ns_dists.append(no_skill_median)

        pers_rho = spearman_rho(pers_dists, dists)
        clim_rho = spearman_rho(clim_dists, dists)
        ns_rho = spearman_rho(ns_dists, dists)

        # Skill gain
        skill_gain = alg_rho - pers_rho

        # Calibration: decile bins
        deciles = defaultdict(list)
        for c in band_tests:
            decile = min(9, c["composite_score"] // 10)
            deciles[decile].append(c["dist"])

        calibration = []
        for d in range(10):
            if deciles[d]:
                sd = sorted(deciles[d])
                median = sd[len(sd) // 2]
                p90 = percentile(sd, 90)
                calibration.append({
                    "decile": d,
                    "score_range": f"{d*10}-{min(100, (d+1)*10)}",
                    "n": len(sd),
                    "median_km": round(median, 1),
                    "p90_km": round(p90, 1),
                })

        per_band[band_mhz] = {
            "n": n,
            "alg_rho": round(alg_rho, 4),
            "pers_rho": round(pers_rho, 4),
            "clim_rho": round(clim_rho, 4),
            "no_skill_rho": round(ns_rho, 4),
            "skill_gain": round(skill_gain, 4),
            "calibration": calibration,
            "score_min": min(scores),
            "score_max": max(scores),
            "score_mean": round(sum(scores) / n, 1),
            "score_median": sorted(scores)[n // 2],
            "dist_median_km": round(sorted(dists)[n // 2], 1),
            "dist_p90_km": round(percentile(sorted(dists), 90), 1),
        }

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

    json_payload = {
        "generated_at": generated_at.isoformat(timespec="seconds"),
        "generated_by": "scripts/validate_algo.py",
        "schema_version": 1,
        "summary": {
            "total_contacts": len(all_contacts),
            "fit_contacts": len(fit_contacts),
            "test_contacts": len(test_contacts),
            "bands_with_50_test": sum(1 for v in per_band.values() if v["n"] >= 50),
            "total_bands": len(per_band),
        },
        "per_band": {str(k): v for k, v in sorted(per_band.items())},
    }

    # ── Write JSON ─────────────────────────────────────────────────────────
    json_path = report_dir / f"validation-{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 ────────────────────────────────────────────────────
    md = []
    md.append(f"# Algorithm Skill Validation — {today}\n")
    md.append(
        "> Auto-generated by `scripts/validate_algo.py`. Measures how "
        "accurately the propagation scoring algorithm predicts contact "
        "distances on held-out data (contacts from 2025-01-01 onward).\n"
    )
    md.append(
        f"- **Total contacts**: {len(all_contacts):,} "
        f"(fit: {len(fit_contacts):,}, test: {len(test_contacts):,})"
    )
    md.append(f"- **Generated**: {generated_at.isoformat(timespec='seconds')}\n")

    # Summary table
    md.append("## Per-Band Spearman ρ\n")
    md.append(
        "Higher ρ = the algorithm score better predicts contact distance. "
        "Positive skill gain means the algorithm beats the persistence "
        "baseline (per-band per-month median distance).\n"
    )
    md.append("| Band | n | ρ(alg) | ρ(pers) | ρ(clim) | ρ(no-skill) | Skill Gain |")
    md.append("|------|--:|-------:|--------:|--------:|------------:|-----------:|")

    key_bands = [222, 432, 902, 1296, 2304, 5760, 10000, 24000]
    for band_mhz in sorted(per_band.keys()):
        info = per_band[band_mhz]
        if info["n"] < 50 and band_mhz not in key_bands:
            continue
        md.append(
            f"| {band_mhz} MHz | {info['n']} | {info['alg_rho']:+.4f} | "
            f"{info['pers_rho']:+.4f} | {info['clim_rho']:+.4f} | "
            f"{info['no_skill_rho']:+.4f} | {info['skill_gain']:+.4f} |"
        )

    if any(info["n"] >= 50 for info in per_band.values()):
        md.append(f"\n_Bands with <50 test contacts omitted from table._\n")
    else:
        md.append(f"\n_No band has ≥50 test contacts._\n")

    # Calibration curves for top 4 bands by contact count
    top_bands = sorted(per_band.keys(), key=lambda b: per_band[b]["n"], reverse=True)[:4]
    if top_bands:
        md.append("## Calibration Curves (Top 4 Bands)\n")
        md.append(
            "A well-calibrated scorer should show monotonically increasing "
            "median and P90 distances as the score decile increases. "
            "Flat or inverted curves indicate the score is not capturing "
            "distance information.\n"
        )
        for band_mhz in top_bands:
            info = per_band[band_mhz]
            cal = info["calibration"]
            if not cal or len(cal) < 3:
                continue
            md.append(f"### {band_mhz} MHz (n={info['n']})\n")
            md.append("| Decile | Score Range | n | Median km | P90 km |")
            md.append("|--------|------------|--:|----------:|-------:|")
            for row in cal:
                md.append(
                    f"| {row['decile']} | {row['score_range']} | {row['n']} | "
                    f"{row['median_km']:.1f} | {row['p90_km']:.1f} |"
                )
            md.append("")

    # Key findings
    md.append("## Key Findings\n")
    positive = sum(1 for info in per_band.values() if info["skill_gain"] > 0 and info["n"] >= 50)
    negative = sum(1 for info in per_band.values() if info["skill_gain"] <= 0 and info["n"] >= 50)
    md.append(f"- **Positive skill gain** (algorithm beats persistence): {positive} bands")
    md.append(f"- **Negative/zero skill gain**: {negative} bands")
    md.append(f"- **Total bands with ≥50 test contacts**: {positive + negative}")

    # Score vs distance summary
    overall_alg_rhos = [info["alg_rho"] for info in per_band.values() if info["n"] >= 50]
    if overall_alg_rhos:
        avg_alg_rho = sum(overall_alg_rhos) / len(overall_alg_rhos)
        md.append(f"- **Mean ρ(alg) across bands**: {avg_alg_rho:+.4f}")
    overall_gains = [info["skill_gain"] for info in per_band.values() if info["n"] >= 50]
    if overall_gains:
        avg_gain = sum(overall_gains) / len(overall_gains)
        md.append(f"- **Mean skill gain**: {avg_gain:+.4f}")

    # Caveats
    md.append("")
    md.append("## Caveats\n")
    md.append(
        "- **Simplified scorer**: Sky, wind, and rain factors use default "
        "values (50/50/100) because HRRR `cloud_cover_pct`, wind, and "
        "precip fields are not available in the production DB schema. "
        "The composite score reflects 7 of 10 factors; correlation values "
        "are a lower bound on what the full production scorer would achieve."
    )
    md.append(
        "- **Nearest-neighbour join**: Contacts are matched to the HRRR "
        "grid point within 0.07° (≈7 km) and ±1 hour. A single grid point "
        "approximates the full path's conditions — the production scorer "
        "uses multi-point path integration."
    )
    md.append(
        "- **Test set size**: 2025+ contacts (n={:,}) may be insufficient "
        "for mm-wave bands (47+ GHz). Treat those results as indicative.".format(
            len(test_contacts)
        )
    )
    md.append(
        "- **Censored data**: Contacts are confirmed success events only. "
        "The algorithm scores propagation quality on the path that was "
        "actually used. There is no counterfactual (paths where propagation "
        "was poor and no contact occurred). Correlation ρ measures "
        "discrimination among successful contacts, not an AUC over all "
        "possible paths."
    )
    md.append("")

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

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

    return 0


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