prop/lib/microwaveprop/backtest.ex

defmodule Microwaveprop.Backtest do
  @moduledoc """
  Evaluates a propagation feature function against the historical QSO
  corpus and a random-time baseline.

  A feature function has the shape `(lat, lon, valid_time) -> float | nil`.
  `nil` means "no data" and is excluded from the distribution stats.

  Given a feature, `evaluate/2` pulls QSOs (up to `:sample_size`) with a
  known `pos1`, applies the feature at the station1 location and the
  QSO timestamp, and then does the same for `:baseline_size` random
  (lat, lon, time) samples drawn in the same geographic and temporal
  neighborhood. Comparing the two distributions tells us whether the
  feature carries information about when propagation is good.

  The baseline is deliberately a matched sample: each random draw
  picks a real QSO location and perturbs its timestamp by a uniform
  ±30 days, so seasonal and diurnal effects are not a free variable.
  """

  import Ecto.Query

  alias Microwaveprop.Backtest.Distribution
  alias Microwaveprop.Radio.Contact
  alias Microwaveprop.Repo

  @type feature :: (float, float, DateTime.t() -> float | nil)

  defmodule Distribution do
    @moduledoc "Summary statistics for a collection of feature values."

    defstruct count: 0, mean: nil, stddev: nil, p50: nil, p90: nil, min: nil, max: nil

    @type t :: %__MODULE__{
            count: non_neg_integer(),
            mean: float | nil,
            stddev: float | nil,
            p50: float | nil,
            p90: float | nil,
            min: float | nil,
            max: float | nil
          }

    def from_values([]), do: %__MODULE__{count: 0}

    def from_values(values) when is_list(values) do
      sorted = Enum.sort(values)
      n = length(sorted)
      mean = Enum.sum(sorted) / n

      variance =
        if n > 1 do
          sorted
          |> Enum.reduce(0.0, fn v, acc -> acc + :math.pow(v - mean, 2) end)
          |> Kernel./(n - 1)
        else
          0.0
        end

      %__MODULE__{
        count: n,
        mean: mean,
        stddev: :math.sqrt(variance),
        p50: percentile(sorted, 0.5),
        p90: percentile(sorted, 0.9),
        min: List.first(sorted),
        max: List.last(sorted)
      }
    end

    defp percentile(sorted, fraction) do
      n = length(sorted)
      idx = trunc(Float.floor(fraction * (n - 1)))
      Enum.at(sorted, idx)
    end
  end

  defmodule Report do
    @moduledoc "Top-level result of Backtest.evaluate/2."

    defstruct feature_name: "anonymous",
              qso_count: 0,
              baseline_count: 0,
              qso_distribution: %Distribution{},
              baseline_distribution: %Distribution{}

    @type t :: %__MODULE__{
            feature_name: String.t(),
            qso_count: non_neg_integer(),
            baseline_count: non_neg_integer(),
            qso_distribution: Distribution.t(),
            baseline_distribution: Distribution.t()
          }
  end

  @doc """
  Generates `n` matched random baseline samples.

  Each sample picks a real QSO location and perturbs its timestamp by
  a uniform ±30 days. This controls for the seasonal and diurnal
  distribution of contacts so the baseline isn't trivially distinguishable
  from the QSO sample by time-of-year effects.

  Returns a list of `{lat, lon, %DateTime{}}` triples. Returns `[]` if
  there are no contacts with a known position.

  ## Options

    * `:sample_size` - maximum number of contacts to draw from (default: 5000).
  """
  @spec random_baseline(non_neg_integer(), keyword) :: [{float, float, DateTime.t()}]
  def random_baseline(n, opts \\ []) when is_integer(n) and n >= 0 do
    sample_size = Keyword.get(opts, :sample_size, 5000)

    sample_size
    |> load_contacts()
    |> baseline_samples(n)
  end

  @distance_bins [
    {"0-100", 0.0, 100.0},
    {"100-250", 100.0, 250.0},
    {"250-500", 250.0, 500.0},
    {"500-1000", 500.0, 1000.0},
    {"1000+", 1000.0, :infinity}
  ]

  @doc """
  Bins QSOs by `distance_km` and reports the feature distribution within each bin.

  A feature that carries information about propagation quality should
  show a monotonically increasing mean across distance bins: longer
  contacts imply better propagation, so the feature should be higher
  where the contacts reach further.

  Returns a map keyed by bin label (`"0-100"`, `"100-250"`, …) whose
  values are `Distribution` structs.
  """
  @spec lift_by_distance(feature, keyword) :: %{String.t() => Distribution.t()}
  def lift_by_distance(feature, opts \\ []) when is_function(feature, 3) do
    sample_size = Keyword.get(opts, :sample_size, 5000)

    contacts = load_contacts(sample_size)

    contacts
    |> Enum.map(fn contact ->
      {distance_bin(contact.distance_km), eval_feature_for_contact(feature, contact)}
    end)
    |> Enum.reject(fn {bin, value} -> is_nil(bin) or is_nil(value) end)
    |> Enum.group_by(fn {bin, _} -> bin end, fn {_, value} -> value end)
    |> Map.new(fn {bin, values} -> {bin, Distribution.from_values(values)} end)
    |> fill_missing_bins()
  end

  defp distance_bin(nil), do: nil

  defp distance_bin(%Decimal{} = d), do: distance_bin(Decimal.to_float(d))

  defp distance_bin(km) when is_number(km) do
    Enum.find_value(@distance_bins, fn {label, lo, hi} ->
      cond do
        hi == :infinity and km >= lo -> label
        km >= lo and km < hi -> label
        true -> nil
      end
    end)
  end

  defp fill_missing_bins(bins) do
    Enum.reduce(@distance_bins, bins, fn {label, _lo, _hi}, acc ->
      Map.put_new(acc, label, %Distribution{})
    end)
  end

  @doc """
  Groups feature values by `band` and reports the distribution per band.

  Useful to spot features whose lift is band-dependent: a duct-geometry
  feature should carry more information on 24+ GHz than on 10 GHz, and a
  humidity feature should do the opposite.

  Returns a map keyed by `Decimal` band (MHz).
  """
  @spec lift_by_band(feature, keyword) :: %{Decimal.t() => Distribution.t()}
  def lift_by_band(feature, opts \\ []) when is_function(feature, 3) do
    sample_size = Keyword.get(opts, :sample_size, 5000)

    contacts = load_contacts(sample_size)

    contacts
    |> Enum.map(fn contact -> {contact.band, eval_feature_for_contact(feature, contact)} end)
    |> Enum.reject(fn {band, value} -> is_nil(band) or is_nil(value) end)
    |> Enum.group_by(fn {band, _} -> band end, fn {_, value} -> value end)
    |> Map.new(fn {band, values} -> {band, Distribution.from_values(values)} end)
  end

  @doc """
  Runs a feature function over QSOs and a random baseline, returning
  a `Report` struct with matched distribution statistics.

  ## Options

    * `:sample_size` - maximum number of QSOs to evaluate (default: 5000).
    * `:baseline_size` - number of random-time baseline samples (default: same as `:sample_size`).
    * `:feature_name` - label for the report (default: "anonymous").
  """
  @spec evaluate(feature, keyword) :: Report.t()
  def evaluate(feature, opts \\ []) when is_function(feature, 3) do
    sample_size = Keyword.get(opts, :sample_size, 5000)
    baseline_size = Keyword.get(opts, :baseline_size, sample_size)
    feature_name = Keyword.get(opts, :feature_name, "anonymous")

    contacts = load_contacts(sample_size)

    qso_values =
      contacts
      |> Enum.map(fn c -> eval_feature_for_contact(feature, c) end)
      |> Enum.reject(&is_nil/1)

    baseline_values =
      contacts
      |> baseline_samples(baseline_size)
      |> Enum.map(fn {lat, lon, time} -> feature.(lat, lon, time) end)
      |> Enum.reject(&is_nil/1)

    %Report{
      feature_name: feature_name,
      qso_count: length(contacts),
      baseline_count: baseline_size,
      qso_distribution: Distribution.from_values(qso_values),
      baseline_distribution: Distribution.from_values(baseline_values)
    }
  end

  @doc """
  Runs multiple features against the same QSO sample and returns a list
  of summary maps — one per feature — suitable for a consolidated table.

  `features` is a map of `%{name => fun/3}`.

  Each entry includes the feature name, QSO and baseline distribution
  stats, and a pass/fail gate based on whether the feature has data.
  """
  @spec consolidated_report(%{String.t() => feature}, keyword) :: [map()]
  def consolidated_report(features, opts \\ []) when is_map(features) do
    sample_size = Keyword.get(opts, :sample_size, 5000)
    baseline_size = Keyword.get(opts, :baseline_size, sample_size)

    contacts = load_contacts(sample_size)

    baseline_samples = baseline_samples(contacts, baseline_size)

    features
    |> Enum.map(fn {name, fun} ->
      qso_values =
        contacts
        |> Enum.map(&eval_feature_for_contact(fun, &1))
        |> Enum.reject(&is_nil/1)

      baseline_values =
        baseline_samples
        |> Enum.map(fn {lat, lon, time} -> fun.(lat, lon, time) end)
        |> Enum.reject(&is_nil/1)

      qso_dist = Distribution.from_values(qso_values)
      baseline_dist = Distribution.from_values(baseline_values)

      %{
        feature_name: name,
        qso_distribution: qso_dist,
        baseline_distribution: baseline_dist,
        gate: if(qso_dist.count > 0, do: :pass, else: :no_data)
      }
    end)
    |> Enum.sort_by(& &1.feature_name)
  end

  @doc """
  Renders the output of `consolidated_report/2` as a single Markdown table.
  """
  @spec to_consolidated_markdown([map()]) :: String.t()
  def to_consolidated_markdown(results) when is_list(results) do
    header = "| Feature | QSO N | QSO Mean | QSO p50 | Baseline N | Baseline Mean | Baseline p50 | Gate |\n"
    sep = "|---|---|---|---|---|---|---|---|\n"

    body =
      Enum.map_join(results, "", fn r ->
        q = r.qso_distribution
        b = r.baseline_distribution
        gate_str = if r.gate == :pass, do: "PASS", else: "NO DATA"

        "| #{r.feature_name} | #{q.count} | #{fmt(q.mean)} | #{fmt(q.p50)} | #{b.count} | #{fmt(b.mean)} | #{fmt(b.p50)} | #{gate_str} |\n"
      end)

    IO.iodata_to_binary([
      "# Consolidated Backtest Report\n\n",
      header,
      sep,
      body,
      "\n"
    ])
  end

  @doc """
  Renders a `Report` (plus optional distance and band tables) as Markdown.

  This is the shape the `mix backtest` task prints to stdout and the
  shape we write into `priv/backtest_reports/` for version control.
  """
  @spec to_markdown(Report.t(), keyword) :: String.t()
  def to_markdown(%Report{} = report, opts \\ []) do
    distance_bins = Keyword.get(opts, :distance_bins)
    band_stats = Keyword.get(opts, :band_stats)

    iodata = [
      "# Backtest: #{report.feature_name}\n\n",
      "QSO sample: #{report.qso_count}  \n",
      "Baseline sample: #{report.baseline_count}\n\n",
      "## Matched distribution\n\n",
      distribution_table([
        {"QSO times", report.qso_distribution},
        {"Random baseline", report.baseline_distribution}
      ]),
      distance_section(distance_bins),
      band_section(band_stats)
    ]

    IO.iodata_to_binary(iodata)
  end

  defp distribution_table(rows) do
    header = "| Set | N | Mean | Stddev | p50 | p90 | Min | Max |\n"
    sep = "|---|---|---|---|---|---|---|---|\n"

    body =
      Enum.map_join(rows, "", fn {label, dist} ->
        "| #{label} | #{dist.count} | #{fmt(dist.mean)} | #{fmt(dist.stddev)} | #{fmt(dist.p50)} | #{fmt(dist.p90)} | #{fmt(dist.min)} | #{fmt(dist.max)} |\n"
      end)

    [header, sep, body, "\n"]
  end

  defp distance_section(nil), do: []

  defp distance_section(bins) do
    rows =
      Enum.map(@distance_bins, fn {label, _lo, _hi} ->
        {label, Map.get(bins, label, %Distribution{})}
      end)

    ["## Lift by distance (km)\n\n", distribution_table(rows)]
  end

  defp band_section(nil), do: []

  defp band_section(bands) do
    rows =
      bands
      |> Enum.sort_by(fn {band, _} -> Decimal.to_float(band) end)
      |> Enum.map(fn {band, dist} -> {"#{band} MHz", dist} end)

    ["## Lift by band\n\n", distribution_table(rows)]
  end

  defp fmt(nil), do: "—"
  defp fmt(n) when is_float(n), do: :erlang.float_to_binary(n, decimals: 3)
  defp fmt(n), do: to_string(n)

  # Load contacts with a known position, newest first.
  defp load_contacts(limit) do
    Contact
    |> where([c], not is_nil(c.pos1))
    |> order_by([c], desc: c.qso_timestamp)
    |> limit(^limit)
    |> Repo.all()
  end

  defp eval_feature_for_contact(feature, %Contact{pos1: pos, qso_timestamp: ts}) do
    case pos_to_latlon(pos) do
      {lat, lon} -> feature.(lat, lon, ts)
      nil -> nil
    end
  end

  defp eval_feature_for_contact(_feature, _contact), do: nil

  # pos1 is sometimes stored with "lon" and sometimes with "lng".
  defp pos_to_latlon(%{"lat" => lat, "lon" => lon}) when is_number(lat) and is_number(lon),
    do: {lat * 1.0, lon * 1.0}

  defp pos_to_latlon(%{"lat" => lat, "lng" => lon}) when is_number(lat) and is_number(lon),
    do: {lat * 1.0, lon * 1.0}

  defp pos_to_latlon(_), do: nil

  # Matched random baseline: pick a random contact, perturb its
  # timestamp ±30 days. Empty contact list → empty baseline.
  defp baseline_samples([], _n), do: []

  defp baseline_samples(_contacts, 0), do: []

  defp baseline_samples(contacts, n) do
    usable = Enum.filter(contacts, fn c -> pos_to_latlon(c.pos1) end)

    if usable == [] do
      []
    else
      for _ <- 1..n do
        contact = Enum.random(usable)
        {lat, lon} = pos_to_latlon(contact.pos1)
        offset_seconds = :rand.uniform(60 * 24 * 3600) - 30 * 24 * 3600
        time = DateTime.add(contact.qso_timestamp, offset_seconds, :second)
        {lat, lon, time}
      end
    end
  end
end