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

20-feature model with solar indices, sounding stability, and ducting

Browse source 
Add SFI, Kp max (solar), K-index, lifted index (sounding stability),
and ducting_detected (HRRR) as model features. Training now joins to
solar_indices and nearest sounding (within 6 hours) for both phases.
Model can learn solar/geomagnetic effects if they exist in the data.
This commit is contained in:
Graham McIntire 2026-04-01 09:31:54 -05:00
parent 07558d17eb
commit 9537c97d1d
No known key found for this signature in database
GPG key ID: F4ABF488E6029E59
3 changed files with 168 additions and 99 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

22
lib/microwaveprop/propagation/model.ex
View file

@ -31,7 +31,7 @@ defmodule Microwaveprop.Propagation.Model do
## Architecture
Feed-forward neural network with:
- Input: 15 features (10 atmospheric + 4 cyclical temporal + 1 frequency)
- Input: 20 features (12 atmospheric + 2 solar + 1 geographic + 4 cyclical temporal + 1 frequency)
- Hidden: 2 dense layers with ReLU activation
- Output: 1 sigmoid unit (score 0-1, scaled to 0-100)
@ -40,7 +40,7 @@ defmodule Microwaveprop.Propagation.Model do
baseline performance metrics.
"""
@feature_count 15
@feature_count 20
@models_dir Path.join(:code.priv_dir(:microwaveprop), "models")
@default_path Path.join(@models_dir, "propagation_v1.nx")
@ -225,6 +225,11 @@ defmodule Microwaveprop.Propagation.Model do
:pwat_mm,
:surface_refractivity,
:latitude,
:sfi,
:kp_max,
:ducting_detected,
:k_index,
:lifted_index,
:solar_hour_sin,
:solar_hour_cos,
:month_sin,
@ -248,6 +253,11 @@ defmodule Microwaveprop.Propagation.Model do
pwat = conditions[:pwat_mm] || 20.0
refractivity = conditions[:surface_refractivity] || 320.0
latitude = conditions[:latitude] || 37.0
sfi = conditions[:sfi] || 120.0
kp_max = conditions[:kp_max] || 2.0
ducting = if conditions[:ducting_detected], do: 1.0, else: 0.0
k_index = conditions[:k_index] || 20.0
lifted_index = conditions[:lifted_index] || 0.0
utc_hour = conditions[:utc_hour] || 12
longitude = conditions[:longitude] || -97.0
month = conditions[:month] || 6
@ -261,8 +271,7 @@ defmodule Microwaveprop.Propagation.Model do
hour_rad = 2 * :math.pi() * local_hour / 24
month_rad = 2 * :math.pi() * (month - 1) / 12
# Normalize atmospheric features to ~[0, 1] using physical bounds.
# Without this, features like pressure (~1013) dominate and cause NaN gradients.
# Normalize all features to ~[0, 1] using physical bounds.
[
(temp_c + 40) / 80,
(dewpoint_c + 40) / 80,
@ -274,6 +283,11 @@ defmodule Microwaveprop.Propagation.Model do
pwat / 60.0,
(refractivity - 250) / 120,
(latitude - 25) / 25,
(sfi - 60) / 200,
kp_max / 9.0,
ducting,
(k_index + 10) / 50,
(lifted_index + 10) / 20,
:math.sin(hour_rad),
:math.cos(hour_rad),
:math.sin(month_rad),

221
lib/mix/tasks/propagation_train.ex
View file

@ -2,15 +2,19 @@ defmodule Mix.Tasks.PropagationTrain do
@shortdoc "Train ML model from QSO-HRRR data with algorithm pre-training"
@moduledoc """
Two-phase ML training for propagation prediction:
Two-phase ML training for propagation prediction using all available data:
Phase 1 (pre-train): Learn broad atmospheric patterns from 500K+ algorithm
scores covering all seasons, times, and locations across CONUS.
scores covering all seasons, times, and locations across CONUS. Joined to
solar indices and nearest sounding for stability indices.
Phase 2 (fine-tune): Calibrate against real QSO outcomes (57K+ contacts
matched to HRRR conditions). Target is within-band distance percentile.
matched to HRRR + solar + sounding data). Target is within-band distance percentile.
Features: 15 inputs (10 atmospheric + latitude + 4 cyclical temporal + frequency).
Features (20): atmospheric (temp, dewpoint, pressure, abs humidity, td depression,
refractivity gradient, HPBL, PWAT, surface refractivity), geographic (latitude),
solar (SFI, Kp max), HRRR-derived (ducting), sounding-derived (K-index,
lifted index), temporal (solar hour sin/cos, month sin/cos), frequency (log MHz).
## Usage
@ -26,7 +30,6 @@ defmodule Mix.Tasks.PropagationTrain do
@query_timeout 600_000
# Per-band normalization for QSO distance percentile computation
@band_max_km %{
10_000 => 800.0,
24_000 => 350.0,
@ -67,7 +70,7 @@ defmodule Mix.Tasks.PropagationTrain do
batch_size = Keyword.get(opts, :batch_size, 256)
"=" |> String.duplicate(70) |> IO.puts()
IO.puts("PROPAGATION MODEL TRAINING (Two-Phase)")
IO.puts("PROPAGATION MODEL TRAINING (Two-Phase, 20 features)")
"=" |> String.duplicate(70) |> IO.puts()
IO.puts("Phase 1: #{pretrain_epochs} epochs on #{pretrain_sample} algorithm scores")
IO.puts("Phase 2: #{finetune_epochs} epochs on real QSO data")
@ -75,12 +78,12 @@ defmodule Mix.Tasks.PropagationTrain do
# ── Phase 1: Pre-train on algorithm scores ──────────────────────
IO.puts("=" <> String.duplicate("─", 50))
IO.puts("PHASE 1: Pre-training on algorithm scores")
IO.puts("PHASE 1: Pre-training on algorithm scores + solar + soundings")
IO.puts("=" <> String.duplicate("─", 50))
{pt_features, pt_targets, pt_counts} = load_pretrain_data(pretrain_sample)
pt_n = elem(Nx.shape(pt_features), 0)
IO.puts("Loaded #{pt_n} algorithm score samples:")
IO.puts("Loaded #{pt_n} samples:")
print_band_counts(pt_counts)
{pt_features, pt_targets} = shuffle(pt_features, pt_targets, pt_n)
@ -96,7 +99,7 @@ defmodule Mix.Tasks.PropagationTrain do
# ── Phase 2: Fine-tune on real QSO data ─────────────────────────
IO.puts("\n" <> "=" <> String.duplicate("─", 50))
IO.puts("PHASE 2: Fine-tuning on real QSO-HRRR data")
IO.puts("PHASE 2: Fine-tuning on QSO + HRRR + solar + soundings")
IO.puts("=" <> String.duplicate("─", 50))
{ft_features, ft_targets, ft_counts} = load_qso_data()
@ -130,7 +133,7 @@ defmodule Mix.Tasks.PropagationTrain do
IO.puts("Finished: #{DateTime.to_string(DateTime.utc_now())}")
end
# ── Phase 1 data: algorithm scores from propagation_scores ────────
# ── Phase 1 data: algorithm scores + solar + nearest sounding ─────
defp load_pretrain_data(sample_size) do
bands = BandConfig.all_freqs()
@ -138,6 +141,7 @@ defmodule Mix.Tasks.PropagationTrain do
IO.puts(" Loading algorithm scores (stratified #{per_band}/band)...")
# Join propagation_scores → hrrr_profiles → solar_indices → nearest sounding
band_queries =
Enum.map_join(bands, "\nUNION ALL\n", fn band_mhz ->
"""
@ -145,14 +149,24 @@ defmodule Mix.Tasks.PropagationTrain do
ps.score, ps.band_mhz,
h.surface_temp_c, h.surface_dewpoint_c, h.surface_pressure_mb,
h.min_refractivity_gradient, h.hpbl_m, h.pwat_mm,
h.surface_refractivity,
ps.lat,
h.surface_refractivity, h.ducting_detected,
ps.lat, ps.lon,
EXTRACT(HOUR FROM ps.valid_time) AS utc_hour,
EXTRACT(MONTH FROM ps.valid_time) AS month,
ps.lon
si.sfi, COALESCE((SELECT MAX(v) FROM unnest(si.kp_values) AS v), 2.0) AS kp_max,
snd.k_index, snd.lifted_index
FROM propagation_scores ps
JOIN hrrr_profiles h
ON h.lat = ps.lat AND h.lon = ps.lon AND h.valid_time = ps.valid_time
LEFT JOIN solar_indices si ON si.date = ps.valid_time::date
LEFT JOIN LATERAL (
SELECT s.k_index, s.lifted_index
FROM soundings s
WHERE s.observed_at BETWEEN ps.valid_time - interval '6 hours'
AND ps.valid_time + interval '6 hours'
ORDER BY ABS(EXTRACT(EPOCH FROM s.observed_at - ps.valid_time))
LIMIT 1
) snd ON true
WHERE ps.band_mhz = #{band_mhz}
AND h.surface_temp_c IS NOT NULL
AND h.surface_dewpoint_c IS NOT NULL
@ -167,48 +181,60 @@ defmodule Mix.Tasks.PropagationTrain do
{feature_rows, target_rows} =
rows
|> Enum.map(fn [
score,
band_mhz,
temp_c,
dewpoint_c,
pressure_mb,
grad,
hpbl,
pwat,
refractivity,
lat,
utc_hour,
month,
lon
] ->
features =
Model.encode_features(%{
surface_temp_c: to_float(temp_c),
surface_dewpoint_c: to_float(dewpoint_c),
surface_pressure_mb: to_float(pressure_mb),
min_refractivity_gradient: to_float(grad),
hpbl_m: to_float(hpbl),
pwat_mm: to_float(pwat),
surface_refractivity: to_float(refractivity),
latitude: to_float(lat),
utc_hour: to_float(utc_hour),
month: trunc(to_float(month)),
longitude: to_float(lon),
freq_mhz: band_mhz
})
{features, [score / 100.0]}
end)
|> Enum.map(&encode_pretrain_row/1)
|> Enum.unzip()
{Nx.tensor(feature_rows, type: :f32), Nx.tensor(target_rows, type: :f32), band_counts}
end
# ── Phase 2 data: real QSO-HRRR matches ──────────────────────────
defp encode_pretrain_row([
score,
band_mhz,
temp_c,
dewpoint_c,
pressure_mb,
grad,
hpbl,
pwat,
refractivity,
ducting,
lat,
lon,
utc_hour,
month,
sfi,
kp_max,
k_index,
lifted_index
]) do
features =
Model.encode_features(%{
surface_temp_c: to_float(temp_c),
surface_dewpoint_c: to_float(dewpoint_c),
surface_pressure_mb: to_float(pressure_mb),
min_refractivity_gradient: to_float(grad),
hpbl_m: to_float(hpbl),
pwat_mm: to_float(pwat),
surface_refractivity: to_float(refractivity),
latitude: to_float(lat),
sfi: to_float(sfi),
kp_max: to_float(kp_max),
ducting_detected: ducting,
k_index: to_float(k_index),
lifted_index: to_float(lifted_index),
utc_hour: to_float(utc_hour),
month: trunc(to_float(month)),
longitude: to_float(lon),
freq_mhz: band_mhz
})
{features, [score / 100.0]}
end
# ── Phase 2 data: QSO-HRRR + solar + sounding ────────────────────
defp load_qso_data do
IO.puts(" Running QSO-HRRR join query...")
IO.puts(" Running QSO-HRRR-solar-sounding join query...")
sql = """
SELECT
@ -231,7 +257,12 @@ defmodule Mix.Tasks.PropagationTrain do
(COALESCE(h1.pwat_mm, h2.pwat_mm) +
COALESCE(h2.pwat_mm, h1.pwat_mm)) / 2.0 AS avg_pwat,
(COALESCE(h1.surface_refractivity, h2.surface_refractivity) +
COALESCE(h2.surface_refractivity, h1.surface_refractivity)) / 2.0 AS avg_refractivity
COALESCE(h2.surface_refractivity, h1.surface_refractivity)) / 2.0 AS avg_refractivity,
COALESCE(h1.ducting_detected, false) OR COALESCE(h2.ducting_detected, false) AS ducting_either,
si.sfi,
COALESCE((SELECT MAX(v) FROM unnest(si.kp_values) AS v), 2.0) AS kp_max,
snd.k_index,
snd.lifted_index
FROM qsos q
INNER JOIN hrrr_profiles h1
ON h1.lat = ROUND((q.pos1->>'lat')::numeric * 8) / 8
@ -241,6 +272,15 @@ defmodule Mix.Tasks.PropagationTrain do
ON h2.lat = ROUND((q.pos2->>'lat')::numeric * 8) / 8
AND h2.lon = ROUND((q.pos2->>'lng')::numeric * 8) / 8
AND h2.valid_time = date_trunc('hour', q.qso_timestamp)
LEFT JOIN solar_indices si ON si.date = q.qso_timestamp::date
LEFT JOIN LATERAL (
SELECT s.k_index, s.lifted_index
FROM soundings s
WHERE s.observed_at BETWEEN q.qso_timestamp - interval '6 hours'
AND q.qso_timestamp + interval '6 hours'
ORDER BY ABS(EXTRACT(EPOCH FROM s.observed_at - q.qso_timestamp))
LIMIT 1
) snd ON true
WHERE q.distance_km > 0
AND q.distance_km < 3000
AND q.qso_timestamp >= '2016-06-30'
@ -249,7 +289,7 @@ defmodule Mix.Tasks.PropagationTrain do
"""
%{rows: rows} = Repo.query!(sql, [], timeout: @query_timeout)
IO.puts(" Matched #{length(rows)} QSOs to HRRR conditions")
IO.puts(" Matched #{length(rows)} QSOs to HRRR+solar+sounding")
band_counts = Enum.frequencies_by(rows, fn row -> Enum.at(row, 0) end)
@ -272,39 +312,7 @@ defmodule Mix.Tasks.PropagationTrain do
{feature_rows, target_rows} =
percentile_rows
|> Enum.map(fn {[
band_mhz,
_dist,
utc_hour,
month,
lon,
lat,
temp,
dewpoint,
pressure,
grad,
hpbl,
pwat,
refractivity
], percentile} ->
features =
Model.encode_features(%{
surface_temp_c: to_float(temp),
surface_dewpoint_c: to_float(dewpoint),
surface_pressure_mb: to_float(pressure),
min_refractivity_gradient: to_float(grad),
hpbl_m: to_float(hpbl),
pwat_mm: to_float(pwat),
surface_refractivity: to_float(refractivity),
latitude: to_float(lat),
utc_hour: to_float(utc_hour),
month: trunc(to_float(month)),
longitude: to_float(lon),
freq_mhz: band_mhz
})
{features, [percentile]}
end)
|> Enum.map(fn {row, percentile} -> encode_qso_row(row, percentile) end)
|> Enum.unzip()
IO.puts(" After filtering: #{length(feature_rows)} training samples")
@ -312,6 +320,53 @@ defmodule Mix.Tasks.PropagationTrain do
{Nx.tensor(feature_rows, type: :f32), Nx.tensor(target_rows, type: :f32), band_counts}
end
defp encode_qso_row(
[
band_mhz,
_dist,
utc_hour,
month,
lon,
lat,
temp,
dewpoint,
pressure,
grad,
hpbl,
pwat,
refractivity,
ducting,
sfi,
kp_max,
k_index,
lifted_index
],
percentile
) do
features =
Model.encode_features(%{
surface_temp_c: to_float(temp),
surface_dewpoint_c: to_float(dewpoint),
surface_pressure_mb: to_float(pressure),
min_refractivity_gradient: to_float(grad),
hpbl_m: to_float(hpbl),
pwat_mm: to_float(pwat),
surface_refractivity: to_float(refractivity),
latitude: to_float(lat),
sfi: to_float(sfi),
kp_max: to_float(kp_max),
ducting_detected: ducting,
k_index: to_float(k_index),
lifted_index: to_float(lifted_index),
utc_hour: to_float(utc_hour),
month: trunc(to_float(month)),
longitude: to_float(lon),
freq_mhz: band_mhz
})
{features, [percentile]}
end
# ── Helpers ───────────────────────────────────────────────────────
defp shuffle(features, targets, n) do

24
test/microwaveprop/propagation/model_test.exs
View file

@ -33,7 +33,7 @@ defmodule Microwaveprop.Propagation.ModelTest do
assert {:ok, loaded} = Model.load(path)
# Predict with both and compare
input = Nx.broadcast(0.5, {1, 15})
input = Nx.broadcast(0.5, {1, 20})
original = Model.predict(params, input)
reloaded = Model.predict(loaded, input)
@ -53,7 +53,7 @@ defmodule Microwaveprop.Propagation.ModelTest do
Model.save(params, path)
loaded = Model.load_or_init(path)
input = Nx.broadcast(0.5, {1, 15})
input = Nx.broadcast(0.5, {1, 20})
assert Nx.to_flat_list(Model.predict(params, input)) == Nx.to_flat_list(Model.predict(loaded, input))
end
@ -66,8 +66,8 @@ defmodule Microwaveprop.Propagation.ModelTest do
describe "predict/2" do
test "produces output in [0, 1] range" do
params = Model.init()
# Single sample, 15 features
input = Nx.broadcast(0.5, {1, 15})
# Single sample, 20 features
input = Nx.broadcast(0.5, {1, 20})
output = Model.predict(params, input)
assert Nx.shape(output) == {1, 1}
@ -78,7 +78,7 @@ defmodule Microwaveprop.Propagation.ModelTest do
test "handles batch of multiple samples" do
params = Model.init()
input = Nx.broadcast(0.5, {4, 15})
input = Nx.broadcast(0.5, {4, 20})
output = Model.predict(params, input)
assert Nx.shape(output) == {4, 1}
@ -86,7 +86,7 @@ defmodule Microwaveprop.Propagation.ModelTest do
end
describe "encode_features/1" do
test "returns 15-element list" do
test "returns 20-element list" do
features =
Model.encode_features(%{
surface_temp_c: 25.0,
@ -100,13 +100,13 @@ defmodule Microwaveprop.Propagation.ModelTest do
freq_mhz: 10_000
})
assert length(features) == 15
assert length(features) == 20
assert Enum.all?(features, &is_float/1)
end
test "uses defaults for missing values" do
features = Model.encode_features(%{})
assert length(features) == 15
assert length(features) == 20
assert Enum.all?(features, &is_float/1)
end
@ -129,7 +129,7 @@ defmodule Microwaveprop.Propagation.ModelTest do
test "trains on synthetic data and returns state with metrics" do
# Generate small synthetic dataset: 64 samples, 13 features
key = Nx.Random.key(42)
{features, key} = Nx.Random.uniform(key, shape: {64, 15}, type: :f32)
{features, key} = Nx.Random.uniform(key, shape: {64, 20}, type: :f32)
{targets, _key} = Nx.Random.uniform(key, shape: {64, 1}, type: :f32)
{trained_state, metrics} = Model.train(features, targets, epochs: 2, batch_size: 32)
@ -144,7 +144,7 @@ defmodule Microwaveprop.Propagation.ModelTest do
test "returns rmse and r_squared" do
params = Model.init()
key = Nx.Random.key(42)
{features, key} = Nx.Random.uniform(key, shape: {32, 15}, type: :f32)
{features, key} = Nx.Random.uniform(key, shape: {32, 20}, type: :f32)
{targets, _key} = Nx.Random.uniform(key, shape: {32, 1}, type: :f32)
metrics = Model.evaluate(params, features, targets)
@ -213,8 +213,8 @@ defmodule Microwaveprop.Propagation.ModelTest do
end
describe "feature_names/0" do
test "returns 15 feature names" do
assert length(Model.feature_names()) == 15
test "returns 20 feature names" do
assert length(Model.feature_names()) == 20
end
end
end