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.
This commit is contained in:
parent
6b9ef24c81
commit
af5a48b797
1 changed files with 410 additions and 0 deletions
410
notebooks/algorithm_analysis.livemd
Normal file
410
notebooks/algorithm_analysis.livemd
Normal file
|
|
@ -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(", ")}")
|
||||
```
|
||||
Loading…
Add table
Reference in a new issue