From 4bf9d0b4b0e0cbfeca0e5d5a552655ac67c4a876 Mon Sep 17 00:00:00 2001 From: Graham McIntire Date: Wed, 1 Apr 2026 08:34:04 -0500 Subject: [PATCH] Add mix propagation.analyze task for QSO-HRRR correlation analysis --- lib/mix/tasks/propagation_analyze.ex | 612 +++++++++++++++++++++++++++ 1 file changed, 612 insertions(+) create mode 100644 lib/mix/tasks/propagation_analyze.ex diff --git a/lib/mix/tasks/propagation_analyze.ex b/lib/mix/tasks/propagation_analyze.ex new file mode 100644 index 00000000..1744a542 --- /dev/null +++ b/lib/mix/tasks/propagation_analyze.ex @@ -0,0 +1,612 @@ +defmodule Mix.Tasks.PropagationAnalyze do + @shortdoc "Analyze QSO-HRRR correlations and factor effects on propagation distance" + + @moduledoc """ + Builds a dataset matching QSOs to HRRR atmospheric conditions at both endpoints, + then runs statistical analysis to identify which factors drive propagation distance. + + Analyses performed: + 1. Spearman rank correlation per band (10G, 24G, 47G, 75G) + 2. Binned factor analysis with median/p25/p75 distances (10G, 24G) + 3. Interaction effects for 10 GHz (gradient×time, HPBL×season, ducting) + + mix propagation.analyze + """ + + use Mix.Task + + alias Microwaveprop.Repo + + @bands [10_000, 24_000, 47_000, 75_000] + @timeout 300_000 + + @impl Mix.Task + def run(_args) do + Application.put_env( + :microwaveprop, + Oban, + Keyword.put( + Application.get_env(:microwaveprop, Oban, []), + :queues, + false + ) + ) + + Mix.Task.run("app.start") + + "=" |> String.duplicate(80) |> IO.puts() + IO.puts("PROPAGATION ANALYSIS: QSO-HRRR Correlation Study") + "=" |> String.duplicate(80) |> IO.puts() + IO.puts("Started: #{DateTime.to_string(DateTime.utc_now())}") + IO.puts("") + + dataset = build_dataset() + IO.puts("Dataset: #{length(dataset)} matched QSO-HRRR rows\n") + + print_dataset_summary(dataset) + run_correlation_analysis(dataset) + run_factor_analysis(dataset) + run_interaction_analysis(dataset) + + IO.puts("\n" <> String.duplicate("=", 80)) + IO.puts("Analysis complete: #{DateTime.to_string(DateTime.utc_now())}") + end + + # --------------------------------------------------------------------------- + # Dataset builder + # --------------------------------------------------------------------------- + + defp build_dataset do + IO.puts("Building dataset (joining QSOs × HRRR at both endpoints)...") + + sql = """ + SELECT + q.id AS qso_id, + q.band::float AS band, + q.distance_km::float AS distance_km, + q.qso_timestamp, + EXTRACT(MONTH FROM q.qso_timestamp)::int AS month, + EXTRACT(HOUR FROM q.qso_timestamp)::int AS utc_hour, + + -- Endpoint 1 atmospheric + h1.surface_temp_c AS h1_temp, + h1.surface_dewpoint_c AS h1_dewpoint, + h1.surface_pressure_mb AS h1_pressure, + h1.surface_refractivity AS h1_refractivity, + h1.min_refractivity_gradient AS h1_gradient, + h1.hpbl_m AS h1_hpbl, + h1.pwat_mm AS h1_pwat, + h1.ducting_detected AS h1_ducting, + + -- Endpoint 2 atmospheric + h2.surface_temp_c AS h2_temp, + h2.surface_dewpoint_c AS h2_dewpoint, + h2.surface_pressure_mb AS h2_pressure, + h2.surface_refractivity AS h2_refractivity, + h2.min_refractivity_gradient AS h2_gradient, + h2.hpbl_m AS h2_hpbl, + h2.pwat_mm AS h2_pwat, + h2.ducting_detected AS h2_ducting, + + -- Terrain + tp.verdict AS terrain_verdict + + FROM qsos q + + INNER JOIN hrrr_profiles h1 + ON h1.lat = ROUND((q.pos1->>'lat')::numeric * 8) / 8 + AND h1.lon = ROUND((q.pos1->>'lng')::numeric * 8) / 8 + AND h1.valid_time = date_trunc('hour', q.qso_timestamp) + + INNER JOIN hrrr_profiles h2 + 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 terrain_profiles tp ON tp.qso_id = q.id + + WHERE q.distance_km > 0 + AND q.distance_km < 3000 + AND q.qso_timestamp >= '2016-06-30' + AND q.pos1 IS NOT NULL + AND q.pos2 IS NOT NULL + """ + + %{rows: rows, columns: columns} = Repo.query!(sql, [], timeout: @timeout) + + Enum.map(rows, fn row -> + columns + |> Enum.zip(row) + |> Map.new(fn {col, val} -> {String.to_atom(col), val} end) + |> derive_averages() + end) + end + + defp derive_averages(row) do + row + |> Map.put(:avg_temp, safe_avg(row[:h1_temp], row[:h2_temp])) + |> Map.put(:avg_dewpoint, safe_avg(row[:h1_dewpoint], row[:h2_dewpoint])) + |> Map.put(:avg_pressure, safe_avg(row[:h1_pressure], row[:h2_pressure])) + |> Map.put(:avg_refractivity, safe_avg(row[:h1_refractivity], row[:h2_refractivity])) + |> Map.put(:avg_gradient, safe_avg(row[:h1_gradient], row[:h2_gradient])) + |> Map.put(:avg_hpbl, safe_avg(row[:h1_hpbl], row[:h2_hpbl])) + |> Map.put(:avg_pwat, safe_avg(row[:h1_pwat], row[:h2_pwat])) + |> Map.put(:either_ducting, row[:h1_ducting] == true or row[:h2_ducting] == true) + end + + defp safe_avg(nil, nil), do: nil + defp safe_avg(nil, b), do: b + defp safe_avg(a, nil), do: a + defp safe_avg(a, b), do: (a + b) / 2.0 + + # --------------------------------------------------------------------------- + # Dataset summary + # --------------------------------------------------------------------------- + + defp print_dataset_summary(dataset) do + IO.puts(String.duplicate("-", 80)) + IO.puts("DATASET SUMMARY") + IO.puts(String.duplicate("-", 80)) + + by_band = + dataset + |> Enum.group_by(& &1.band) + |> Enum.sort_by(fn {band, _} -> band end) + + IO.puts("\nQSOs per band:") + + for {band, rows} <- by_band do + distances = Enum.map(rows, & &1.distance_km) + + IO.puts( + " #{format_band(band)}: n=#{length(rows)}, median=#{distances |> median() |> format_num(1)} km, " <> + "p25=#{distances |> percentile(25) |> format_num(1)} km, p75=#{distances |> percentile(75) |> format_num(1)} km" + ) + end + + terrain_counts = + dataset + |> Enum.frequencies_by(& &1.terrain_verdict) + |> Enum.sort_by(fn {_, count} -> -count end) + + IO.puts("\nTerrain verdicts:") + + for {verdict, count} <- terrain_counts do + IO.puts(" #{verdict || "NULL"}: #{count}") + end + + IO.puts("") + end + + # --------------------------------------------------------------------------- + # Spearman rank correlation + # --------------------------------------------------------------------------- + + defp run_correlation_analysis(dataset) do + IO.puts(String.duplicate("-", 80)) + IO.puts("SPEARMAN RANK CORRELATION WITH DISTANCE (per band)") + IO.puts(String.duplicate("-", 80)) + + variables = [ + {:avg_temp, "Temperature (°C)"}, + {:avg_dewpoint, "Dewpoint (°C)"}, + {:avg_pressure, "Pressure (mb)"}, + {:avg_refractivity, "Surface Refractivity"}, + {:avg_gradient, "Refractivity Gradient"}, + {:avg_hpbl, "HPBL (m)"}, + {:avg_pwat, "PWAT (mm)"}, + {:month, "Month"}, + {:utc_hour, "UTC Hour"} + ] + + for band <- @bands do + band_data = Enum.filter(dataset, &(trunc(&1.band) == band)) + + IO.puts("\n#{format_band(band)} (n=#{length(band_data)}):") + + IO.puts( + " #{String.pad_trailing("Variable", 28)} #{String.pad_leading("rho", 8)} #{String.pad_leading("n_valid", 8)}" + ) + + IO.puts(" " <> String.duplicate("-", 46)) + + correlations = + variables + |> Enum.map(fn {key, label} -> + pairs = + band_data + |> Enum.filter(&(&1[key] != nil)) + |> Enum.map(&{to_float(&1[key]), to_float(&1.distance_km)}) + + rho = spearman(pairs) + {label, rho, length(pairs)} + end) + |> Enum.sort_by(fn {_, rho, _} -> -abs(rho || 0) end) + + for {label, rho, n} <- correlations do + rho_str = if rho, do: format_num(rho, 4), else: "N/A" + + IO.puts( + " #{String.pad_trailing(label, 28)} #{String.pad_leading(rho_str, 8)} #{String.pad_leading(Integer.to_string(n), 8)}" + ) + end + end + + IO.puts("") + end + + # --------------------------------------------------------------------------- + # Binned factor analysis + # --------------------------------------------------------------------------- + + defp run_factor_analysis(dataset) do + IO.puts(String.duplicate("-", 80)) + IO.puts("BINNED FACTOR ANALYSIS (median distance per bin)") + IO.puts(String.duplicate("-", 80)) + + factor_configs = [ + {:avg_gradient, "Refractivity Gradient", + [ + {nil, -500, "< -500 (strong)"}, + {-500, -300, "-500 to -300 (moderate)"}, + {-300, -200, "-300 to -200 (enhanced)"}, + {-200, -100, "-200 to -100 (mild)"}, + {-100, 0, "-100 to 0 (normal)"}, + {0, nil, "> 0"} + ]}, + {:avg_temp, "Temperature (°C)", + [ + {nil, 0, "< 0"}, + {0, 10, "0-10"}, + {10, 20, "10-20"}, + {20, 30, "20-30"}, + {30, nil, "> 30"} + ]}, + {:avg_hpbl, "HPBL (m)", + [ + {nil, 200, "< 200"}, + {200, 500, "200-500"}, + {500, 1000, "500-1000"}, + {1000, 2000, "1000-2000"}, + {2000, nil, "> 2000"} + ]}, + {:avg_pressure, "Pressure (mb)", + [ + {nil, 1005, "< 1005"}, + {1005, 1013, "1005-1013"}, + {1013, 1020, "1013-1020"}, + {1020, nil, "> 1020"} + ]}, + {:avg_pwat, "PWAT (mm)", + [ + {nil, 10, "< 10"}, + {10, 20, "10-20"}, + {20, 30, "20-30"}, + {30, 40, "30-40"}, + {40, nil, "> 40"} + ]}, + {:utc_hour, "UTC Hour", + [ + {0, 3, "00-02 UTC"}, + {3, 6, "03-05 UTC"}, + {6, 9, "06-08 UTC"}, + {9, 12, "09-11 UTC"}, + {12, 15, "12-14 UTC"}, + {15, 18, "15-17 UTC"}, + {18, 21, "18-20 UTC"}, + {21, 24, "21-23 UTC"} + ]}, + {:month, "Month", + [ + {1, 2, "Jan"}, + {2, 3, "Feb"}, + {3, 4, "Mar"}, + {4, 5, "Apr"}, + {5, 6, "May"}, + {6, 7, "Jun"}, + {7, 8, "Jul"}, + {8, 9, "Aug"}, + {9, 10, "Sep"}, + {10, 11, "Oct"}, + {11, 12, "Nov"}, + {12, 13, "Dec"} + ]} + ] + + terrain_config = + {:terrain_verdict, "Terrain Verdict", + [ + {:eq, "CLEAR", "CLEAR"}, + {:eq, "FRESNEL_PARTIAL", "FRESNEL_PARTIAL"}, + {:eq, "BLOCKED", "BLOCKED"} + ]} + + for band <- [10_000, 24_000] do + band_data = Enum.filter(dataset, &(trunc(&1.band) == band)) + + IO.puts("\n#{format_band(band)} (n=#{length(band_data)}):") + + for {key, label, bins} <- factor_configs do + print_binned_factor(band_data, key, label, bins) + end + + print_categorical_factor(band_data, terrain_config) + end + + IO.puts("") + end + + defp print_binned_factor(data, key, label, bins) do + IO.puts("\n #{label}:") + + IO.puts( + " #{String.pad_trailing("Bin", 28)} #{String.pad_leading("n", 6)} #{String.pad_leading("median", 8)} #{String.pad_leading("p25", 8)} #{String.pad_leading("p75", 8)}" + ) + + IO.puts(" " <> String.duplicate("-", 60)) + + for {low, high, bin_label} <- bins do + rows = + Enum.filter(data, fn row -> + val = row[key] + val != nil and in_bin?(val, low, high) + end) + + print_bin_row(bin_label, rows) + end + end + + defp print_categorical_factor(data, {key, label, bins}) do + IO.puts("\n #{label}:") + + IO.puts( + " #{String.pad_trailing("Bin", 28)} #{String.pad_leading("n", 6)} #{String.pad_leading("median", 8)} #{String.pad_leading("p25", 8)} #{String.pad_leading("p75", 8)}" + ) + + IO.puts(" " <> String.duplicate("-", 60)) + + for {:eq, match_val, bin_label} <- bins do + rows = Enum.filter(data, &(&1[key] == match_val)) + print_bin_row(bin_label, rows) + end + end + + defp print_bin_row(bin_label, rows) do + if rows == [] do + IO.puts( + " #{String.pad_trailing(bin_label, 28)} #{String.pad_leading("0", 6)} #{String.pad_leading("-", 8)} #{String.pad_leading("-", 8)} #{String.pad_leading("-", 8)}" + ) + else + distances = Enum.map(rows, & &1.distance_km) + + IO.puts( + " #{String.pad_trailing(bin_label, 28)} " <> + "#{String.pad_leading(Integer.to_string(length(rows)), 6)} " <> + "#{String.pad_leading(format_num(median(distances), 1), 8)} " <> + "#{String.pad_leading(format_num(percentile(distances, 25), 1), 8)} " <> + "#{String.pad_leading(format_num(percentile(distances, 75), 1), 8)}" + ) + end + end + + defp in_bin?(val, nil, high), do: val < high + defp in_bin?(val, low, nil), do: val >= low + defp in_bin?(val, low, high), do: val >= low and val < high + + # --------------------------------------------------------------------------- + # Interaction effects (10 GHz only) + # --------------------------------------------------------------------------- + + defp run_interaction_analysis(dataset) do + IO.puts(String.duplicate("-", 80)) + IO.puts("INTERACTION EFFECTS (10 GHz)") + IO.puts(String.duplicate("-", 80)) + + data_10g = Enum.filter(dataset, &(trunc(&1.band) == 10_000)) + + IO.puts("\nDataset: #{length(data_10g)} QSOs at 10 GHz\n") + + # Gradient × Time of day + IO.puts("1. Refractivity Gradient × Time of Day") + IO.puts(" (strong gradient = avg < -100, weak = avg >= -100)") + IO.puts("") + + time_blocks = [ + {0, 6, "00-05 UTC (night)"}, + {6, 12, "06-11 UTC (morning)"}, + {12, 18, "12-17 UTC (afternoon)"}, + {18, 24, "18-23 UTC (evening)"} + ] + + IO.puts( + " #{String.pad_trailing("Time Block", 24)} #{String.pad_trailing("Gradient", 10)} #{String.pad_leading("n", 6)} #{String.pad_leading("median km", 10)}" + ) + + IO.puts(" " <> String.duplicate("-", 54)) + + for {t_low, t_high, t_label} <- time_blocks, gradient_class <- ["strong", "weak"] do + rows = + Enum.filter(data_10g, fn row -> + row.utc_hour >= t_low and row.utc_hour < t_high and + row.avg_gradient != nil and + if gradient_class == "strong", + do: row.avg_gradient < -100, + else: row.avg_gradient >= -100 + end) + + n = length(rows) + med = if n > 0, do: format_num(median(Enum.map(rows, & &1.distance_km)), 1), else: "-" + + IO.puts( + " #{String.pad_trailing(t_label, 24)} #{String.pad_trailing(gradient_class, 10)} #{String.pad_leading(Integer.to_string(n), 6)} #{String.pad_leading(med, 10)}" + ) + end + + # HPBL × Season + IO.puts("\n2. HPBL × Season") + IO.puts(" (shallow = avg < 500m, mid = 500-999m, deep = avg >= 1000m)") + IO.puts("") + + seasons = [ + {[12, 1, 2], "Winter (Dec-Feb)"}, + {[3, 4, 5], "Spring (Mar-May)"}, + {[6, 7, 8], "Summer (Jun-Aug)"}, + {[9, 10, 11], "Fall (Sep-Nov)"} + ] + + hpbl_classes = [ + {"shallow", fn h -> h < 500 end}, + {"mid", fn h -> h >= 500 and h < 1000 end}, + {"deep", fn h -> h >= 1000 end} + ] + + IO.puts( + " #{String.pad_trailing("Season", 22)} #{String.pad_trailing("HPBL", 10)} #{String.pad_leading("n", 6)} #{String.pad_leading("median km", 10)}" + ) + + IO.puts(" " <> String.duplicate("-", 52)) + + for {months, s_label} <- seasons, {h_label, h_filter} <- hpbl_classes do + rows = + Enum.filter(data_10g, fn row -> + row.month in months and row.avg_hpbl != nil and h_filter.(row.avg_hpbl) + end) + + n = length(rows) + med = if n > 0, do: format_num(median(Enum.map(rows, & &1.distance_km)), 1), else: "-" + + IO.puts( + " #{String.pad_trailing(s_label, 22)} #{String.pad_trailing(h_label, 10)} #{String.pad_leading(Integer.to_string(n), 6)} #{String.pad_leading(med, 10)}" + ) + end + + # Ducting detection vs distance + IO.puts("\n3. Ducting Detection vs Distance") + IO.puts("") + + for ducting_val <- [true, false] do + rows = Enum.filter(data_10g, &(&1.either_ducting == ducting_val)) + n = length(rows) + + if n > 0 do + distances = Enum.map(rows, & &1.distance_km) + + IO.puts( + " Ducting #{if ducting_val, do: "YES", else: "NO "}: n=#{n}, " <> + "median=#{format_num(median(distances), 1)} km, " <> + "p25=#{format_num(percentile(distances, 25), 1)} km, " <> + "p75=#{format_num(percentile(distances, 75), 1)} km, " <> + "max=#{format_num(Enum.max(distances), 1)} km" + ) + else + IO.puts(" Ducting #{if ducting_val, do: "YES", else: "NO "}: n=0") + end + end + + IO.puts("") + end + + # --------------------------------------------------------------------------- + # Statistics helpers + # --------------------------------------------------------------------------- + + defp spearman(pairs) when length(pairs) < 3, do: nil + + defp spearman(pairs) do + n = length(pairs) + + x_vals = Enum.map(pairs, &elem(&1, 0)) + y_vals = Enum.map(pairs, &elem(&1, 1)) + + x_ranks = rank(x_vals) + y_ranks = rank(y_vals) + + d_squared = + x_ranks + |> Enum.zip(y_ranks) + |> Enum.map(fn {rx, ry} -> + d = rx - ry + d * d + end) + |> Enum.sum() + + # Standard Spearman formula: 1 - (6 * sum(d^2)) / (n * (n^2 - 1)) + 1.0 - 6.0 * d_squared / (n * (n * n - 1)) + end + + defp rank(values) do + indexed = + values + |> Enum.with_index() + |> Enum.sort_by(&elem(&1, 0)) + + # Assign average ranks for ties + ranked = + indexed + |> Enum.chunk_by(&elem(&1, 0)) + |> Enum.flat_map(fn group -> + start_rank = length(ranked_before(indexed, group)) + avg_rank = start_rank + (length(group) - 1) / 2.0 + 1 + + Enum.map(group, fn {_val, orig_idx} -> {orig_idx, avg_rank} end) + end) + + ranked + |> Enum.sort_by(&elem(&1, 0)) + |> Enum.map(&elem(&1, 1)) + end + + defp ranked_before(all_sorted, group) do + first_val = group |> hd() |> elem(0) + Enum.take_while(all_sorted, fn {v, _} -> v < first_val end) + end + + defp median([]), do: 0.0 + + defp median(values) do + sorted = Enum.sort(values) + n = length(sorted) + mid = div(n, 2) + + if rem(n, 2) == 0 do + (Enum.at(sorted, mid - 1) + Enum.at(sorted, mid)) / 2.0 + else + Enum.at(sorted, mid) / 1.0 + end + end + + defp percentile([], _p), do: 0.0 + + defp percentile(values, p) do + sorted = Enum.sort(values) + n = length(sorted) + k = p / 100.0 * (n - 1) + f = trunc(k) + c = k - f + + if f + 1 < n do + Enum.at(sorted, f) * (1 - c) + Enum.at(sorted, f + 1) * c + else + Enum.at(sorted, f) / 1.0 + end + end + + defp to_float(val) when is_float(val), do: val + defp to_float(val) when is_integer(val), do: val / 1.0 + defp to_float(%Decimal{} = val), do: Decimal.to_float(val) + defp to_float(val), do: val + + # --------------------------------------------------------------------------- + # Formatting helpers + # --------------------------------------------------------------------------- + + defp format_band(band) when is_float(band), do: format_band(trunc(band)) + defp format_band(10_000), do: "10 GHz" + defp format_band(24_000), do: "24 GHz" + defp format_band(47_000), do: "47 GHz" + defp format_band(75_000), do: "75 GHz" + defp format_band(band), do: "#{band} MHz" + + defp format_num(nil, _decimals), do: "N/A" + defp format_num(val, decimals), do: :erlang.float_to_binary(val / 1.0, decimals: decimals) +end