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(", ")}") +```