diff --git a/lib/microwaveprop/propagation/model.ex b/lib/microwaveprop/propagation/model.ex index 2c505c6f..326fa9bb 100644 --- a/lib/microwaveprop/propagation/model.ex +++ b/lib/microwaveprop/propagation/model.ex @@ -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), diff --git a/lib/mix/tasks/propagation_train.ex b/lib/mix/tasks/propagation_train.ex index 31abd829..e33ab2a0 100644 --- a/lib/mix/tasks/propagation_train.ex +++ b/lib/mix/tasks/propagation_train.ex @@ -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 diff --git a/test/microwaveprop/propagation/model_test.exs b/test/microwaveprop/propagation/model_test.exs index a90a474c..c3426b7b 100644 --- a/test/microwaveprop/propagation/model_test.exs +++ b/test/microwaveprop/propagation/model_test.exs @@ -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