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,

      -- Longitude for solar time
      ((q.pos1->>'lng')::float + (q.pos2->>'lng')::float) / 2.0 AS avg_longitude

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

  defp derive_solar_hour(row) do
    case {row[:utc_hour], row[:avg_longitude]} do
      {h, lng} when is_number(h) and is_number(lng) ->
        solar = :math.fmod(h + lng / 15 + 24, 24)
        Map.put(row, :solar_hour, solar)

      _ ->
        Map.put(row, :solar_hour, nil)
    end
  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"},
      {:solar_hour, "Solar 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"}
       ]},
      {:solar_hour, "Solar Hour",
       [
         {0, 3, "00-02 Solar (night)"},
         {3, 6, "03-05 Solar (pre-dawn)"},
         {6, 9, "06-08 Solar (dawn)"},
         {9, 12, "09-11 Solar (morning)"},
         {12, 15, "12-14 Solar (midday)"},
         {15, 18, "15-17 Solar (afternoon)"},
         {18, 21, "18-20 Solar (evening)"},
         {21, 24, "21-23 Solar (night)"}
       ]},
      {: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