prop/notebooks/algorithm_analysis.livemd

Propagation Algorithm Analysis

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

{:ok, conn} = Postgrex.start_link(
  hostname: "localhost",
  database: "microwaveprop_dev",
  username: "postgres",
  password: "postgres",
  port: 5432
)

Helper Functions

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

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()

# 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()

# 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)

# 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.

# 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?

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)

# 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)

# 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)

# 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

# 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?

# 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)")

# 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

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")

# 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?

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")

# 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.

# 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.

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