From 5af4142e28a4a12d683cbf25a63d1ef99d340ac1 Mon Sep 17 00:00:00 2001 From: Graham McIntire Date: Mon, 13 Apr 2026 12:13:54 -0500 Subject: [PATCH] Extend recalibrate_algo.py with NEXRAD, native-duct, commercial sections MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three new analysis passes in scripts/recalibrate_algo.py that match the scoring factors landed in commit 86364f4: * NEXRAD composite reflectivity vs distance per band — checks that the Marshall-Palmer rain rate in Scorer.dbz_to_rain_rate_mmhr/1 moves the needle in the expected direction at 24+ GHz. * hrrr_native_profiles.best_duct_band_ghz vs distance per band — validates the 1.15× boost in Scorer.score_refractivity/4. Current corpus shows essentially zero contacts with a native duct supporting ≥5 GHz at the target band, consistent with the 0.12% rate in Part 2c; the boost is physics-correct but almost never fires today. * Commercial-link rx_power degradation vs contemporaneous DFW contacts — pearson + binned distance. Empty today because commercial_samples only covers Mar 30–Apr 13 and contests don't overlap. Scaffolding is in place so the moment contest-season overlap exists, the signal shows. Report rebuilt at docs/algo-reports/2026-04-13-recalibration.md. --- docs/algo-reports/2026-04-13-recalibration.md | 85 +++++- scripts/recalibrate_algo.py | 252 ++++++++++++++++++ 2 files changed, 324 insertions(+), 13 deletions(-) diff --git a/docs/algo-reports/2026-04-13-recalibration.md b/docs/algo-reports/2026-04-13-recalibration.md index e52f9972..5cf44e0b 100644 --- a/docs/algo-reports/2026-04-13-recalibration.md +++ b/docs/algo-reports/2026-04-13-recalibration.md @@ -15,13 +15,13 @@ Statement timeout: 20min | contacts | 58560 | 1991-05-04 | 2026-03-07 | | era5_profiles | 0 | | | | hrrr_native_profiles | 11472 | 2019-03-09 | 2025-09-21 | -| hrrr_profiles | 76058705 | 2016-09-17 | 2026-04-13 | +| hrrr_profiles | 76151193 | 2016-09-17 | 2026-04-13 | | iemre_observations | 15745 | | | | nexrad_observations | 7286 | 2019-08-17 | 2024-09-22 | -| propagation_scores | 9618752 | 2026-04-13 | 2026-04-13 | +| propagation_scores | 10543632 | 2026-04-13 | 2026-04-13 | | rtma_observations | 0 | | | | soundings | 16864 | 1990-06-07 | 2026-03-07 | -| surface_observations | 392442 | 1991-05-03 | 2026-04-13 | +| surface_observations | 392444 | 1991-05-03 | 2026-04-13 | | terrain_profiles | 58554 | | | @@ -93,20 +93,20 @@ Statement timeout: 20min | band_mhz | n | rho_pr | rho_dpc | rho_pwat | rho_sref | rho_grad | rho_tc | rho_hpbl | |-----------:|--------:|---------:|----------:|-----------:|-----------:|-----------:|---------:|-----------:| | 1296.000 | 63.000 | -0.446 | 0.150 | -0.004 | 0.321 | -0.308 | -0.313 | -0.383 | -| 10000.000 | 680.000 | -0.104 | 0.046 | 0.003 | 0.181 | -0.127 | -0.108 | -0.200 | -| 24000.000 | 147.000 | -0.290 | -0.179 | -0.220 | 0.096 | -0.034 | -0.301 | -0.226 | -| 47000.000 | 75.000 | -0.111 | 0.015 | 0.006 | 0.218 | -0.167 | 0.086 | -0.154 | +| 10000.000 | 680.000 | -0.109 | 0.044 | 0.002 | 0.182 | -0.127 | -0.107 | -0.198 | +| 24000.000 | 147.000 | -0.290 | -0.179 | -0.219 | 0.096 | -0.036 | -0.300 | -0.223 | +| 47000.000 | 75.000 | -0.111 | 0.016 | 0.007 | 0.217 | -0.175 | 0.089 | -0.154 | ## 10 GHz: HPBL bin distance distribution | bin | n | avg_km | p50_km | |:----------|----:|---------:|---------:| -| <200 | 55 | 222.600 | 174.000 | -| 200-500 | 174 | 183.200 | 154.500 | -| 500-1000 | 181 | 163.700 | 126.000 | -| 1000-1500 | 103 | 160.900 | 136.000 | -| 1500-2000 | 93 | 140.500 | 128.000 | +| <200 | 54 | 217.600 | 174.000 | +| 200-500 | 177 | 184.800 | 160.000 | +| 500-1000 | 177 | 163.800 | 126.000 | +| 1000-1500 | 104 | 158.800 | 136.000 | +| 1500-2000 | 94 | 142.800 | 136.000 | | >=2000 | 74 | 100.300 | 90.000 | @@ -114,8 +114,67 @@ Statement timeout: 20min | bin | n | avg_km | p50_km | |:----------|----:|---------:|---------:| -| <990 | 101 | 211.900 | 134.000 | -| 990-1000 | 538 | 155.400 | 128.000 | +| <990 | 103 | 214.200 | 136.000 | +| 990-1000 | 536 | 154.800 | 128.000 | | 1000-1010 | 28 | 164.600 | 115.000 | | 1010-1020 | 2 | 23.000 | 23.000 | | >=1020 | 11 | 103.100 | 23.000 | + + +## NEXRAD composite reflectivity vs distance (matched n=102) + +Sanity check for the rain-attenuation direction: at rain-sensitive bands higher `max_dbz` should map to shorter contacts. At 10 GHz the effect should be barely detectable. + + +**10 GHz:** + +| bin | n | avg_km | p50_km | +|:-----------------|----:|---------:|---------:| +| none (<5) | 59 | 155.700 | 143.000 | +| drizzle (5-20) | 7 | 56.400 | 49.000 | +| light (20-30) | 13 | 95.200 | 40.000 | +| moderate (30-45) | 4 | 97.200 | 92.000 | +| heavy (>=45) | 1 | 213.000 | 213.000 | + + +**24 GHz:** + +| bin | n | avg_km | p50_km | +|:-----------------|----:|---------:|---------:| +| none (<5) | 5 | 29.800 | 23.000 | +| drizzle (5-20) | 1 | 49.000 | 49.000 | +| light (20-30) | 5 | 34.000 | 40.000 | +| moderate (30-45) | 1 | 60.000 | 60.000 | + + +## hrrr_native_profiles.best_duct_band_ghz vs distance (matched n=619) + +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. + + +**10 GHz:** + +| bin | n | avg_km | p50_km | +|:-------|--------:|---------:|---------:| +| none | 447.000 | 140.800 | 128.000 | +| <5 GHz | 38.000 | 150.000 | 114.500 | + + +**24 GHz:** + +| bin | n | avg_km | p50_km | +|:-------|-------:|---------:|---------:| +| none | 97.000 | 47.600 | 32.000 | +| <5 GHz | 3.000 | 56.300 | 60.000 | + + +**47 GHz:** + +| bin | n | avg_km | p50_km | +|:------|-------:|---------:|---------:| +| none | 16.000 | 54.700 | 42.000 | + + +## Commercial-link rx_power degradation vs contemporaneous DFW contacts (matched n=0) + +_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._ diff --git a/scripts/recalibrate_algo.py b/scripts/recalibrate_algo.py index 2603d9df..c62dc995 100755 --- a/scripts/recalibrate_algo.py +++ b/scripts/recalibrate_algo.py @@ -12,6 +12,9 @@ and produces a Markdown report covering: * native HRRR profile duct stats * per-band Pearson correlations of contact distance vs HRRR fields * binned distance distributions for HPBL and pressure at 10 GHz + * 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. @@ -148,6 +151,102 @@ WITH joined AS ( 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 @@ -226,6 +325,98 @@ def pressure_bins(df: pd.DataFrame, band_mhz: int = 10000) -> pd.DataFrame: 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 @@ -320,6 +511,15 @@ def main() -> int: print("• per-band correlations (this is the slow one)", file=sys.stderr) joined = fetch_df(conn, PER_BAND_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: @@ -335,6 +535,58 @@ def main() -> int: sections.append("\n## 10 GHz: pressure bin distance distribution\n") sections.append(md_table(pressure_bins(joined))) + 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( + "Sanity check for the rain-attenuation direction: at rain-sensitive " + "bands higher `max_dbz` should map to shorter contacts. At 10 GHz " + "the effect should be barely detectable.\n\n" + ) + sections.append("**10 GHz:**\n") + sections.append(md_table(nexrad_bins(nexrad, 10_000))) + sections.append("\n**24 GHz:**\n") + sections.append(md_table(nexrad_bins(nexrad, 24_000))) + + 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.\n\n" + ) + sections.append("**10 GHz:**\n") + sections.append(md_table(native_duct_bins(native, 10_000))) + sections.append("\n**24 GHz:**\n") + sections.append(md_table(native_duct_bins(native, 24_000))) + sections.append("\n**47 GHz:**\n") + sections.append(md_table(native_duct_bins(native, 47_000))) + + 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