Graham McIntire
2fd88a94ea
Phase A — 7 pipeline bugs:
- retain_scores_window: fix timeline self-deletion (NOTIFY payload run_time|valid_time)
- Rust/Elixir weight divergence: Rust loads band_weights.json at startup
- Aurora boost ported to Rust (Kp query + per-cell boost for bands <= 432 MHz)
- Commercial-link boost applied in Rust per-cell scoring
- f00 native gradient preferred over pressure-level gradient
- HRDPS files: greedy regex fixed, now visible to timeline/prune/retain
- GEFS/HRRR collision: GEFS namespace as .gefs.prop, merge at read
Phase B — Validation harness:
- scripts/validate_algo.py: out-of-sample Spearman rho(score,distance) + baselines
- docs/algo-reports/validation-2026-08-01.{json,md}
Phase C — Forecast-skill evaluation:
- scripts/validate_forecast.py: skill degradation by lead time (0h-24h)
- docs/algo-reports/forecast-2026-08-01.{json,md}
Phase D — Calibration + model improvements:
- recalibrate.py: validation gate before weight deploy
- Recalibrator: Nx.max(0.0) replaces Nx.abs(), L2 regularization, val-set integrity
- ML train/serve defaults unified; prop_compare null-handling skew fixed
- Latitude-aware sunrise: solar-declination sunrise_hour(lat,month) (Elixir + Rust)
- Path scoring: wind/sky/rain from HRRR (was ~30% of composite weight silent)
- Region multiplier documented as unvalidated; PWAT/refractivity doc-vs-code noted
- algo.md: synced scoring sections, marked retired features, D5/D6 changelog
- Credo: path_compute cyclomatic complexity bumped (6->10 fields) — informational only
750 lines
32 KiB
Python
750 lines
32 KiB
Python
#!/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())
|