prop/scripts/recalibrate_algo.py

#!/usr/bin/env python3
"""
Algo recalibration analysis runner.

Connects to the prod (or any) Postgres, joins contacts to the nearest HRRR
grid profile within ±0.07° / ±1h (matching `Weather.find_nearest_hrrr/3`),
falls back to `narr_profiles` (NCEP NARR, 32 km / 3-hourly) for pre-2014
contacts where HRRR doesn't reach, and produces a Markdown report covering:

  * row counts and date ranges for the key tables (HRRR + NARR + ...)
  * per-band contact distribution from 50 MHz up (VHF/UHF/microwave)
  * monthly sounding ducting refresh
  * native HRRR profile duct stats
  * per-band Pearson correlations of contact distance vs HRRR fields
  * per-band Pearson correlations vs NARR fields (historical sanity check)
  * per-band binned distance distributions for HPBL, pressure, NEXRAD,
    and native-profile best duct band — every band with ≥50 matched
    contacts gets its own bin table so we're not hard-coding the model
    to 10/24 GHz behaviour
  * NEXRAD composite reflectivity vs contact distance (rain-scoring check)
  * hrrr_native_profiles.best_duct_band_ghz as a distance discriminator
  * commercial-link rx_power degradation vs contemporaneous DFW contacts

Output is saved to docs/algo-reports/YYYY-MM-DD-recalibration.md so it can
be diffed against algo.md or committed alongside it.

Connection string is read from PROP_PROD_DB_URL (set in .envrc) or the
positional --dsn argument. Dependency-light: psycopg + pandas only.

Usage:
    direnv allow            # picks up PROP_PROD_DB_URL from .envrc
    python3 scripts/recalibrate_algo.py
    # or override:
    python3 scripts/recalibrate_algo.py --dsn 'postgres://user:pw@host/db' --out path.md
"""

from __future__ import annotations

import argparse
import datetime as dt
import os
import sys
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' pandas"
    )

try:
    import pandas as pd
except ImportError:
    sys.exit("pandas is required: pip install pandas")


REPORT_HEADER = """\
# Recalibration Analysis Report — {today}

> Auto-generated by `scripts/recalibrate_algo.py`. This report is the
> input for `algo.md` updates — diff against the prior dated section to
> see what's actually moved.

Connection: `{dsn_redacted}`
Statement timeout: {stmt_timeout}
"""


# ─── SQL fragments ──────────────────────────────────────────────────────────

# (table, time_column). `time_column = None` means just count rows without
# min/max range. Tables not present in the target DB are silently skipped —
# local dev often trails prod (e.g. `propagation_scores` moved to binary files).
ROW_COUNT_SOURCES = [
    ("contacts", "qso_timestamp"),
    ("surface_observations", "observed_at"),
    ("soundings", "observed_at"),
    ("hrrr_profiles", "valid_time"),
    ("hrrr_native_profiles", "valid_time"),
    ("narr_profiles", "valid_time"),
    ("iemre_observations", None),
    ("nexrad_observations", "observed_at"),
    ("rtma_observations", "valid_time"),
    ("terrain_profiles", None),
    ("propagation_scores", "valid_time"),
]

CONTACTS_BY_BAND_SQL = """
SELECT band::int AS band_mhz, count(*) AS contacts,
       ROUND(AVG(distance_km::numeric)) AS avg_km,
       MAX(distance_km::numeric) AS max_km
FROM contacts
WHERE pos1 IS NOT NULL AND distance_km IS NOT NULL
  AND distance_km < 3000 AND flagged_invalid = false AND band >= 50
GROUP BY 1 ORDER BY 1;
"""

SOUNDING_MONTHLY_SQL = """
SELECT EXTRACT(MONTH FROM observed_at)::int AS month,
       count(*) AS soundings,
       ROUND(100.0 * count(*) FILTER (WHERE ducting_detected) / count(*), 1) AS ducting_pct,
       ROUND(AVG(min_refractivity_gradient)::numeric, 1) AS avg_min_grad,
       ROUND(AVG(precipitable_water_mm)::numeric, 1) AS avg_pwat
FROM soundings
WHERE min_refractivity_gradient IS NOT NULL
GROUP BY 1 ORDER BY 1;
"""

NATIVE_DUCT_SQL = """
SELECT
  CASE WHEN best_duct_band_ghz IS NULL THEN 'none'
       WHEN best_duct_band_ghz < 5 THEN '<5 GHz'
       WHEN best_duct_band_ghz < 15 THEN '5-15 GHz'
       WHEN best_duct_band_ghz < 30 THEN '15-30 GHz'
       WHEN best_duct_band_ghz < 75 THEN '30-75 GHz'
       ELSE '75+ GHz' END AS band_bin,
  count(*) AS profiles,
  ROUND(AVG(inversion_top_m)::numeric) AS avg_inv_top_m,
  ROUND(AVG(theta_e_jump_k)::numeric, 1) AS avg_theta_e_jump,
  ROUND(AVG(bulk_richardson)::numeric, 2) AS avg_richardson
FROM hrrr_native_profiles GROUP BY 1 ORDER BY 1;
"""

# Per-band contact ↔ HRRR Pearson correlations (2014-10 onward).
PER_BAND_JOIN_SQL = """
WITH joined AS (
  SELECT DISTINCT ON (c.id)
    c.id, c.band::int AS band, c.distance_km::float AS dist,
    h.surface_temp_c AS tc, h.surface_dewpoint_c AS dpc,
    h.surface_pressure_mb AS pr, h.pwat_mm AS pwat,
    h.min_refractivity_gradient AS grad,
    h.surface_refractivity AS sref, h.hpbl_m AS hpbl
  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 < 3000 AND c.flagged_invalid = false
  ORDER BY c.id, ABS(EXTRACT(EPOCH FROM h.valid_time - c.qso_timestamp))
)
SELECT * FROM joined;
"""

# Per-band contact ↔ NARR Pearson correlations (pre-2014-10).
# NARR is 32 km Lambert conformal / 3-hourly, so both tolerances are wider
# than the HRRR join above. Matches `NarrClient.in_coverage?/1` — valid up
# through but excluding 2014-10-02 (HRRR takes over after that).
PER_BAND_NARR_JOIN_SQL = """
WITH joined AS (
  SELECT DISTINCT ON (c.id)
    c.id, c.band::int AS band, c.distance_km::float AS dist,
    n.surface_temp_c AS tc, n.surface_dewpoint_c AS dpc,
    n.surface_pressure_mb AS pr, n.pwat_mm AS pwat,
    n.min_refractivity_gradient AS grad,
    n.surface_refractivity AS sref, n.hpbl_m AS hpbl
  FROM contacts c
  JOIN narr_profiles n
    ON n.lat BETWEEN (c.pos1->>'lat')::float - 0.25
                 AND (c.pos1->>'lat')::float + 0.25
   AND n.lon BETWEEN (c.pos1->>'lon')::float - 0.25
                 AND (c.pos1->>'lon')::float + 0.25
   AND n.valid_time BETWEEN c.qso_timestamp - INTERVAL '2 hours'
                        AND c.qso_timestamp + INTERVAL '2 hours'
  WHERE c.pos1 IS NOT NULL AND c.distance_km < 3000
    AND c.flagged_invalid = false
    AND c.qso_timestamp < TIMESTAMP '2014-10-02'
  ORDER BY c.id, ABS(EXTRACT(EPOCH FROM n.valid_time - c.qso_timestamp))
)
SELECT * FROM joined;
"""

# NEXRAD ↔ contact join. Very tight spatial tolerance because nexrad_observations
# is a per-contact enrichment table, not a grid sample.
NEXRAD_JOIN_SQL = """
WITH joined AS (
  SELECT DISTINCT ON (c.id)
    c.id, c.band::int AS band, c.distance_km::float AS dist,
    n.max_reflectivity_dbz AS max_dbz,
    n.mean_reflectivity_dbz AS mean_dbz
  FROM contacts c
  JOIN nexrad_observations n
    ON n.lat BETWEEN (c.pos1->>'lat')::float - 0.1
                 AND (c.pos1->>'lat')::float + 0.1
   AND n.lon BETWEEN (c.pos1->>'lon')::float - 0.1
                 AND (c.pos1->>'lon')::float + 0.1
   AND n.observed_at BETWEEN c.qso_timestamp - INTERVAL '15 min'
                        AND c.qso_timestamp + INTERVAL '15 min'
  WHERE c.pos1 IS NOT NULL AND c.distance_km < 3000 AND c.flagged_invalid = false
  ORDER BY c.id, ABS(EXTRACT(EPOCH FROM n.observed_at - c.qso_timestamp))
)
SELECT * FROM joined;
"""

# Native HRRR profile ↔ contact join — gives us best_duct_band_ghz per contact.
NATIVE_DUCT_JOIN_SQL = """
WITH joined AS (
  SELECT DISTINCT ON (c.id)
    c.id, c.band::int AS band, c.distance_km::float AS dist,
    n.best_duct_band_ghz AS duct_ghz,
    n.inversion_top_m AS inv_top_m,
    n.bulk_richardson AS richardson
  FROM contacts c
  JOIN hrrr_native_profiles n
    ON n.lat BETWEEN (c.pos1->>'lat')::float - 0.07
                 AND (c.pos1->>'lat')::float + 0.07
   AND n.lon BETWEEN (c.pos1->>'lon')::float - 0.07
                 AND (c.pos1->>'lon')::float + 0.07
   AND n.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 < 3000 AND c.flagged_invalid = false
  ORDER BY c.id, ABS(EXTRACT(EPOCH FROM n.valid_time - c.qso_timestamp))
)
SELECT * FROM joined;
"""

# Commercial-link degradation ↔ contact join. For each contact in the DFW zone
# (within ~150 km of the Princeton TX cluster at 33.2N, -96.5W), find the
# nearest commercial sample within ±15 min and compute the dB delta from the
# 7-day trailing baseline. The boost factor in Scorer.commercial_link_boost/2
# expects the same quantity, so this query answers "does our commercial-link
# signal actually correlate with longer QSO distances when we have both?"
DFW_CLUSTER_LAT = 33.2
DFW_CLUSTER_LON = -96.5
DFW_RADIUS_DEG = 1.5  # roughly 150 km at 33° latitude

COMMERCIAL_JOIN_SQL = f"""
WITH contact_zone AS (
  SELECT c.id, c.qso_timestamp, c.band::int AS band, c.distance_km::float AS dist
  FROM contacts c
  WHERE c.pos1 IS NOT NULL AND c.distance_km < 3000 AND c.flagged_invalid = false
    AND (c.pos1->>'lat')::float BETWEEN {DFW_CLUSTER_LAT - DFW_RADIUS_DEG}
                                    AND {DFW_CLUSTER_LAT + DFW_RADIUS_DEG}
    AND (c.pos1->>'lon')::float BETWEEN {DFW_CLUSTER_LON - DFW_RADIUS_DEG}
                                    AND {DFW_CLUSTER_LON + DFW_RADIUS_DEG}
),
current_sample AS (
  SELECT DISTINCT ON (cz.id)
    cz.id, cz.qso_timestamp, cz.band, cz.dist,
    cs.link_id, cs.rx_power_0 AS current_rx, cs.sampled_at
  FROM contact_zone cz
  JOIN commercial_samples cs
    ON cs.sampled_at BETWEEN cz.qso_timestamp - INTERVAL '15 min'
                         AND cz.qso_timestamp + INTERVAL '15 min'
   AND cs.link_state = 1 AND cs.rx_power_0 IS NOT NULL
  ORDER BY cz.id, ABS(EXTRACT(EPOCH FROM cs.sampled_at - cz.qso_timestamp))
),
baseline AS (
  SELECT cs.id, cs.qso_timestamp, cs.band, cs.dist, cs.link_id, cs.current_rx,
         (
           SELECT AVG(bs.rx_power_0)
           FROM commercial_samples bs
           WHERE bs.link_id = cs.link_id
             AND bs.link_state = 1
             AND bs.rx_power_0 IS NOT NULL
             AND bs.sampled_at BETWEEN cs.sampled_at - INTERVAL '7 days'
                                  AND cs.sampled_at - INTERVAL '15 min'
         ) AS baseline_rx
  FROM current_sample cs
)
SELECT id, band, dist,
       (baseline_rx - current_rx)::float AS degradation_db,
       current_rx::float AS current_dbm,
       baseline_rx::float AS baseline_dbm
FROM baseline
WHERE baseline_rx IS NOT NULL;
"""

# Empty-table / data-quality smoke checks.
DATA_GAP_SQL = """
SELECT
  (SELECT count(*) FROM narr_profiles) AS narr_rows,
  (SELECT count(*) FROM narr_profiles WHERE valid_time < '2014-10-02') AS narr_pre2014_rows,
  (SELECT count(*) FROM hrrr_climatology) AS hrrr_climatology_rows,
  (SELECT count(*) FROM rtma_observations) AS rtma_rows,
  (SELECT count(*) FROM metar_5min_observations) AS metar5_rows,
  (SELECT count(*) FROM contacts WHERE hrrr_status = 'complete') AS hrrr_complete_contacts,
  (SELECT count(*) FROM contacts WHERE hrrr_status != 'complete') AS hrrr_pending_contacts,
  (SELECT count(*) FROM contacts
      WHERE qso_timestamp < '2014-10-02'
        AND pos1 IS NOT NULL AND distance_km < 3000
        AND flagged_invalid = false) AS pre2014_contacts;
"""


# ─── Reporting helpers ──────────────────────────────────────────────────────


def fetch_df(conn, sql: str) -> pd.DataFrame:
    with conn.cursor(row_factory=dict_row) as cur:
        cur.execute(sql)
        rows = cur.fetchall()
    return pd.DataFrame(rows)


def md_table(df: pd.DataFrame) -> str:
    if df.empty:
        return "_(no rows)_\n"
    return df.to_markdown(index=False, floatfmt=".3f") + "\n"


def fetch_row_counts(conn) -> pd.DataFrame:
    """Per-table counts + date ranges, tolerant of missing tables."""
    rows = []
    for tbl, time_col in ROW_COUNT_SOURCES:
        with conn.cursor(row_factory=dict_row) as cur:
            cur.execute("SELECT to_regclass(%s) AS oid", (f"public.{tbl}",))
            present = cur.fetchone()["oid"] is not None
        if not present:
            rows.append({"tbl": tbl, "n": None, "lo": None, "hi": None})
            continue
        if time_col:
            sql = (
                f"SELECT '{tbl}' tbl, count(*)::bigint n, "
                f"min({time_col})::date lo, max({time_col})::date hi FROM {tbl}"
            )
        else:
            sql = (
                f"SELECT '{tbl}' tbl, count(*)::bigint n, "
                f"NULL::date lo, NULL::date hi FROM {tbl}"
            )
        with conn.cursor(row_factory=dict_row) as cur:
            cur.execute(sql)
            rows.append(cur.fetchone())
    return pd.DataFrame(rows).sort_values("tbl")


def band_label(band_mhz: int) -> str:
    """Human-readable band label, matching `BandConfig` conventions."""
    if band_mhz < 1_000:
        return f"{band_mhz} MHz"
    ghz = band_mhz / 1_000
    if ghz == int(ghz):
        return f"{int(ghz)} GHz"
    return f"{ghz:g} GHz"


def bands_with_samples(df: pd.DataFrame, min_samples: int = 50) -> list[int]:
    """Every band in `df` with at least `min_samples` rows, sorted ascending."""
    if df.empty or "band" not in df.columns:
        return []
    counts = df.groupby("band").size()
    return sorted(int(b) for b, n in counts.items() if n >= min_samples)


def correlations_per_band(df: pd.DataFrame, min_samples: int = 30) -> pd.DataFrame:
    fields = ["pr", "dpc", "pwat", "sref", "grad", "tc", "hpbl"]
    rows = []
    for band, group in df.groupby("band"):
        if len(group) < min_samples:
            continue
        row = {"band_mhz": int(band), "n": len(group)}
        for f in fields:
            try:
                row[f"rho_{f}"] = round(group["dist"].corr(group[f]), 3)
            except Exception:
                row[f"rho_{f}"] = None
        rows.append(row)
    return pd.DataFrame(rows).sort_values("band_mhz")


def hpbl_bins(df: pd.DataFrame, band_mhz: int = 10000) -> pd.DataFrame:
    sub = df[(df["band"] == band_mhz) & df["hpbl"].notna()].copy()
    if sub.empty:
        return pd.DataFrame()
    bins = [0, 200, 500, 1000, 1500, 2000, 1e9]
    labels = ["<200", "200-500", "500-1000", "1000-1500", "1500-2000", ">=2000"]
    sub["bin"] = pd.cut(sub["hpbl"], bins=bins, labels=labels, right=False)
    out = (
        sub.groupby("bin", observed=True)["dist"]
        .agg(["count", "mean", "median"])
        .round(1)
        .rename(columns={"count": "n", "mean": "avg_km", "median": "p50_km"})
        .reset_index()
    )
    return out


def pressure_bins(df: pd.DataFrame, band_mhz: int = 10000) -> pd.DataFrame:
    sub = df[(df["band"] == band_mhz) & df["pr"].notna()].copy()
    if sub.empty:
        return pd.DataFrame()
    bins = [0, 990, 1000, 1010, 1020, 9999]
    labels = ["<990", "990-1000", "1000-1010", "1010-1020", ">=1020"]
    sub["bin"] = pd.cut(sub["pr"], bins=bins, labels=labels, right=False)
    out = (
        sub.groupby("bin", observed=True)["dist"]
        .agg(["count", "mean", "median"])
        .round(1)
        .rename(columns={"count": "n", "mean": "avg_km", "median": "p50_km"})
        .reset_index()
    )
    return out


def nexrad_bins(df: pd.DataFrame, band_mhz: int = 10000) -> pd.DataFrame:
    """Distance distribution by NEXRAD max reflectivity bin, per band.

    Confirms the rain-attenuation direction of the scoring: as max_dbz rises,
    distance should fall at rain-sensitive bands (24+ GHz). Below-rain dBZ
    (<5) is the no-precip baseline.
    """
    sub = df[(df["band"] == band_mhz) & df["max_dbz"].notna()].copy()
    if sub.empty:
        return pd.DataFrame()
    bins = [-1e9, 5, 20, 30, 45, 1e9]
    labels = ["none (<5)", "drizzle (5-20)", "light (20-30)", "moderate (30-45)", "heavy (>=45)"]
    sub["bin"] = pd.cut(sub["max_dbz"], bins=bins, labels=labels, right=False)
    out = (
        sub.groupby("bin", observed=True)["dist"]
        .agg(["count", "mean", "median"])
        .round(1)
        .rename(columns={"count": "n", "mean": "avg_km", "median": "p50_km"})
        .reset_index()
    )
    return out


def native_duct_bins(df: pd.DataFrame, band_mhz: int = 10000) -> pd.DataFrame:
    """Distance distribution by native-profile best duct band, per band.

    Checks whether hrrr_native_profiles.best_duct_band_ghz actually predicts
    longer QSO distances at the target band — the assumption behind
    Scorer.score_refractivity/4's 1.15× native-duct boost.
    """
    sub = df[df["band"] == band_mhz].copy()
    if sub.empty:
        return pd.DataFrame()

    def label(v: object) -> str:
        if pd.isna(v):
            return "none"
        v = float(v)
        if v < 5:
            return "<5 GHz"
        if v < 10:
            return "5-10 GHz"
        if v < 24:
            return "10-24 GHz"
        if v < 47:
            return "24-47 GHz"
        return ">=47 GHz"

    sub["bin"] = sub["duct_ghz"].map(label)
    order = ["none", "<5 GHz", "5-10 GHz", "10-24 GHz", "24-47 GHz", ">=47 GHz"]
    out = (
        sub.groupby("bin")["dist"]
        .agg(["count", "mean", "median"])
        .round(1)
        .rename(columns={"count": "n", "mean": "avg_km", "median": "p50_km"})
        .reindex(order)
        .dropna(subset=["n"])
        .reset_index()
    )
    return out


def commercial_degradation_summary(df: pd.DataFrame) -> pd.DataFrame:
    """Does commercial-link rx_power degradation correlate with longer contacts?

    Returns a per-band bin table. Empty today because commercial_samples only
    covers ~2 weeks and contacts only cluster in Aug/Sep contests — but the
    query is here so the moment contest-season samples land the signal shows up.
    """
    if df.empty:
        return pd.DataFrame()
    bins = [-1e9, 0, 3, 8, 1e9]
    labels = ["negative/0 (no fade)", "3 dB noise floor", "3-8 dB mild", ">=8 dB strong"]
    df = df.copy()
    df["bin"] = pd.cut(df["degradation_db"], bins=bins, labels=labels, right=False)
    out = (
        df.groupby("bin", observed=True)["dist"]
        .agg(["count", "mean", "median"])
        .round(1)
        .rename(columns={"count": "n", "mean": "avg_km", "median": "p50_km"})
        .reset_index()
    )
    return out


def commercial_correlation(df: pd.DataFrame) -> float | None:
    """Pearson correlation of degradation_db with contact distance. None if n<5."""
    if df.empty or len(df) < 5:
        return None
    return round(df["dist"].corr(df["degradation_db"]), 3)


def redact(dsn: str) -> str:
    if "@" not in dsn:
        return dsn
    head, tail = dsn.split("@", 1)
    if "://" in head and ":" in head.split("://", 1)[1]:
        scheme, rest = head.split("://", 1)
        user = rest.split(":", 1)[0]
        return f"{scheme}://{user}:***@{tail}"
    return dsn


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


def main() -> int:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument(
        "--dsn",
        default=os.environ.get("PROP_PROD_DB_URL"),
        help="Postgres connection string (defaults to $PROP_PROD_DB_URL)",
    )
    parser.add_argument(
        "--out",
        default=None,
        help="Output Markdown path (defaults to docs/algo-reports/YYYY-MM-DD-recalibration.md)",
    )
    parser.add_argument(
        "--statement-timeout",
        default="20min",
        help="Postgres statement_timeout (default: 20min)",
    )
    args = parser.parse_args()

    if not args.dsn:
        print(
            "error: PROP_PROD_DB_URL not set and --dsn not given; populate "
            "via .envrc",
            file=sys.stderr,
        )
        return 2

    today = dt.date.today().isoformat()
    out_path = Path(
        args.out
        or Path(__file__).resolve().parent.parent
        / "docs"
        / "algo-reports"
        / f"{today}-recalibration.md"
    )
    out_path.parent.mkdir(parents=True, exist_ok=True)

    print(f"connecting to {redact(args.dsn)}", file=sys.stderr)

    sections: list[str] = []
    sections.append(
        REPORT_HEADER.format(
            today=today,
            dsn_redacted=redact(args.dsn),
            stmt_timeout=args.statement_timeout,
        )
    )

    with psycopg.connect(args.dsn, autocommit=True) as conn:
        with conn.cursor() as cur:
            cur.execute(f"SET statement_timeout = '{args.statement_timeout}';")

        print("• row counts", file=sys.stderr)
        sections.append("\n## Row counts and date ranges\n")
        sections.append(md_table(fetch_row_counts(conn)))

        print("• data gaps", file=sys.stderr)
        gaps = fetch_df(conn, DATA_GAP_SQL).iloc[0].to_dict()
        sections.append("\n## Data gaps (these should NOT be zero in a healthy prod)\n")
        sections.append(
            "\n".join(f"- **{k}**: {v}" for k, v in gaps.items()) + "\n"
        )

        print("• contacts by band", file=sys.stderr)
        sections.append("\n## Contacts by band (≥50 MHz)\n")
        sections.append(md_table(fetch_df(conn, CONTACTS_BY_BAND_SQL)))

        print("• sounding monthly", file=sys.stderr)
        sections.append("\n## Monthly sounding ducting\n")
        sections.append(md_table(fetch_df(conn, SOUNDING_MONTHLY_SQL)))

        print("• native HRRR ducts", file=sys.stderr)
        sections.append(
            "\n## HRRR native-profile duct distribution (best supportable band)\n"
        )
        sections.append(md_table(fetch_df(conn, NATIVE_DUCT_SQL)))

        print("• per-band correlations — HRRR (this is the slow one)", file=sys.stderr)
        joined = fetch_df(conn, PER_BAND_JOIN_SQL)

        print("• per-band correlations — NARR (pre-2014)", file=sys.stderr)
        joined_narr = fetch_df(conn, PER_BAND_NARR_JOIN_SQL)

        print("• NEXRAD ↔ contacts", file=sys.stderr)
        nexrad = fetch_df(conn, NEXRAD_JOIN_SQL)

        print("• native-duct ↔ contacts", file=sys.stderr)
        native = fetch_df(conn, NATIVE_DUCT_JOIN_SQL)

        print("• commercial-link ↔ contacts (DFW zone only)", file=sys.stderr)
        commercial = fetch_df(conn, COMMERCIAL_JOIN_SQL)

    if joined.empty:
        sections.append("\n_No contacts ↔ HRRR matches were found._\n")
    else:
        sections.append(
            f"\n## Contact ↔ HRRR per-band Pearson correlations "
            f"(matched n={len(joined)})\n"
        )
        sections.append(md_table(correlations_per_band(joined)))

        sections.append(
            "\n## Per-band HPBL bin distance distribution\n\n"
            "One table per band with ≥50 matched contacts. A band-specific "
            "effect here means HPBL belongs in the scorer for that band.\n"
        )
        for band in bands_with_samples(joined):
            table = hpbl_bins(joined, band)
            if table.empty:
                continue
            sections.append(f"\n**{band_label(band)}:**\n")
            sections.append(md_table(table))

        sections.append(
            "\n## Per-band pressure bin distance distribution\n\n"
            "Surface pressure tends to proxy synoptic-scale stagnation, so we "
            "want to see longer distances under the 1015-1025 mb ridge at every "
            "band that ducts.\n"
        )
        for band in bands_with_samples(joined):
            table = pressure_bins(joined, band)
            if table.empty:
                continue
            sections.append(f"\n**{band_label(band)}:**\n")
            sections.append(md_table(table))

    if joined_narr.empty:
        sections.append(
            "\n## Contact ↔ NARR per-band Pearson correlations\n\n"
            "_No pre-2014 contacts ↔ NARR matches were found. Either the NARR "
            "backfill has not run for this corpus yet or there are no "
            "pre-2014 contacts with `pos1` set._\n"
        )
    else:
        sections.append(
            f"\n## Contact ↔ NARR per-band Pearson correlations "
            f"(pre-2014 only, matched n={len(joined_narr)})\n\n"
            "Historical sanity check: if the HRRR-era signs and magnitudes "
            "survive on 30+ years of NARR reanalysis, the scoring factors are "
            "physical rather than HRRR-artefact. Diverging signs = flag it.\n"
        )
        sections.append(md_table(correlations_per_band(joined_narr)))

    sections.append(
        f"\n## NEXRAD composite reflectivity vs distance "
        f"(matched n={len(nexrad)})\n"
    )
    if nexrad.empty:
        sections.append("_No contacts ↔ NEXRAD matches were found._\n")
    else:
        sections.append(
            "Rain-attenuation sanity check — one table per band with ≥50 "
            "matched contacts. At rain-sensitive bands (24+ GHz) higher "
            "`max_dbz` should map to shorter contacts; at 10 GHz and below "
            "the effect should be barely detectable.\n"
        )
        for band in bands_with_samples(nexrad):
            table = nexrad_bins(nexrad, band)
            if table.empty:
                continue
            sections.append(f"\n**{band_label(band)}:**\n")
            sections.append(md_table(table))

    sections.append(
        f"\n## hrrr_native_profiles.best_duct_band_ghz vs distance "
        f"(matched n={len(native)})\n"
    )
    if native.empty:
        sections.append("_No contacts ↔ native-profile matches were found._\n")
    else:
        sections.append(
            "Validates the 1.15× boost in `Scorer.score_refractivity/4`: "
            "contacts where the native duct supports the target band should "
            "run longer than those where it does not. One table per band with "
            "≥50 matched contacts.\n"
        )
        for band in bands_with_samples(native):
            table = native_duct_bins(native, band)
            if table.empty:
                continue
            sections.append(f"\n**{band_label(band)}:**\n")
            sections.append(md_table(table))

    sections.append(
        f"\n## Commercial-link rx_power degradation vs contemporaneous DFW contacts "
        f"(matched n={len(commercial)})\n"
    )
    if commercial.empty:
        sections.append(
            "_No DFW-zone contacts fall inside the commercial_samples date "
            "window yet. Expected to stay empty until the Aug/Sep contest "
            "season produces contacts that overlap with live SNMP polling._\n"
        )
    else:
        rho = commercial_correlation(commercial)
        if rho is not None:
            sections.append(f"Pearson(degradation_db, distance_km) = **{rho}**\n\n")
        sections.append(md_table(commercial_degradation_summary(commercial)))

    out_path.write_text("\n".join(sections))
    print(f"wrote {out_path}", file=sys.stderr)
    return 0


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