From 7f0e1f41827c7e1678447b7125b23b3f9409fb9c Mon Sep 17 00:00:00 2001 From: Graham McIntire Date: Tue, 7 Apr 2026 08:45:10 -0500 Subject: [PATCH] Add Livebook notebook for propagation algorithm analysis Connects to local dev DB to analyze contacts, HRRR conditions, terrain impact, commercial link correlation, and score sensitivity. Ready for use once prod data is imported locally. --- notebooks/algorithm_analysis.livemd | 410 ++++++++++++++++++++++++++++ 1 file changed, 410 insertions(+) create mode 100644 notebooks/algorithm_analysis.livemd diff --git a/notebooks/algorithm_analysis.livemd b/notebooks/algorithm_analysis.livemd new file mode 100644 index 00000000..19927192 --- /dev/null +++ b/notebooks/algorithm_analysis.livemd @@ -0,0 +1,410 @@ +# Propagation Algorithm Analysis + +```elixir +Mix.install([ + {:postgrex, "~> 0.19"}, + {:ecto_sql, "~> 3.12"}, + {:kino, "~> 0.14"}, + {:kino_vega_lite, "~> 0.1"}, + {:explorer, "~> 0.10"}, + {:statistex, "~> 1.0"} +]) +``` + +## Database Connection + +```elixir +{:ok, conn} = Postgrex.start_link( + hostname: "localhost", + database: "microwaveprop_dev", + username: "postgres", + password: "postgres", + port: 5432 +) +``` + +## Helper Functions + +```elixir +defmodule Q do + @doc "Run a query and return results as a list of maps" + def query!(conn, sql, params \\ []) do + %{columns: cols, rows: rows} = Postgrex.query!(conn, sql, params) + Enum.map(rows, fn row -> Enum.zip(cols, row) |> Map.new() end) + end + + @doc "Run a query and return as an Explorer DataFrame" + def df!(conn, sql, params \\ []) do + %{columns: cols, rows: rows} = Postgrex.query!(conn, sql, params) + cols + |> Enum.with_index() + |> Enum.map(fn {col, i} -> {col, Enum.map(rows, &Enum.at(&1, i))} end) + |> Map.new() + |> Explorer.DataFrame.new() + end +end +``` + +## 1. Contact Overview + +```elixir +Q.query!(conn, """ + SELECT COUNT(*) as total, + COUNT(*) FILTER (WHERE distance_km < 3000) as tropo, + COUNT(*) FILTER (WHERE distance_km >= 3000) as eme, + COUNT(DISTINCT band) as bands, + MIN(qso_timestamp) as earliest, + MAX(qso_timestamp) as latest + FROM contacts +""") +|> hd() +|> Kino.Tree.new() +``` + +```elixir +# Contacts by band +Q.df!(conn, """ + SELECT band::integer as band_mhz, COUNT(*) as count, + ROUND(AVG(distance_km::numeric), 1) as avg_dist_km, + ROUND(MAX(distance_km::numeric), 1) as max_dist_km + FROM contacts + WHERE distance_km < 3000 + GROUP BY band ORDER BY band +""") +|> Kino.DataTable.new() +``` + +```elixir +# Contacts by month — when do contacts happen? +contacts_by_month = Q.df!(conn, """ + SELECT EXTRACT(MONTH FROM qso_timestamp)::integer as month, + COUNT(*) as count + FROM contacts WHERE distance_km < 3000 + GROUP BY 1 ORDER BY 1 +""") + +VegaLite.new(width: 600, height: 300, title: "Tropo Contacts by Month") +|> VegaLite.data_from_values(contacts_by_month |> Explorer.DataFrame.to_rows()) +|> VegaLite.mark(:bar) +|> VegaLite.encode_field(:x, "month", type: :ordinal) +|> VegaLite.encode_field(:y, "count", type: :quantitative) +``` + +```elixir +# Contacts by hour of day +contacts_by_hour = Q.df!(conn, """ + SELECT EXTRACT(HOUR FROM qso_timestamp)::integer as utc_hour, + COUNT(*) as count + FROM contacts WHERE distance_km < 3000 + GROUP BY 1 ORDER BY 1 +""") + +VegaLite.new(width: 600, height: 300, title: "Tropo Contacts by UTC Hour") +|> VegaLite.data_from_values(contacts_by_hour |> Explorer.DataFrame.to_rows()) +|> VegaLite.mark(:bar) +|> VegaLite.encode_field(:x, "utc_hour", type: :ordinal) +|> VegaLite.encode_field(:y, "count", type: :quantitative) +``` + +## 2. HRRR Conditions at Contact Time + +Join contacts with their nearest HRRR profile to see what atmospheric +conditions were present when contacts actually happened. + +```elixir +# Contacts with HRRR data — the core analysis dataset +contact_hrrr = Q.df!(conn, """ + SELECT + c.id, + c.band::integer as band_mhz, + c.distance_km::float as distance_km, + c.qso_timestamp, + EXTRACT(MONTH FROM c.qso_timestamp)::integer as month, + EXTRACT(HOUR FROM c.qso_timestamp)::integer as utc_hour, + h.surface_temp as temp_c, + h.surface_dewpoint as dewpoint_c, + (h.surface_temp - h.surface_dewpoint) as td_depression_c, + h.surface_pressure as pressure_mb, + h.surface_refractivity, + h.min_refractivity_gradient, + h.hpbl_m as bl_depth_m, + h.pwat_mm, + h.ducting_detected + FROM contacts c + JOIN LATERAL ( + SELECT * + FROM hrrr_profiles hp + WHERE hp.lat BETWEEN ROUND((c.pos1->>'lat')::numeric, 2) - 0.07 + AND ROUND((c.pos1->>'lat')::numeric, 2) + 0.07 + AND hp.lon BETWEEN ROUND(COALESCE((c.pos1->>'lon')::numeric, (c.pos1->>'lng')::numeric), 2) - 0.07 + AND ROUND(COALESCE((c.pos1->>'lon')::numeric, (c.pos1->>'lng')::numeric), 2) + 0.07 + AND hp.valid_time BETWEEN c.qso_timestamp - interval '1 hour' + AND c.qso_timestamp + interval '1 hour' + ORDER BY ABS(EXTRACT(EPOCH FROM hp.valid_time - c.qso_timestamp)) + LIMIT 1 + ) h ON true + WHERE c.distance_km < 3000 + AND c.pos1 IS NOT NULL +""") + +Kino.DataTable.new(contact_hrrr, name: "Contacts + HRRR") +``` + +## 3. Factor Distributions — What Does the Atmosphere Look Like During Contacts? + +```elixir +rows = Explorer.DataFrame.to_rows(contact_hrrr) + +# Refractivity gradient distribution +VegaLite.new(width: 600, height: 300, title: "Refractivity Gradient at Contact Time") +|> VegaLite.data_from_values(rows) +|> VegaLite.mark(:bar) +|> VegaLite.encode_field(:x, "min_refractivity_gradient", + type: :quantitative, bin: %{step: 20}, title: "Min Refractivity Gradient (N-units/km)" +) +|> VegaLite.encode(:y, aggregate: :count, type: :quantitative) +``` + +```elixir +# Td depression distribution +VegaLite.new(width: 600, height: 300, title: "Td Depression at Contact Time") +|> VegaLite.data_from_values(rows) +|> VegaLite.mark(:bar) +|> VegaLite.encode_field(:x, "td_depression_c", + type: :quantitative, bin: %{step: 1}, title: "T - Td (C)" +) +|> VegaLite.encode(:y, aggregate: :count, type: :quantitative) +``` + +```elixir +# Boundary layer depth +VegaLite.new(width: 600, height: 300, title: "Boundary Layer Depth at Contact Time") +|> VegaLite.data_from_values(rows) +|> VegaLite.mark(:bar) +|> VegaLite.encode_field(:x, "bl_depth_m", + type: :quantitative, bin: %{step: 100}, title: "BL Depth (m)" +) +|> VegaLite.encode(:y, aggregate: :count, type: :quantitative) +``` + +```elixir +# PWAT distribution +VegaLite.new(width: 600, height: 300, title: "Precipitable Water at Contact Time") +|> VegaLite.data_from_values(rows) +|> VegaLite.mark(:bar) +|> VegaLite.encode_field(:x, "pwat_mm", + type: :quantitative, bin: %{step: 3}, title: "PWAT (mm)" +) +|> VegaLite.encode(:y, aggregate: :count, type: :quantitative) +``` + +## 4. Ducting Analysis + +```elixir +# Ducting rate by month +ducting_monthly = Q.df!(conn, """ + SELECT + EXTRACT(MONTH FROM c.qso_timestamp)::integer as month, + COUNT(*) as total, + COUNT(*) FILTER (WHERE h.ducting_detected = true) as ducting, + ROUND(100.0 * COUNT(*) FILTER (WHERE h.ducting_detected = true) / COUNT(*), 1) as ducting_pct + FROM contacts c + JOIN LATERAL ( + SELECT ducting_detected + FROM hrrr_profiles hp + WHERE hp.lat BETWEEN ROUND((c.pos1->>'lat')::numeric, 2) - 0.07 + AND ROUND((c.pos1->>'lat')::numeric, 2) + 0.07 + AND hp.lon BETWEEN ROUND(COALESCE((c.pos1->>'lon')::numeric, (c.pos1->>'lng')::numeric), 2) - 0.07 + AND ROUND(COALESCE((c.pos1->>'lon')::numeric, (c.pos1->>'lng')::numeric), 2) + 0.07 + AND hp.valid_time BETWEEN c.qso_timestamp - interval '1 hour' + AND c.qso_timestamp + interval '1 hour' + ORDER BY ABS(EXTRACT(EPOCH FROM hp.valid_time - c.qso_timestamp)) + LIMIT 1 + ) h ON true + WHERE c.distance_km < 3000 AND c.pos1 IS NOT NULL + GROUP BY 1 ORDER BY 1 +""") + +VegaLite.new(width: 600, height: 300, title: "Ducting Detection Rate by Month (at contact time)") +|> VegaLite.data_from_values(ducting_monthly |> Explorer.DataFrame.to_rows()) +|> VegaLite.mark(:bar) +|> VegaLite.encode_field(:x, "month", type: :ordinal) +|> VegaLite.encode_field(:y, "ducting_pct", type: :quantitative, title: "% Ducting") +``` + +## 5. Distance vs Atmospheric Conditions + +Do contacts at longer distances correlate with better atmospheric conditions? + +```elixir +# Distance vs refractivity gradient (scatter) +VegaLite.new(width: 600, height: 400, title: "Distance vs Refractivity Gradient") +|> VegaLite.data_from_values(rows) +|> VegaLite.mark(:circle, opacity: 0.3, size: 15) +|> VegaLite.encode_field(:x, "min_refractivity_gradient", + type: :quantitative, title: "Min Refractivity Gradient" +) +|> VegaLite.encode_field(:y, "distance_km", + type: :quantitative, title: "Distance (km)" +) +|> VegaLite.encode_field(:color, "band_mhz", type: :nominal, title: "Band (MHz)") +``` + +```elixir +# Distance vs Td depression +VegaLite.new(width: 600, height: 400, title: "Distance vs Td Depression") +|> VegaLite.data_from_values(rows) +|> VegaLite.mark(:circle, opacity: 0.3, size: 15) +|> VegaLite.encode_field(:x, "td_depression_c", + type: :quantitative, title: "T - Td (C)" +) +|> VegaLite.encode_field(:y, "distance_km", + type: :quantitative, title: "Distance (km)" +) +|> VegaLite.encode_field(:color, "band_mhz", type: :nominal, title: "Band (MHz)") +``` + +## 6. Terrain Impact + +```elixir +terrain_stats = Q.df!(conn, """ + SELECT + tp.verdict, + COUNT(*) as count, + ROUND(AVG(c.distance_km::numeric), 1) as avg_dist_km, + ROUND(AVG(tp.diffraction_db::numeric), 1) as avg_diffraction_db, + ROUND(AVG(tp.min_clearance_m::numeric), 1) as avg_clearance_m + FROM terrain_profiles tp + JOIN contacts c ON c.id = tp.contact_id + WHERE c.distance_km < 3000 + GROUP BY tp.verdict + ORDER BY count DESC +""") + +Kino.DataTable.new(terrain_stats, name: "Terrain Verdicts") +``` + +```elixir +# Diffraction loss vs distance +terrain_scatter = Q.df!(conn, """ + SELECT + c.band::integer as band_mhz, + c.distance_km::float as distance_km, + tp.diffraction_db::float as diffraction_db, + tp.verdict + FROM terrain_profiles tp + JOIN contacts c ON c.id = tp.contact_id + WHERE c.distance_km < 3000 +""") + +VegaLite.new(width: 600, height: 400, title: "Diffraction Loss vs Distance") +|> VegaLite.data_from_values(terrain_scatter |> Explorer.DataFrame.to_rows()) +|> VegaLite.mark(:circle, opacity: 0.2, size: 10) +|> VegaLite.encode_field(:x, "distance_km", type: :quantitative, title: "Distance (km)") +|> VegaLite.encode_field(:y, "diffraction_db", type: :quantitative, title: "Diffraction Loss (dB)") +|> VegaLite.encode_field(:color, "verdict", type: :nominal) +``` + +## 7. Commercial Link Correlation + +Ground truth: do commercial link signal levels correlate with atmospheric conditions? + +```elixir +commercial = Q.df!(conn, """ + SELECT + cl.name, cl.frequency_ghz, + cs.sampled_at, + EXTRACT(HOUR FROM cs.sampled_at)::integer as utc_hour, + EXTRACT(MONTH FROM cs.sampled_at)::integer as month, + cs.rx_power_0::float as rx_power_dbm + FROM commercial_samples cs + JOIN commercial_links cl ON cl.id = cs.link_id + WHERE cs.rx_power_0 IS NOT NULL + ORDER BY cs.sampled_at +""") + +Kino.DataTable.new(commercial, name: "Commercial Samples") +``` + +```elixir +# Signal strength by hour of day per link +VegaLite.new(width: 600, height: 300, title: "Commercial Link Signal by Hour of Day") +|> VegaLite.data_from_values(commercial |> Explorer.DataFrame.to_rows()) +|> VegaLite.mark(:line, point: true) +|> VegaLite.encode_field(:x, "utc_hour", type: :ordinal) +|> VegaLite.encode_field(:y, "rx_power_dbm", + type: :quantitative, aggregate: :mean, title: "Mean RX Power (dBm)" +) +|> VegaLite.encode_field(:color, "name", type: :nominal) +``` + +## 8. Scoring the Historical Contacts + +Recompute scores for all contacts using the current algorithm and compare +factor weights. This section is a template — once you import the full dump +and attach to the running node, you can call `Scorer.composite_score/2` directly. + +```elixir +# For now, compute a simplified composite score from HRRR data +# TODO: attach to running node for full Scorer access + +# Humidity score approximation (beneficial band) +humidity_score = fn temp_c, dewpoint_c -> + if is_nil(temp_c) or is_nil(dewpoint_c), do: 50, + else: min(100, max(0, round(217 * 6.112 * :math.exp(17.67 * dewpoint_c / (dewpoint_c + 243.5)) / (temp_c + 273.15) * 5))) +end + +# Td depression score (beneficial band) +td_score = fn temp_c, dewpoint_c -> + case temp_c - dewpoint_c do + d when d <= 1 -> 95 + d when d <= 2 -> 88 + d when d <= 3 -> 80 + d when d <= 5 -> 65 + d when d <= 8 -> 50 + d when d <= 12 -> 35 + _ -> 20 + end +end + +# Refractivity score +refrac_score = fn grad -> + cond do + is_nil(grad) -> 50 + grad < -200 -> 98 + grad < -150 -> 92 + grad < -100 -> 82 + grad < -75 -> 68 + grad < -55 -> 55 + grad < -40 -> 48 + true -> 42 + end +end + +IO.puts("Scoring functions defined. Use with Explorer.DataFrame.mutate/2 or Enum.map/2") +``` + +## 9. Weight Sensitivity Analysis + +Template for testing different weight configurations against historical data. + +```elixir +# Define weight sets to test +weight_sets = %{ + current: %{humidity: 0.18, time_of_day: 0.10, td_depression: 0.10, + refractivity: 0.08, sky: 0.08, season: 0.08, + wind: 0.05, rain: 0.08, pressure: 0.15, pwat: 0.10}, + + refractivity_heavy: %{humidity: 0.15, time_of_day: 0.08, td_depression: 0.10, + refractivity: 0.18, sky: 0.06, season: 0.08, + wind: 0.04, rain: 0.06, pressure: 0.12, pwat: 0.13}, + + humidity_heavy: %{humidity: 0.25, time_of_day: 0.08, td_depression: 0.12, + refractivity: 0.08, sky: 0.06, season: 0.08, + wind: 0.04, rain: 0.06, pressure: 0.12, pwat: 0.11} +} + +# TODO: score all contacts with each weight set, +# then correlate scores with distance achieved and commercial link signal levels +IO.puts("Weight sets defined: #{Map.keys(weight_sets) |> Enum.join(", ")}") +```