Add mix propagation.analyze task for QSO-HRRR correlation analysis
This commit is contained in:
parent
14e86c4ddb
commit
4bf9d0b4b0
1 changed files with 612 additions and 0 deletions
612
lib/mix/tasks/propagation_analyze.ex
Normal file
612
lib/mix/tasks/propagation_analyze.ex
Normal file
|
|
@ -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
|
||||
Loading…
Add table
Reference in a new issue