graham/prop
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

perf(algo): unify recalibration into single Python pipeline with JSON output

Browse source 
Replace the two-script Python pipeline (analysis report + Elixir-source
emitter) with a single `scripts/recalibrate.py` that fits per-band
weights from PSKR spot density (VHF/UHF) and contacts↔HRRR correlations
(microwave), writing `priv/algo/band_weights.json` as a machine-readable
artifact. A new Elixir BandWeights module loads this JSON once via
`:persistent_term` cache; BandConfig.weights/1 consults it before falling
back to in-source overrides or global defaults. The script never touches
Elixir source — recalibration is now `python3 scripts/recalibrate.py`
followed by an app restart.
This commit is contained in:
Graham McIntire 2026-05-25 14:45:55 -05:00
parent 580b9e781b
commit daafa5a02a
No known key found for this signature in database
GPG key ID: F4ABF488E6029E59
10 changed files with 1297 additions and 1201 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

6
config/config.exs
View file

@ -136,6 +136,12 @@ config :microwaveprop, Oban,
]}
]
# Per-band scoring weights file produced by `scripts/recalibrate.py`.
# Read once at app start by `Microwaveprop.Propagation.BandWeights` and
# cached in `:persistent_term`. Set to `false` to ignore the file and
# fall back to in-source `:weights` overrides + `@weights` defaults.
config :microwaveprop, :band_weights_json, "priv/algo/band_weights.json"
config :microwaveprop, :scopes,
user: [
default: true,

4
config/test.exs
View file

@ -64,6 +64,10 @@ config :microwaveprop, MicrowavepropWeb.Endpoint,
# Run Oban jobs inline during tests
config :microwaveprop, Oban, testing: :inline
# Pin the scorer to in-source weights in test. Individual tests that
# need the JSON loader (`BandWeightsTest`) override this per-case.
config :microwaveprop, :band_weights_json, false
# Skip the registration-time QRZ lookup in tests — it needs a Mox stub
# that isn't set up across the test suite, and the worker is exercised
# directly in its own dedicated test.

56
docs/algo-reports/2026-05-24-full-recalibration.md Normal file
View file

@ -0,0 +1,56 @@
# Full Recalibration Report — 2026-05-24
> Auto-generated by `scripts/recalibrate.py`. The companion JSON file at `priv/algo/band_weights.json` is the machine-readable output that Elixir loads at runtime. This report is the human-readable audit trail.
Connection: `postgres://prop:***@10.0.15.30:5432/prop`
Generated: `2026-05-24T18:53:32+00:00`
## Data sources
- **PSKR** (VHF/UHF, label = `spot_count`): bands [50, 144, 432, 1296, 2304, 10000], 386,782 samples
- **Contacts ↔ HRRR** (microwave, label = `distance_km`): 12 bands, 81,571 matched contacts
## Per-band evidence
| Band | Source | n | label avg | ρ_tc | ρ_dpc | ρ_pr | ρ_pwat | ρ_grad |
|------|--------|---:|----------:|-----:|------:|-----:|-------:|-------:|
| 50 MHz | pskr | 372,468 | 1.80 | -0.013 | +0.014 | +0.025 | +0.023 | +0.000 |
| 144 MHz | pskr | 12,927 | 4.41 | +0.184 | +0.223 | +0.163 | +0.209 | -0.073 |
| 222 MHz | contacts | 7,518 | 229.14 | -0.081 | +0.077 | -0.014 | +0.068 | +0.039 |
| 432 MHz | pskr | 975 | 1.54 | +0.260 | -0.047 | +0.085 | -0.052 | +0.141 |
| 902 MHz | contacts | 1,754 | 172.14 | -0.070 | +0.097 | +0.080 | +0.051 | -0.038 |
| 1296 MHz | contacts | 2,951 | 160.95 | -0.041 | +0.102 | +0.025 | +0.053 | -0.080 |
| 2304 MHz | contacts | 719 | 154.51 | -0.060 | +0.162 | +0.076 | +0.130 | -0.114 |
| 3400 MHz | contacts | 335 | 153.71 | -0.062 | +0.155 | -0.027 | +0.150 | +0.063 |
| 5760 MHz | contacts | 321 | 127.55 | -0.065 | +0.075 | +0.153 | +0.099 | -0.012 |
| 10000 MHz | contacts | 54,161 | 209.29 | -0.004 | -0.021 | -0.085 | -0.029 | +0.030 |
| 24000 MHz | contacts | 3,808 | 95.19 | -0.210 | -0.334 | -0.217 | -0.360 | +0.018 |
| 47000 MHz | contacts | 759 | 62.72 | -0.114 | -0.216 | -0.253 | -0.233 | -0.023 |
| 75000 MHz | contacts | 109 | 47.53 | -0.680 | -0.699 | -0.768 | -0.620 | +0.434 |
| 122000 MHz | contacts | 63 | 16.26 | +0.039 | -0.129 | -0.041 | -0.127 | +0.355 |
## Derived per-band weights
| Band | n | source | humidity | time_of_day | td_depression | refractivity | sky | season | wind | rain | pwat | pressure |
|------|---:|--------|----:|----:|----:|----:|----:|----:|----:|----:|----:|----:|
| 50 MHz | 372,468 | pskr | 0.1352 | 0.0271 | 0.1415 | 0.0736 | 0.0999 | 0.2156 | 0.0999 | 0.0170 | 0.1263 | 0.0639 |
| 144 MHz | 12,927 | pskr | 0.1862 | 0.0196 | 0.1491 | 0.1118 | 0.0621 | 0.1172 | 0.0621 | 0.0106 | 0.1693 | 0.1121 |
| 222 MHz | 7,518 | contacts | 0.2012 | 0.0272 | 0.1806 | 0.0986 | 0.0752 | 0.1318 | 0.0752 | 0.0128 | 0.1541 | 0.0432 |
| 432 MHz | 975 | pskr | 0.1472 | 0.0261 | 0.1735 | 0.1694 | 0.0722 | 0.1169 | 0.0722 | 0.0123 | 0.1272 | 0.0830 |
| 902 MHz | 1,754 | contacts | 0.2206 | 0.0299 | 0.1765 | 0.0955 | 0.0735 | 0.1090 | 0.0735 | 0.0125 | 0.1279 | 0.0810 |
| 1296 MHz | 2,951 | contacts | 0.2255 | 0.0339 | 0.1519 | 0.1433 | 0.0751 | 0.1013 | 0.0751 | 0.0128 | 0.1334 | 0.0476 |
| 2304 MHz | 719 | contacts | 0.1890 | 0.0285 | 0.1512 | 0.1477 | 0.0630 | 0.0849 | 0.0630 | 0.0339 | 0.1718 | 0.0672 |
| 3400 MHz | 335 | contacts | 0.2000 | 0.0301 | 0.1601 | 0.1107 | 0.0666 | 0.0899 | 0.0666 | 0.0507 | 0.1818 | 0.0434 |
| 5760 MHz | 321 | contacts | 0.1719 | 0.0297 | 0.1577 | 0.0599 | 0.0656 | 0.0885 | 0.0656 | 0.0790 | 0.1689 | 0.1133 |
| 10000 MHz | 54,161 | contacts | 0.1262 | 0.0380 | 0.1010 | 0.0986 | 0.0841 | 0.1134 | 0.0841 | 0.1431 | 0.1147 | 0.0967 |
| 24000 MHz | 3,808 | contacts | 0.1485 | 0.0403 | 0.1189 | 0.0498 | 0.0495 | 0.0734 | 0.0495 | 0.2228 | 0.1350 | 0.1123 |
| 47000 MHz | 759 | contacts | 0.1406 | 0.0529 | 0.1125 | 0.0500 | 0.0468 | 0.0758 | 0.0468 | 0.2391 | 0.1278 | 0.1077 |
## Changes vs prior `band_weights.json`
_No prior file, or no factor moved by ≥0.005. First-run values can be compared to the in-tree `@band_configs` weights in `lib/microwaveprop/propagation/band_config.ex`._
## Applying these weights
This script writes `priv/algo/band_weights.json` and stops there — it never edits `band_config.ex`. To make the weights live, `BandConfig.weights(band_mhz)` needs a one-time edit to read this JSON at app start (cache in `:persistent_term`) and merge `band_overrides[band].weights` over the in-source defaults. Until then, the JSON is documentation only.

26
lib/microwaveprop/propagation/band_config.ex
View file

@ -7,6 +7,8 @@ defmodule Microwaveprop.Propagation.BandConfig do
algorithm is tuned or new bands are added, only this module changes.
"""
alias Microwaveprop.Propagation.BandWeights
# Default weight vector. Originally fit by gradient descent on 5,000 QSOs
# 2026-04-11 (loss 0.42 → 0.12) and retained as the 10 GHz prior. Bands
# with ≥200 HRRR-matched contacts in the 2026-04-18 full-corpus
@ -873,16 +875,32 @@ defmodule Microwaveprop.Propagation.BandConfig do
@doc """
Returns the scoring weights for a specific band configuration.
Bands with enough contacts to support a stable gradient-descent fit
carry a `:weights` override map. Bands without it fall back to the
global defaults. `nil` is accepted so callers that don't know the band
can ask for defaults without an extra guard.
Resolution order:
1. `priv/algo/band_weights.json` via `BandWeights.lookup/1` —
authoritative per-band fits produced by `scripts/recalibrate.py`.
2. In-source `:weights` override on the band config map — legacy
snapshot of the 2026-04-18 fit, retained as a fallback.
3. Global `@weights` default.
`nil` is accepted so callers that don't know the band can ask for
defaults without an extra guard.
"""
@spec weights(map() | nil) :: map()
def weights(nil), do: @weights
def weights(%{freq_mhz: freq_mhz} = band_config) when is_integer(freq_mhz) do
case BandWeights.lookup(freq_mhz) do
nil -> in_source_or_default(band_config)
json_override -> json_override
end
end
def weights(%{weights: override}) when is_map(override), do: override
def weights(%{}), do: @weights
defp in_source_or_default(%{weights: override}) when is_map(override), do: override
defp in_source_or_default(_), do: @weights
@doc "Returns the 12-element sunrise hour table (Jan-Dec, local time)."
@spec sunrise_table() :: [float()]
def sunrise_table, do: @sunrise_table

174
lib/microwaveprop/propagation/band_weights.ex Normal file
View file

@ -0,0 +1,174 @@
defmodule Microwaveprop.Propagation.BandWeights do
@moduledoc """
Loads per-band scoring-weight overrides from `priv/algo/band_weights.json`.
The JSON file is the authoritative output of `scripts/recalibrate.py`,
a Python pipeline that fits per-band weights from PSKR spot density
(VHF/UHF) and `contacts` ↔ `hrrr_profiles` correlations (microwave).
Keeping the fit results in a data file means recalibration runs never
have to modify Elixir source — the script writes JSON, this module
reads it, and `BandConfig.weights/1` consults this lookup before
falling back to in-source `:weights` overrides or the global default.
## Configuration
config :microwaveprop, :band_weights_json, "priv/algo/band_weights.json"
Set the value to `false` (or `nil`) to disable the loader entirely —
every `lookup/1` then returns `nil`. Useful for tests that want to
pin the scorer to the in-source defaults.
## Caching
Reads happen once per file path and cache in `:persistent_term`. A
missing file, malformed JSON, or disabled config all cache an empty
override map so a single failure doesn't cause repeated I/O.
## File schema (v1)
{
"schema_version": 1,
"band_overrides": {
"10000": {
"weights": {
"humidity": 0.20, "time_of_day": 0.04, ...
},
"n_samples": 54161,
"source": "contacts"
}
}
}
Keys in the inner `weights` map are stringified atoms — the loader
converts them with `String.to_existing_atom/1`, which is safe because
`BandConfig` already references every factor at compile time.
"""
require Logger
@cache_key __MODULE__
@factor_keys ~w(humidity time_of_day td_depression refractivity sky season wind rain pwat pressure)a
@doc """
Returns the override weights map for the given band in MHz, or `nil`
if the band has no override (or the loader is disabled / failed).
"""
@spec lookup(pos_integer()) :: map() | nil
def lookup(freq_mhz) when is_integer(freq_mhz) do
Map.get(ensure_loaded(), freq_mhz)
end
@doc """
Returns every override the file defined, keyed by integer freq_mhz.
Mostly for diagnostics — `BandConfig.weights/1` uses `lookup/1`.
"""
@spec all_overrides() :: %{pos_integer() => map()}
def all_overrides, do: ensure_loaded()
@doc """
Clears the persistent_term cache. Tests call this in setup; production
code should never need it because the JSON file is read-only at runtime.
"""
@spec reset() :: :ok
def reset do
_ = :persistent_term.erase(@cache_key)
:ok
end
# ── Lazy loader ─────────────────────────────────────────────────
defp ensure_loaded do
case :persistent_term.get(@cache_key, :not_loaded) do
:not_loaded ->
overrides = load_from_disk()
:persistent_term.put(@cache_key, overrides)
overrides
cached ->
cached
end
end
defp load_from_disk do
case configured_path() do
nil ->
%{}
path when is_binary(path) ->
read_and_parse(path)
end
end
defp configured_path do
case Application.get_env(:microwaveprop, :band_weights_json, default_path()) do
false -> nil
nil -> nil
path when is_binary(path) -> expand(path)
end
end
defp default_path, do: "priv/algo/band_weights.json"
# Resolve a relative path against the app's priv_dir so the file ships
# with a release. Absolute paths pass through unchanged so tests can
# point at a tmp file.
defp expand(path) do
cond do
Path.type(path) == :absolute ->
path
String.starts_with?(path, "priv/") ->
Path.join(:code.priv_dir(:microwaveprop), String.replace_prefix(path, "priv/", ""))
true ->
path
end
end
defp read_and_parse(path) do
with {:ok, body} <- File.read(path),
{:ok, %{} = json} <- Jason.decode(body) do
parse_overrides(json)
else
{:error, :enoent} ->
Logger.info("BandWeights: #{path} not present, using in-source defaults")
%{}
{:error, reason} ->
Logger.warning("BandWeights: failed to load #{path}: #{inspect(reason)}")
%{}
end
end
defp parse_overrides(%{"band_overrides" => overrides}) when is_map(overrides) do
for {band_str, %{"weights" => weights}} <- overrides,
{freq, ""} <- [Integer.parse(band_str)],
parsed = parse_weights(weights),
not is_nil(parsed),
into: %{} do
{freq, parsed}
end
end
defp parse_overrides(_), do: %{}
defp parse_weights(weights) when is_map(weights) do
parsed =
for k <- @factor_keys, into: %{} do
case Map.fetch(weights, Atom.to_string(k)) do
{:ok, v} when is_number(v) -> {k, v / 1}
_ -> {k, nil}
end
end
if Enum.any?(parsed, fn {_, v} -> is_nil(v) end) do
Logger.warning("BandWeights: skipping override missing factors: #{inspect(Map.keys(weights))}")
nil
else
parsed
end
end
defp parse_weights(_), do: nil
end

241
priv/algo/band_weights.json Normal file
View file

@ -0,0 +1,241 @@
{
"generated_at": "2026-05-25T19:44:49+00:00",
"generated_by": "scripts/recalibrate.py",
"schema_version": 1,
"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
},
"data_sources": {
"pskr": {
"bands": [
50,
144,
432,
1296,
2304,
10000
],
"total_samples": 411619
},
"contacts": {
"bands": [
222,
432,
902,
1296,
2304,
3400,
5760,
10000,
24000,
47000,
75000,
122000
],
"total_samples": 81571
}
},
"band_overrides": {
"50": {
"weights": {
"humidity": 0.1252,
"time_of_day": 0.0273,
"td_depression": 0.1511,
"refractivity": 0.0764,
"sky": 0.1006,
"season": 0.2171,
"wind": 0.1006,
"rain": 0.0171,
"pwat": 0.1197,
"pressure": 0.0649
},
"n_samples": 396474,
"source": "pskr"
},
"144": {
"weights": {
"humidity": 0.1865,
"time_of_day": 0.0197,
"td_depression": 0.1492,
"refractivity": 0.1106,
"sky": 0.0621,
"season": 0.1173,
"wind": 0.0621,
"rain": 0.0106,
"pwat": 0.1695,
"pressure": 0.1125
},
"n_samples": 13710,
"source": "pskr"
},
"222": {
"weights": {
"humidity": 0.2002,
"time_of_day": 0.0271,
"td_depression": 0.1804,
"refractivity": 0.101,
"sky": 0.0751,
"season": 0.1316,
"wind": 0.0751,
"rain": 0.0128,
"pwat": 0.1535,
"pressure": 0.0432
},
"n_samples": 7518,
"source": "contacts"
},
"432": {
"weights": {
"humidity": 0.153,
"time_of_day": 0.0256,
"td_depression": 0.17,
"refractivity": 0.1659,
"sky": 0.0708,
"season": 0.1145,
"wind": 0.0708,
"rain": 0.012,
"pwat": 0.1376,
"pressure": 0.0799
},
"n_samples": 1013,
"source": "pskr"
},
"902": {
"weights": {
"humidity": 0.2204,
"time_of_day": 0.0299,
"td_depression": 0.1764,
"refractivity": 0.0967,
"sky": 0.0734,
"season": 0.1089,
"wind": 0.0734,
"rain": 0.0125,
"pwat": 0.1273,
"pressure": 0.081
},
"n_samples": 1754,
"source": "contacts"
},
"1296": {
"weights": {
"humidity": 0.2259,
"time_of_day": 0.034,
"td_depression": 0.1528,
"refractivity": 0.1429,
"sky": 0.0753,
"season": 0.1015,
"wind": 0.0753,
"rain": 0.0128,
"pwat": 0.1327,
"pressure": 0.0468
},
"n_samples": 2951,
"source": "contacts"
},
"2304": {
"weights": {
"humidity": 0.1886,
"time_of_day": 0.0284,
"td_depression": 0.1503,
"refractivity": 0.1474,
"sky": 0.0629,
"season": 0.0848,
"wind": 0.0629,
"rain": 0.0338,
"pwat": 0.1714,
"pressure": 0.0695
},
"n_samples": 719,
"source": "contacts"
},
"3400": {
"weights": {
"humidity": 0.2001,
"time_of_day": 0.0301,
"td_depression": 0.1601,
"refractivity": 0.1102,
"sky": 0.0667,
"season": 0.0899,
"wind": 0.0667,
"rain": 0.0507,
"pwat": 0.1818,
"pressure": 0.0437
},
"n_samples": 335,
"source": "contacts"
},
"5760": {
"weights": {
"humidity": 0.1753,
"time_of_day": 0.0299,
"td_depression": 0.1589,
"refractivity": 0.0491,
"sky": 0.0661,
"season": 0.0892,
"wind": 0.0661,
"rain": 0.0796,
"pwat": 0.1708,
"pressure": 0.1151
},
"n_samples": 321,
"source": "contacts"
},
"10000": {
"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
},
"n_samples": 54161,
"source": "contacts"
},
"24000": {
"weights": {
"humidity": 0.1483,
"time_of_day": 0.0402,
"td_depression": 0.1187,
"refractivity": 0.0517,
"sky": 0.0494,
"season": 0.0733,
"wind": 0.0494,
"rain": 0.2224,
"pwat": 0.1348,
"pressure": 0.112
},
"n_samples": 3808,
"source": "contacts"
},
"47000": {
"weights": {
"humidity": 0.1403,
"time_of_day": 0.0528,
"td_depression": 0.1123,
"refractivity": 0.0517,
"sky": 0.0468,
"season": 0.0757,
"wind": 0.0468,
"rain": 0.2387,
"pwat": 0.1275,
"pressure": 0.1075
},
"n_samples": 759,
"source": "contacts"
}
}
}

281
scripts/derive_band_weights.py
View file

@ -1,281 +0,0 @@
#!/usr/bin/env python3
"""
Per-band weight derivation from the full-corpus correlation analysis.
Reads Pearson correlations that `recalibrate_algo.py` wrote into the DB
and emits a `@band_weight_overrides` Elixir map ready to drop into
`lib/microwaveprop/propagation/band_config.ex`.
The derivation rule is deliberately simple and legible, because the
correlations are noisy (|r| typically 0.05–0.25) and we'd rather preserve
the globally-fit default weights as a prior than over-fit per-band.
Rule:
1. Start from the default weight vector (gradient-descent fit, global).
2. Compute a per-factor "relative signal strength" ratio
s_band,factor = |r_band,factor| / |r_10GHz,factor|
with a small ε floor so zero-correlation factors don't collapse.
3. Clamp the ratio to [0.4, 2.5] — noisy per-band correlations shouldn't
move a factor more than 2.5× or less than 0.4× the prior.
4. Factors we don't have direct correlations for (sky, wind, rain,
season, time_of_day) carry a physics-informed multiplier per band:
* rain scales with `rain_k` (ITU-R P.838 specific attenuation
coefficient) normalized to 10 GHz's rain_k.
* season scales up at VHF/low-UHF where Es-like physics amplifies
monthly variation (we don't model Es directly).
* time_of_day scales with frequency (Finding 8 in algo.md).
* sky/wind stay flat — no per-band evidence either direction.
5. Multiply default weights by the multipliers, then normalize to 1.0.
The output goes to stdout as Elixir source ready to paste into BandConfig.
"""
from __future__ import annotations
import argparse
import os
import sys
from dataclasses import dataclass
try:
import psycopg
from psycopg.rows import dict_row
except ImportError:
sys.exit("psycopg is required: pip install 'psycopg[binary]>=3.1'")
# Default weights (must match `@weights` in band_config.ex). These are the
# globally-fit gradient-descent priors — the per-band deltas modulate them.
DEFAULT_WEIGHTS = {
"rain": 0.1362,
"humidity": 0.1243,
"pwat": 0.1128,
"season": 0.1112,
"refractivity": 0.1049,
"pressure": 0.1032,
"td_depression": 0.0978,
"sky": 0.08,
"wind": 0.08,
"time_of_day": 0.0496,
}
# Correlation-fed factor → source-field mapping. These factors get data-
# driven per-band multipliers; the rest get physics-informed multipliers.
FACTOR_SOURCE = {
"humidity": "dpc", # dewpoint carries the moisture signal
"td_depression": "tc", # temperature dominates T-Td variance
"refractivity": "grad", # refractivity gradient
"pressure": "pr",
"pwat": "pwat",
}
# ITU-R P.838-3 rain_k values (per BandConfig). Used for `rain` weight
# scaling — rain attenuation grows by ~4 orders of magnitude from VHF
# to sub-mm, so `rain` weight should track that.
BAND_RAIN_K = {
50: 0.0, 144: 0.0, 222: 0.0, 432: 0.0,
902: 0.0, 1_296: 0.0,
2_304: 0.001, 3_400: 0.002, 5_760: 0.005,
10_000: 0.010, 24_000: 0.070, 47_000: 0.187,
68_000: 0.310, 75_000: 0.345, 122_000: 0.498,
134_000: 0.520, 142_000: 0.530, 145_000: 0.535,
241_000: 0.550, 288_000: 0.560, 322_000: 0.570,
403_000: 0.580, 411_000: 0.580,
}
# Season weight multiplier by band. VHF has large seasonal swings driven
# by tropo-mixing physics (summer peak). Microwave seasonal swings are
# contest-schedule artefacts but the physics (summer convective mixing
# hurting 24+ GHz) still applies. No per-band evidence in the correlation
# report so we use physics priors.
BAND_SEASON_MULT = {
50: 1.6, 144: 1.4, 222: 1.3, 432: 1.2,
902: 1.1, 1_296: 1.0,
2_304: 1.0, 3_400: 1.0, 5_760: 1.0,
10_000: 1.0, 24_000: 1.1, 47_000: 1.2,
68_000: 1.3, 75_000: 1.4, 122_000: 1.5, 134_000: 1.5,
142_000: 1.5, 145_000: 1.5, 241_000: 1.6, 288_000: 1.6,
322_000: 1.7, 403_000: 1.7, 411_000: 1.7,
}
# Time-of-day multiplier by band (Finding 8 in algo.md — effect scales
# with frequency). 10 GHz is the baseline.
BAND_TOD_MULT = {
50: 0.6, 144: 0.7, 222: 0.8, 432: 0.8,
902: 0.9, 1_296: 1.0,
2_304: 1.0, 3_400: 1.0, 5_760: 1.0,
10_000: 1.0, 24_000: 1.8, 47_000: 2.5,
68_000: 3.5, 75_000: 4.0,
122_000: 5.0, 134_000: 5.0, 142_000: 5.0, 145_000: 5.0,
241_000: 5.0, 288_000: 5.0, 322_000: 5.0, 403_000: 5.0,
411_000: 5.0,
}
# EPS sets an absolute floor on |r| so two near-zero correlations don't
# produce a giant ratio. Set at the noise floor we can reliably distinguish
# from zero given typical n. Tightened along with the clamp so small
# corpora don't drag the weights around.
EPS = 0.05
CLAMP_LO = 0.5
CLAMP_HI = 2.0
MIN_N_FOR_FIT = 200 # bands with fewer matched contacts use physics only
@dataclass
class BandCorr:
band_mhz: int
n: int
rho: dict[str, float]
PER_BAND_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
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 band,
count(*) AS n,
CORR(dist, tc) AS rho_tc,
CORR(dist, dpc) AS rho_dpc,
CORR(dist, pr) AS rho_pr,
CORR(dist, pwat) AS rho_pwat,
CORR(dist, grad) AS rho_grad
FROM joined
WHERE band >= 50
GROUP BY band
HAVING count(*) >= 50
ORDER BY band;
"""
def load_correlations(dsn: str) -> dict[int, BandCorr]:
with psycopg.connect(dsn) as conn:
with conn.cursor(row_factory=dict_row) as cur:
cur.execute("SET statement_timeout = '30min'")
cur.execute(PER_BAND_SQL)
rows = cur.fetchall()
out = {}
for r in rows:
out[int(r["band"])] = BandCorr(
band_mhz=int(r["band"]),
n=int(r["n"]),
rho={
"tc": r["rho_tc"] or 0.0,
"dpc": r["rho_dpc"] or 0.0,
"pr": r["rho_pr"] or 0.0,
"pwat": r["rho_pwat"] or 0.0,
"grad": r["rho_grad"] or 0.0,
},
)
return out
def derive_weights(
band_mhz: int, corrs: dict[int, BandCorr]
) -> dict[str, float] | None:
"""Returns a normalized weight map or None if the band should use defaults."""
band_corr = corrs.get(band_mhz)
ref_corr = corrs.get(10_000)
if band_corr is None or band_corr.n < MIN_N_FOR_FIT or ref_corr is None:
return None
multipliers = {}
# Data-driven multipliers for correlation-backed factors.
for factor, field in FACTOR_SOURCE.items():
band_r = abs(band_corr.rho[field]) + EPS
ref_r = abs(ref_corr.rho[field]) + EPS
ratio = band_r / ref_r
multipliers[factor] = max(CLAMP_LO, min(CLAMP_HI, ratio))
# Physics-driven multipliers for the rest.
# Rain: weight tracks rain_k but with sqrt-dampening. A 7× rain_k jump
# (10→24 GHz) should not become a 7× weight jump — the score itself
# already scales with rain_k via ITU-R P.838 in `score_rain`, so the
# weight only needs to capture how often rain drives the outcome.
rain_k = BAND_RAIN_K.get(band_mhz, 0.010)
rain_k_10g = BAND_RAIN_K[10_000]
if rain_k > 0:
ratio = rain_k / rain_k_10g
rain_mult = max(0.2, min(3.0, ratio**0.5))
else:
rain_mult = 0.1
multipliers["rain"] = rain_mult
multipliers["season"] = BAND_SEASON_MULT.get(band_mhz, 1.0)
multipliers["time_of_day"] = BAND_TOD_MULT.get(band_mhz, 1.0)
multipliers["sky"] = 1.0
multipliers["wind"] = 1.0
# Apply and normalize.
raw = {k: DEFAULT_WEIGHTS[k] * multipliers[k] for k in DEFAULT_WEIGHTS}
total = sum(raw.values())
return {k: round(v / total, 4) for k, v in raw.items()}
def format_elixir(band_mhz: int, weights: dict[str, float]) -> str:
# Preserve the factor ordering used in band_config.ex comments.
order = [
"humidity", "time_of_day", "td_depression", "refractivity",
"sky", "season", "wind", "rain", "pwat", "pressure",
]
lines = [f" {k}: {weights[k]:.4f}" for k in order]
inner = ",\n".join(lines)
return f" weights: %{{\n{inner}\n }}"
def main() -> int:
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--dsn", default=os.environ.get("PROP_DB_URL"))
args = parser.parse_args()
if not args.dsn:
print("error: PROP_DB_URL not set and --dsn not given", file=sys.stderr)
return 2
print(f"loading correlations from {args.dsn}", file=sys.stderr)
corrs = load_correlations(args.dsn)
print(
f"correlations available for {len(corrs)} bands: "
f"{sorted(corrs.keys())}",
file=sys.stderr,
)
# Every band in BandConfig — we emit an entry for each so the reader
# sees at a glance which bands got fit and which inherited defaults.
all_bands = sorted(set(BAND_RAIN_K.keys()) | set(corrs.keys()))
sections = []
for band in all_bands:
weights = derive_weights(band, corrs)
band_corr = corrs.get(band)
n = band_corr.n if band_corr else 0
if weights is None:
sections.append(
f"# {band} MHz — n={n}, uses default @weights (insufficient signal)"
)
else:
sections.append(
f"# {band} MHz — n={n}, derived from corpus correlations\n"
f"{format_elixir(band, weights)}"
)
print("\n\n".join(sections))
return 0
if __name__ == "__main__":
raise SystemExit(main())

629
scripts/recalibrate.py Executable file
View file

@ -0,0 +1,629 @@
#!/usr/bin/env python3
"""
Full algorithm recalibration — unified Python pipeline.
Recomputes per-band weight overrides for the 10-factor scorer from real
propagation evidence in the prod DB:
* VHF / UHF (50, 144, 432 MHz): Pearson correlations between PSKR
`spot_count` density and HRRR features in `pskr_calibration_samples`
(~520k samples since 2026-05-04 in prod).
* Microwave (>=902 MHz): Pearson correlations between QSO `distance_km`
and HRRR features in `contacts` joined to `hrrr_profiles`.
Weight derivation (per band):
1. Start from the globally-fit default weight vector (`DEFAULT_WEIGHTS`).
2. Per correlation-backed factor, build a "relative signal strength"
multiplier = (|band_r| + EPS) / (|ref_r| + EPS), clamped to
[CLAMP_LO, CLAMP_HI]. The reference is the 10 GHz contacts fit
because it has the densest evidence in the corpus.
3. Factors without per-band correlations (rain, season, time_of_day,
sky, wind) get physics-informed multipliers — same priors used by
the legacy `derive_band_weights.py`.
4. Multiply defaults by multipliers, normalize to sum = 1.0.
Outputs (no Elixir source ever modified):
* priv/algo/band_weights.json — machine-readable per-band overrides
that BandConfig can load at runtime via Jason.
* docs/algo-reports/{today}-full-recalibration.md — human report
with per-band correlations, sample counts, derived weights, and
a "what moved" diff against the previous JSON file (if present).
Usage:
direnv allow # picks up PROP_PROD_DB_URL
python3 scripts/recalibrate.py # writes JSON + report
python3 scripts/recalibrate.py --dry-run # prints JSON to stdout
The script is self-contained: psycopg + stdlib only (no pandas, no
numpy). Run it whenever the contact or PSKR corpus grows materially.
"""
from __future__ import annotations
import argparse
import datetime as dt
import json
import os
import sys
from dataclasses import dataclass, field
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
DEFAULT_JSON_OUT = REPO_ROOT / "priv" / "algo" / "band_weights.json"
DEFAULT_REPORT_DIR = REPO_ROOT / "docs" / "algo-reports"
# ── Constants ────────────────────────────────────────────────────────────────
# Must match `@weights` in lib/microwaveprop/propagation/band_config.ex.
# These are the globally-fit gradient-descent priors (2026-04-11) that
# every per-band fit modulates.
DEFAULT_WEIGHTS: dict[str, float] = {
"humidity": 0.1262,
"time_of_day": 0.0380,
"td_depression": 0.1010,
"refractivity": 0.0986,
"sky": 0.0841,
"season": 0.1134,
"wind": 0.0841,
"rain": 0.1431,
"pwat": 0.1147,
"pressure": 0.0967,
}
# Correlation-fed factor → HRRR field. These factors get data-driven
# multipliers per band. The rest get physics priors.
FACTOR_SOURCE: dict[str, str] = {
"humidity": "dpc", # dewpoint carries the moisture signal
"td_depression": "tc", # temperature dominates T-Td variance
"refractivity": "grad", # refractivity gradient
"pressure": "pr",
"pwat": "pwat",
}
# ITU-R P.838-3 rain_k (per BandConfig). Used to scale the `rain` weight
# by a sqrt-dampened ratio against 10 GHz (the score itself already
# scales by rain_k internally).
BAND_RAIN_K: dict[int, float] = {
50: 0.0, 144: 0.0, 222: 0.0, 432: 0.0,
902: 0.0, 1_296: 0.0,
2_304: 0.001, 3_400: 0.002, 5_760: 0.005,
10_000: 0.010, 24_000: 0.070, 47_000: 0.187,
68_000: 0.310, 75_000: 0.345,
122_000: 0.498, 134_000: 0.520, 142_000: 0.530, 145_000: 0.535,
241_000: 0.550, 288_000: 0.560, 322_000: 0.570,
403_000: 0.580, 411_000: 0.580,
}
# Season weight multiplier (physics prior — see algo.md and the legacy
# derive_band_weights.py). VHF has large seasonal swings (Es / tropo
# mixing); 24+ GHz peaks in winter due to dry-air absorption.
BAND_SEASON_MULT: dict[int, float] = {
50: 1.6, 144: 1.4, 222: 1.3, 432: 1.2,
902: 1.1, 1_296: 1.0,
2_304: 1.0, 3_400: 1.0, 5_760: 1.0,
10_000: 1.0, 24_000: 1.1, 47_000: 1.2,
68_000: 1.3, 75_000: 1.4,
122_000: 1.5, 134_000: 1.5, 142_000: 1.5, 145_000: 1.5,
241_000: 1.6, 288_000: 1.6, 322_000: 1.7, 403_000: 1.7, 411_000: 1.7,
}
# Time-of-day multiplier — Finding 8 in algo.md (effect scales with
# frequency). 10 GHz is the baseline at 1.0.
BAND_TOD_MULT: dict[int, float] = {
50: 0.6, 144: 0.7, 222: 0.8, 432: 0.8,
902: 0.9, 1_296: 1.0,
2_304: 1.0, 3_400: 1.0, 5_760: 1.0,
10_000: 1.0, 24_000: 1.8, 47_000: 2.5,
68_000: 3.5, 75_000: 4.0,
122_000: 5.0, 134_000: 5.0, 142_000: 5.0, 145_000: 5.0,
241_000: 5.0, 288_000: 5.0, 322_000: 5.0, 403_000: 5.0, 411_000: 5.0,
}
# PSKR `band` strings → MHz (PSKR scraper stores ham band labels).
PSKR_BAND_TO_MHZ: dict[str, int] = {
"6m": 50,
"2m": 144,
"70cm": 432,
"23cm": 1_296,
"13cm": 2_304,
"3cm": 10_000,
}
# Per-band data source routing. PSKR for VHF/UHF where spot density is
# rich and contacts are sparse; contacts for microwave where PSKR
# activity is rare and the contact corpus is dense.
PSKR_BANDS = {50, 144, 432}
# Floor on |r| (noise level we can distinguish from zero given typical n).
# Same as legacy derive_band_weights.py for continuity.
EPS = 0.05
CLAMP_LO = 0.5
CLAMP_HI = 2.0
MIN_N_FOR_FIT = 200
# Factor display order in JSON + Markdown. Matches band_config.ex
# comment ordering for diffability.
FACTOR_ORDER = [
"humidity", "time_of_day", "td_depression", "refractivity",
"sky", "season", "wind", "rain", "pwat", "pressure",
]
# ── Data structures ─────────────────────────────────────────────────────────
@dataclass
class BandCorr:
"""Per-band correlation snapshot, source-agnostic."""
band_mhz: int
source: str # "pskr" or "contacts"
n: int
rho: dict[str, float] = field(default_factory=dict)
# Optional per-band diagnostics for the report.
total_label: float | None = None
avg_label: float | None = None
# ── SQL ─────────────────────────────────────────────────────────────────────
# Contacts ↔ HRRR (microwave): per-band Pearson(distance_km, feature).
# Same nearest-neighbour join as the legacy recalibrate_algo.py.
CONTACTS_CORR_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
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 band,
count(*)::int AS n,
CORR(dist, tc) AS rho_tc,
CORR(dist, dpc) AS rho_dpc,
CORR(dist, pr) AS rho_pr,
CORR(dist, pwat) AS rho_pwat,
CORR(dist, grad) AS rho_grad,
avg(dist)::float AS avg_dist
FROM joined
WHERE band >= 50
GROUP BY band
HAVING count(*) >= 50
ORDER BY band;
"""
# PSKR samples already carry HRRR fields (the PSKR sampler joins them at
# ingestion). Label is `spot_count` per (hour, midpoint, band) cell.
# Filter to HRRR-complete rows so all five correlations are computed
# from the same denominator.
PSKR_CORR_SQL = """
SELECT band,
count(*)::int AS n,
CORR(spot_count, surface_temp_c) AS rho_tc,
CORR(spot_count, surface_dewpoint_c) AS rho_dpc,
CORR(spot_count, surface_pressure_mb) AS rho_pr,
CORR(spot_count, pwat_mm) AS rho_pwat,
CORR(spot_count, min_refractivity_gradient) AS rho_grad,
sum(spot_count)::bigint AS total_spots,
avg(spot_count)::float AS avg_spots
FROM pskr_calibration_samples
WHERE surface_pressure_mb IS NOT NULL
AND min_refractivity_gradient IS NOT NULL
GROUP BY band
ORDER BY band;
"""
# ── Loaders ─────────────────────────────────────────────────────────────────
def _row_to_rho(row: dict) -> dict[str, float]:
return {
"tc": row["rho_tc"] or 0.0,
"dpc": row["rho_dpc"] or 0.0,
"pr": row["rho_pr"] or 0.0,
"pwat": row["rho_pwat"] or 0.0,
"grad": row["rho_grad"] or 0.0,
}
def load_contacts_correlations(conn) -> dict[int, BandCorr]:
print("• loading contacts ↔ HRRR correlations", file=sys.stderr)
with conn.cursor(row_factory=dict_row) as cur:
cur.execute(CONTACTS_CORR_SQL)
rows = cur.fetchall()
out: dict[int, BandCorr] = {}
for r in rows:
band = int(r["band"])
out[band] = BandCorr(
band_mhz=band,
source="contacts",
n=int(r["n"]),
rho=_row_to_rho(r),
avg_label=r.get("avg_dist"),
)
return out
def load_pskr_correlations(conn) -> dict[int, BandCorr]:
print("• loading PSKR ↔ HRRR correlations", file=sys.stderr)
with conn.cursor(row_factory=dict_row) as cur:
cur.execute(PSKR_CORR_SQL)
rows = cur.fetchall()
out: dict[int, BandCorr] = {}
for r in rows:
band_label = r["band"]
mhz = PSKR_BAND_TO_MHZ.get(band_label)
if mhz is None:
# Unknown PSKR band string — log and skip rather than guess.
print(f" (skipping unknown PSKR band {band_label!r})", file=sys.stderr)
continue
out[mhz] = BandCorr(
band_mhz=mhz,
source="pskr",
n=int(r["n"]),
rho=_row_to_rho(r),
total_label=float(r["total_spots"]) if r["total_spots"] is not None else None,
avg_label=r.get("avg_spots"),
)
return out
def merge_correlations(
pskr: dict[int, BandCorr],
contacts: dict[int, BandCorr],
) -> dict[int, BandCorr]:
"""Per-band routing: PSKR for VHF/UHF, contacts for everything else."""
out: dict[int, BandCorr] = {}
for band in set(pskr) | set(contacts):
if band in PSKR_BANDS and band in pskr:
out[band] = pskr[band]
elif band in contacts:
out[band] = contacts[band]
elif band in pskr:
# No contacts but PSKR has it (e.g. 23cm/13cm if we expand routing).
out[band] = pskr[band]
return out
# ── Weight derivation ──────────────────────────────────────────────────────
def derive_weights(
band_mhz: int,
corrs: dict[int, BandCorr],
) -> tuple[dict[str, float] | None, dict[str, float]]:
"""Returns (weights, multipliers). weights=None means inherit defaults.
`multipliers` is always populated so the report can show why a band
moved (or why it didn't qualify for a fit).
"""
band_corr = corrs.get(band_mhz)
# Reference is always the 10 GHz contacts fit — densest, most trusted.
ref_corr = corrs.get(10_000)
multipliers: dict[str, float] = {}
# Correlation-fed factors first (skipped if either side missing).
if band_corr is not None and ref_corr is not None:
for factor, field_name in FACTOR_SOURCE.items():
band_r = abs(band_corr.rho[field_name]) + EPS
ref_r = abs(ref_corr.rho[field_name]) + EPS
ratio = band_r / ref_r
multipliers[factor] = max(CLAMP_LO, min(CLAMP_HI, ratio))
else:
for factor in FACTOR_SOURCE:
multipliers[factor] = 1.0
# Physics priors for the rest.
rain_k = BAND_RAIN_K.get(band_mhz, 0.010)
rain_k_ref = BAND_RAIN_K[10_000]
if rain_k > 0:
multipliers["rain"] = max(0.2, min(3.0, (rain_k / rain_k_ref) ** 0.5))
else:
multipliers["rain"] = 0.1
multipliers["season"] = BAND_SEASON_MULT.get(band_mhz, 1.0)
multipliers["time_of_day"] = BAND_TOD_MULT.get(band_mhz, 1.0)
multipliers["sky"] = 1.0
multipliers["wind"] = 1.0
# No fit if we don't have enough data on this band.
if band_corr is None or band_corr.n < MIN_N_FOR_FIT:
return None, multipliers
raw = {k: DEFAULT_WEIGHTS[k] * multipliers[k] for k in DEFAULT_WEIGHTS}
total = sum(raw.values())
weights = {k: round(v / total, 4) for k, v in raw.items()}
return weights, multipliers
# ── Output ─────────────────────────────────────────────────────────────────
def build_json(
corrs: dict[int, BandCorr],
pskr_only: dict[int, BandCorr],
contacts_only: dict[int, BandCorr],
generated_at: dt.datetime,
) -> dict:
"""Stable JSON shape — keys sorted, numbers rounded for diffability."""
all_bands = sorted(set(BAND_RAIN_K.keys()) | set(corrs.keys()))
overrides: dict[str, dict] = {}
for band in all_bands:
weights, _mults = derive_weights(band, corrs)
if weights is None:
continue
bc = corrs.get(band)
overrides[str(band)] = {
"weights": {k: weights[k] for k in FACTOR_ORDER},
"n_samples": bc.n if bc else 0,
"source": bc.source if bc else "default",
}
return {
"generated_at": generated_at.isoformat(timespec="seconds"),
"generated_by": "scripts/recalibrate.py",
"schema_version": 1,
"global_weights": {k: DEFAULT_WEIGHTS[k] for k in FACTOR_ORDER},
"data_sources": {
"pskr": {
"bands": sorted(pskr_only.keys()),
"total_samples": sum(bc.n for bc in pskr_only.values()),
},
"contacts": {
"bands": sorted(contacts_only.keys()),
"total_samples": sum(bc.n for bc in contacts_only.values()),
},
},
"band_overrides": overrides,
}
def load_prior_json(path: Path) -> dict | None:
if not path.exists():
return None
try:
return json.loads(path.read_text())
except (OSError, json.JSONDecodeError) as exc:
print(f" (could not read prior JSON at {path}: {exc})", file=sys.stderr)
return None
def diff_weights(
prior: dict | None, current: dict,
) -> list[tuple[str, str, float, float, float]]:
"""Returns rows (band, factor, prior_w, new_w, delta) sorted by |delta|."""
if not prior:
return []
prior_overrides = prior.get("band_overrides", {})
rows: list[tuple[str, str, float, float, float]] = []
for band, entry in current.get("band_overrides", {}).items():
prior_entry = prior_overrides.get(band, {}).get("weights", {})
new_entry = entry["weights"]
for factor in FACTOR_ORDER:
old = float(prior_entry.get(factor, DEFAULT_WEIGHTS[factor]))
new = float(new_entry[factor])
delta = new - old
if abs(delta) >= 0.005: # 0.5pp threshold to keep noise out
rows.append((band, factor, old, new, delta))
rows.sort(key=lambda r: abs(r[4]), reverse=True)
return rows
def render_report(
payload: dict,
corrs: dict[int, BandCorr],
diffs: list[tuple[str, str, float, float, float]],
dsn_redacted: str,
) -> str:
today = payload["generated_at"][:10]
lines: list[str] = []
lines.append(f"# Full Recalibration Report — {today}\n")
lines.append(
"> Auto-generated by `scripts/recalibrate.py`. The companion "
"JSON file at `priv/algo/band_weights.json` is the machine-"
"readable output that Elixir loads at runtime. This report is "
"the human-readable audit trail.\n"
)
lines.append(f"Connection: `{dsn_redacted}`")
lines.append(f"Generated: `{payload['generated_at']}`\n")
# ── Data sources
lines.append("## Data sources\n")
src = payload["data_sources"]
lines.append(
f"- **PSKR** (VHF/UHF, label = `spot_count`): bands "
f"{src['pskr']['bands']}, {src['pskr']['total_samples']:,} samples"
)
lines.append(
f"- **Contacts ↔ HRRR** (microwave, label = `distance_km`): "
f"{len(src['contacts']['bands'])} bands, "
f"{src['contacts']['total_samples']:,} matched contacts"
)
lines.append("")
# ── Per-band evidence summary
lines.append("## Per-band evidence\n")
lines.append("| Band | Source | n | label avg | ρ_tc | ρ_dpc | ρ_pr | ρ_pwat | ρ_grad |")
lines.append("|------|--------|---:|----------:|-----:|------:|-----:|-------:|-------:|")
for band in sorted(corrs):
bc = corrs[band]
avg = f"{bc.avg_label:.2f}" if bc.avg_label is not None else "—"
lines.append(
f"| {band} MHz | {bc.source} | {bc.n:,} | {avg} | "
f"{bc.rho['tc']:+.3f} | {bc.rho['dpc']:+.3f} | {bc.rho['pr']:+.3f} | "
f"{bc.rho['pwat']:+.3f} | {bc.rho['grad']:+.3f} |"
)
lines.append("")
# ── Derived overrides
lines.append("## Derived per-band weights\n")
overrides = payload["band_overrides"]
if not overrides:
lines.append("_No band had ≥{} samples — nothing to override._\n".format(MIN_N_FOR_FIT))
else:
header = "| Band | n | source | " + " | ".join(FACTOR_ORDER) + " |"
sep = "|------|---:|--------|" + "|".join(["----:"] * len(FACTOR_ORDER)) + "|"
lines.append(header)
lines.append(sep)
for band_str in sorted(overrides, key=int):
entry = overrides[band_str]
ws = entry["weights"]
cells = " | ".join(f"{ws[f]:.4f}" for f in FACTOR_ORDER)
lines.append(
f"| {band_str} MHz | {entry['n_samples']:,} | "
f"{entry['source']} | {cells} |"
)
lines.append("")
# ── Diff against prior
lines.append("## Changes vs prior `band_weights.json`\n")
if not diffs:
lines.append(
"_No prior file, or no factor moved by ≥0.005. First-run "
"values can be compared to the in-tree `@band_configs` "
"weights in `lib/microwaveprop/propagation/band_config.ex`._\n"
)
else:
lines.append("Largest weight moves (≥0.5 pp), sorted by magnitude:\n")
lines.append("| Band | Factor | Prior | New | Δ |")
lines.append("|------|--------|------:|----:|--:|")
for band, factor, old, new, delta in diffs[:30]:
lines.append(
f"| {band} MHz | {factor} | {old:.4f} | {new:.4f} | {delta:+.4f} |"
)
if len(diffs) > 30:
lines.append(f"\n_(+{len(diffs) - 30} smaller moves omitted)_")
lines.append("")
# ── How to apply
lines.append("## Applying these weights\n")
lines.append(
"This script writes `priv/algo/band_weights.json` and stops "
"there — it never edits `band_config.ex`. To make the weights "
"live, `BandConfig.weights(band_mhz)` needs a one-time edit to "
"read this JSON at app start (cache in `:persistent_term`) and "
"merge `band_overrides[band].weights` over the in-source "
"defaults. Until then, the JSON is documentation only."
)
return "\n".join(lines)
# ── Plumbing ───────────────────────────────────────────────────────────────
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
def main() -> int:
parser = argparse.ArgumentParser(
description=__doc__,
formatter_class=argparse.RawDescriptionHelpFormatter,
)
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-json",
type=Path,
default=DEFAULT_JSON_OUT,
help=f"JSON output path (default: {DEFAULT_JSON_OUT.relative_to(REPO_ROOT)})",
)
parser.add_argument(
"--out-report",
type=Path,
default=None,
help="Markdown report path (default: docs/algo-reports/{today}-full-recalibration.md)",
)
parser.add_argument(
"--dry-run",
action="store_true",
help="Print JSON to stdout, do not write files",
)
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",
file=sys.stderr,
)
return 2
today = dt.date.today().isoformat()
out_report = args.out_report or (DEFAULT_REPORT_DIR / f"{today}-full-recalibration.md")
print(f"connecting to {redact(args.dsn)}", file=sys.stderr)
with psycopg.connect(args.dsn, autocommit=True) as conn:
with conn.cursor() as cur:
cur.execute(f"SET statement_timeout = '{args.statement_timeout}'")
pskr = load_pskr_correlations(conn)
contacts = load_contacts_correlations(conn)
corrs = merge_correlations(pskr, contacts)
print(
f"• merged: {len(corrs)} bands "
f"({sum(1 for b in corrs if corrs[b].source == 'pskr')} PSKR, "
f"{sum(1 for b in corrs if corrs[b].source == 'contacts')} contacts)",
file=sys.stderr,
)
generated_at = dt.datetime.now(dt.timezone.utc).replace(microsecond=0)
payload = build_json(corrs, pskr, contacts, generated_at)
prior = load_prior_json(args.out_json)
diffs = diff_weights(prior, payload)
if args.dry_run:
print(json.dumps(payload, indent=2, sort_keys=False))
print(f"\n(dry-run: would write {len(payload['band_overrides'])} overrides to {args.out_json})", file=sys.stderr)
return 0
args.out_json.parent.mkdir(parents=True, exist_ok=True)
args.out_json.write_text(json.dumps(payload, indent=2, sort_keys=False) + "\n")
print(f"• wrote {args.out_json} ({len(payload['band_overrides'])} overrides)", file=sys.stderr)
out_report.parent.mkdir(parents=True, exist_ok=True)
out_report.write_text(render_report(payload, corrs, diffs, redact(args.dsn)) + "\n")
print(f"• wrote {out_report}", file=sys.stderr)
return 0
if __name__ == "__main__":
raise SystemExit(main())

916
scripts/recalibrate_algo.py
View file

@ -1,916 +0,0 @@
#!/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"),
("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. Each metric is classified in
# DATA_GAP_RULES below so the report can flag BROKEN vs WARN vs OK with
# concrete remediation hints instead of dumping a flat key/value list.
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 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;
"""
# Coverage breakdown for tables whose total row count alone hides the
# real gap — e.g. narr_profiles with 358 rows could be one good year or
# 358 stragglers spread across 30. Year-bucketed counts let the report
# flag whether the backfill ever made meaningful progress.
NARR_COVERAGE_BY_YEAR_SQL = """
SELECT EXTRACT(YEAR FROM valid_time)::int AS year, count(*)::bigint AS profiles
FROM narr_profiles
GROUP BY 1 ORDER BY 1;
"""
# Where is the HRRR enrichment backlog? A 36k pending count could mean
# "the backfill stalled in 2018" or "every contact submitted in the last
# week is queued". The breakdown drives a different remediation: rerun
# HrrrFetchWorker.enqueue_for_qso/1 over a year range vs. unstick the
# Oban queue.
HRRR_STATUS_BY_YEAR_SQL = """
SELECT EXTRACT(YEAR FROM qso_timestamp)::int AS year,
count(*) FILTER (WHERE hrrr_status = 'complete')::bigint AS complete,
count(*) FILTER (WHERE hrrr_status = 'queued')::bigint AS queued,
count(*) FILTER (WHERE hrrr_status = 'failed')::bigint AS failed,
count(*) FILTER (WHERE hrrr_status IS NULL)::bigint AS unset
FROM contacts
WHERE pos1 IS NOT NULL AND distance_km < 3000 AND flagged_invalid = false
AND qso_timestamp >= '2014-10-02'
GROUP BY 1 ORDER BY 1;
"""
# Severity rules for each metric in DATA_GAP_SQL.
# key metric name returned by DATA_GAP_SQL
# broken_when lambda(value) → True means BROKEN (red)
# warn_when lambda(value) → True means WARN (yellow). Skipped if BROKEN.
# remediation one-line action the operator should take
# A metric not flagged by either lambda is OK (green).
DATA_GAP_RULES: list[dict] = [
{
"key": "narr_rows",
# NARR feeds pre-2014 enrichment via BackfillEnqueueWorker's
# 30-min cron. NarrFetchWorker dedupes by 0.13° grid, so the
# equilibrium count is ~3 path points × pre2014_contacts after
# spatial collapse — not the contact count itself. Threshold
# follows: BROKEN at zero, WARN below half the expected
# post-dedup ceiling.
"broken_when": lambda v: v == 0,
"warn_when": lambda _v: False,
"remediation": (
"BackfillEnqueueWorker dispatches narr jobs every 30 min; "
"if rows aren't climbing, check NarrFetchWorker for upstream "
"NCEI errors in the Oban dashboard."
),
},
{
"key": "narr_pre2014_rows",
"broken_when": lambda v: v == 0,
"warn_when": lambda _v: False,
"remediation": (
"NARR's purpose is filling the pre-2014-10 gap; without rows here "
"the historical-sanity correlation table is empty. Same self-heal "
"as narr_rows."
),
},
{
"key": "hrrr_climatology_rows",
"broken_when": lambda v: v == 0,
"warn_when": lambda _v: False,
"remediation": (
"Should self-heal within 24h via the daily 02:30Z cron that "
"runs AdminTaskWorker(task=climatology). If still 0 after a "
"day, check Oban dashboard for failed admin jobs."
),
},
{
"key": "hrrr_complete_contacts",
# Healthy is the fraction complete, not the absolute number — handled
# separately via _check_hrrr_completeness below.
"broken_when": lambda _v: False,
"warn_when": lambda _v: False,
"remediation": "",
},
{
"key": "hrrr_pending_contacts",
"broken_when": lambda _v: False,
# Anything over 5k pending is a real backlog worth chasing.
"warn_when": lambda v: v > 5_000,
"remediation": (
"See `HRRR enrichment status by year` table below to locate the "
"backlog, then `mix hrrr.backfill --year YYYY` (or rerun "
"HrrrFetchWorker.enqueue_for_qso/1 over the affected ids)."
),
},
{
"key": "pre2014_contacts",
# pre-2014 contacts existing is informational, not a gap on its own.
"broken_when": lambda _v: False,
"warn_when": lambda _v: False,
"remediation": "",
},
]
# ─── 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 classify_data_gaps(gaps: dict) -> pd.DataFrame:
"""Apply DATA_GAP_RULES to the raw gap counts.
Returns a dataframe with columns: metric, value, severity, remediation.
`severity` is one of "BROKEN", "WARN", "OK" (informational metrics with
no rule are also "OK"). The dataframe order matches DATA_GAP_RULES so
related metrics group naturally in the rendered table.
"""
rows = []
for rule in DATA_GAP_RULES:
key = rule["key"]
if key not in gaps:
continue
value = gaps[key]
if rule["broken_when"](value):
severity = "BROKEN"
elif rule["warn_when"](value):
severity = "WARN"
else:
severity = "OK"
rows.append(
{
"metric": key,
"value": int(value) if value is not None else None,
"severity": severity,
"remediation": rule["remediation"] if severity != "OK" else "",
}
)
# Append any metric DATA_GAP_SQL returned that we don't have a rule for,
# so the report stays exhaustive even if SQL changes ahead of the rules.
handled = {r["key"] for r in DATA_GAP_RULES}
for key, value in gaps.items():
if key in handled:
continue
rows.append(
{
"metric": key,
"value": int(value) if value is not None else None,
"severity": "OK",
"remediation": "",
}
)
return pd.DataFrame(rows)
def hrrr_completeness_pct(gaps: dict) -> float | None:
"""Fraction of eligible contacts with hrrr_status = 'complete'."""
complete = gaps.get("hrrr_complete_contacts") or 0
pending = gaps.get("hrrr_pending_contacts") or 0
total = complete + pending
if total == 0:
return None
return round(100.0 * complete / total, 1)
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()
gap_table = classify_data_gaps(gaps)
completeness = hrrr_completeness_pct(gaps)
broken = gap_table[gap_table["severity"] == "BROKEN"]
for _, row in broken.iterrows():
print(
f" ⚠ BROKEN: {row['metric']}={row['value']} — {row['remediation']}",
file=sys.stderr,
)
sections.append("\n## Data gaps and ingestion health\n")
sections.append(
"Severity legend: **BROKEN** = table is empty and feeds an "
"active scoring path; **WARN** = populated but coverage is "
"low enough to skew per-band statistics; **OK** = within "
"spec or informational. Each non-OK row carries a one-line "
"remediation hint.\n"
)
sections.append(md_table(gap_table))
if completeness is not None:
sections.append(
f"\nHRRR enrichment completeness: **{completeness}%** "
f"({gaps.get('hrrr_complete_contacts', 0):,} complete / "
f"{gaps.get('hrrr_pending_contacts', 0):,} pending).\n"
)
print("• narr coverage by year", file=sys.stderr)
narr_year_df = fetch_df(conn, NARR_COVERAGE_BY_YEAR_SQL)
sections.append("\n### NARR coverage by year\n")
if narr_year_df.empty:
sections.append(
"_narr_profiles is empty. Run `mix narr.backfill` to "
"populate the pre-2014 corpus._\n"
)
else:
sections.append(
"Year-bucketed counts. NARR is supposed to span 1979 → "
"2014-10 at 3-hourly cadence (~2,920 timestamps per year "
"per grid cell). Years far below the others are partial "
"backfills the worker hasn't finished.\n"
)
sections.append(md_table(narr_year_df))
print("• hrrr status by year", file=sys.stderr)
hrrr_year_df = fetch_df(conn, HRRR_STATUS_BY_YEAR_SQL)
sections.append("\n### HRRR enrichment status by year\n")
if hrrr_year_df.empty:
sections.append("_(no contacts with `pos1` set after 2014-10-02)_\n")
else:
sections.append(
"Locate the pending backlog: a year with high `queued`/"
"`failed`/`unset` counts is where HrrrFetchWorker needs "
"to be rerun. Recent years dominating `queued` usually "
"means the live Oban queue is stuck rather than a "
"historical backfill gap.\n"
)
sections.append(md_table(hrrr_year_df))
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())

165
test/microwaveprop/propagation/band_weights_test.exs Normal file
View file

@ -0,0 +1,165 @@
defmodule Microwaveprop.Propagation.BandWeightsTest do
# async: false because BandWeights caches in :persistent_term, which is
# global to the BEAM. Each test resets the cache in setup so ordering
# is irrelevant, but concurrent runs would still see each other.
use ExUnit.Case, async: false
alias Microwaveprop.Propagation.BandWeights
@valid_json """
{
"generated_at": "2026-05-24T18:53:32+00:00",
"schema_version": 1,
"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": {
"10000": {
"weights": {
"humidity": 0.20, "time_of_day": 0.04, "td_depression": 0.10,
"refractivity": 0.10, "sky": 0.08, "season": 0.11,
"wind": 0.08, "rain": 0.14, "pwat": 0.11, "pressure": 0.04
},
"n_samples": 54161,
"source": "contacts"
},
"144": {
"weights": {
"humidity": 0.18, "time_of_day": 0.02, "td_depression": 0.15,
"refractivity": 0.11, "sky": 0.06, "season": 0.12,
"wind": 0.06, "rain": 0.01, "pwat": 0.17, "pressure": 0.12
},
"n_samples": 12927,
"source": "pskr"
}
}
}
"""
setup do
# Cache is global to the BEAM; reset before every test to avoid
# bleed between cases.
BandWeights.reset()
prior_env = Application.get_env(:microwaveprop, :band_weights_json)
on_exit(fn ->
BandWeights.reset()
if prior_env do
Application.put_env(:microwaveprop, :band_weights_json, prior_env)
else
Application.delete_env(:microwaveprop, :band_weights_json)
end
end)
:ok
end
describe "lookup/1 with a valid JSON file" do
setup do
path = Path.join(System.tmp_dir!(), "band_weights_#{System.unique_integer([:positive])}.json")
File.write!(path, @valid_json)
Application.put_env(:microwaveprop, :band_weights_json, path)
on_exit(fn -> File.rm(path) end)
%{path: path}
end
test "returns the override weights for a band present in the file" do
result = BandWeights.lookup(10_000)
assert is_map(result)
assert result.humidity == 0.20
assert result.pressure == 0.04
# All 10 keys present.
expected_keys =
MapSet.new([
:humidity,
:time_of_day,
:td_depression,
:refractivity,
:sky,
:season,
:wind,
:rain,
:pwat,
:pressure
])
assert MapSet.new(Map.keys(result)) == expected_keys
end
test "returns the override weights for a second band" do
result = BandWeights.lookup(144)
assert result.humidity == 0.18
end
test "returns nil for a band not present in the file" do
assert BandWeights.lookup(75_000) == nil
end
test "is cheap on repeat calls (cached)" do
first = BandWeights.lookup(10_000)
# Delete the source file — second call must still succeed from cache.
File.rm!(Application.fetch_env!(:microwaveprop, :band_weights_json))
second = BandWeights.lookup(10_000)
assert first == second
end
end
describe "lookup/1 when the file is missing" do
setup do
Application.put_env(:microwaveprop, :band_weights_json, "/nonexistent/path/band_weights.json")
:ok
end
test "returns nil for every band (does not raise)" do
assert BandWeights.lookup(10_000) == nil
assert BandWeights.lookup(50) == nil
end
end
describe "lookup/1 when the loader is disabled in config" do
setup do
Application.put_env(:microwaveprop, :band_weights_json, false)
:ok
end
test "returns nil for every band" do
assert BandWeights.lookup(10_000) == nil
assert BandWeights.lookup(144) == nil
end
end
describe "lookup/1 with malformed JSON" do
setup do
path = Path.join(System.tmp_dir!(), "band_weights_bad_#{System.unique_integer([:positive])}.json")
File.write!(path, "{ not valid json")
Application.put_env(:microwaveprop, :band_weights_json, path)
on_exit(fn -> File.rm(path) end)
%{path: path}
end
test "returns nil for every band (logs but does not raise)" do
assert BandWeights.lookup(10_000) == nil
end
end
describe "all_overrides/0" do
setup do
path = Path.join(System.tmp_dir!(), "band_weights_#{System.unique_integer([:positive])}.json")
File.write!(path, @valid_json)
Application.put_env(:microwaveprop, :band_weights_json, path)
on_exit(fn -> File.rm(path) end)
:ok
end
test "returns a map of all band overrides keyed by freq_mhz integer" do
result = BandWeights.all_overrides()
assert is_map(result)
assert Map.has_key?(result, 10_000)
assert Map.has_key?(result, 144)
assert result[10_000].humidity == 0.20
end
end
end