prop/lib_ml/calibration.ex

defmodule Microwaveprop.Propagation.Calibration do
  @moduledoc """
  Pure analysis functions for the algorithm-vs-ML-vs-reality
  comparison reports written by `mix prop.compare`.

  Every function here takes plain maps and lists, so the Mix task can
  load data however it wants and these helpers stay easy to test
  without a database.

  Sample shape (one map per QSO):

      %{
        algorithm_score: 0..100,    # composite score from Scorer
        ml_score: 0..100,           # ML model prediction (or nil if unavailable)
        distance_km: float(),       # actual contact distance — the empirical signal
        band: integer(),            # MHz
        # optional, copied through for outlier inspection:
        contact_id: term(), conditions: map()
      }
  """

  @score_buckets [
    {"0-20", 0, 20},
    {"20-40", 20, 40},
    {"40-60", 40, 60},
    {"60-80", 60, 80},
    {"80-100", 80, 101}
  ]

  @doc """
  Spearman rank correlation. Returns `nil` for fewer than 3 pairs.
  Tied ranks are averaged. Returns `0` (by convention) when one input
  has zero variance.
  """
  @spec spearman([number()], [number()]) :: float() | nil
  def spearman(xs, ys) when length(xs) != length(ys), do: nil
  def spearman(xs, _) when length(xs) < 3, do: nil

  def spearman(xs, ys) do
    n = length(xs)
    rx = ranks(xs)
    ry = ranks(ys)

    if zero_variance?(rx) or zero_variance?(ry) do
      0.0
    else
      sum_d_sq =
        rx
        |> Enum.zip(ry)
        |> Enum.reduce(0.0, fn {a, b}, acc ->
          d = a - b
          acc + d * d
        end)

      1.0 - 6.0 * sum_d_sq / (n * (n * n - 1))
    end
  end

  @doc "Root-mean-squared error between two equal-length numeric lists. Nil for empty."
  @spec rmse([number()], [number()]) :: float() | nil
  def rmse([], []), do: nil

  def rmse(xs, ys) when length(xs) == length(ys) do
    n = length(xs)

    sum_sq =
      xs
      |> Enum.zip(ys)
      |> Enum.reduce(0.0, fn {x, y}, acc ->
        d = x - y
        acc + d * d
      end)

    :math.sqrt(sum_sq / n)
  end

  @doc """
  Splits samples into the five fixed score buckets ([0-20, 20-40, ...,
  80-100]) and reports `n` + median observed distance for each. Always
  returns 5 entries, even when buckets are empty.

  A well-calibrated scorer shows monotonically increasing median
  distance across the buckets.
  """
  @spec bucket_by_score([map()], atom()) :: [%{bucket: String.t(), n: non_neg_integer(), median_distance_km: float() | nil}]
  def bucket_by_score(samples, score_key) do
    Enum.map(@score_buckets, fn {label, low, high} ->
      in_bucket =
        Enum.filter(samples, fn s ->
          score = Map.get(s, score_key)
          is_number(score) and score >= low and score < high
        end)

      %{
        bucket: label,
        n: length(in_bucket),
        median_distance_km: in_bucket |> Enum.map(& &1.distance_km) |> median()
      }
    end)
  end

  @doc """
  One-band summary: row count, mean scores, algorithm-vs-ML RMSE, and
  Spearman correlations of each score against actual contact distance.
  """
  @spec band_summary([map()]) :: %{
          n: non_neg_integer(),
          mean_algorithm_score: float() | nil,
          mean_ml_score: float() | nil,
          alg_ml_rmse: float() | nil,
          alg_distance_spearman: float() | nil,
          ml_distance_spearman: float() | nil
        }
  def band_summary([]) do
    %{
      n: 0,
      mean_algorithm_score: nil,
      mean_ml_score: nil,
      alg_ml_rmse: nil,
      alg_distance_spearman: nil,
      ml_distance_spearman: nil
    }
  end

  def band_summary(samples) do
    alg = Enum.map(samples, & &1.algorithm_score)
    ml = Enum.map(samples, & &1.ml_score)
    distances = Enum.map(samples, & &1.distance_km)

    %{
      n: length(samples),
      mean_algorithm_score: mean(alg),
      mean_ml_score: mean(ml),
      alg_ml_rmse: rmse(alg, ml),
      alg_distance_spearman: spearman(alg, distances),
      ml_distance_spearman: spearman(ml, distances)
    }
  end

  @doc "Top `n` samples by absolute (algorithm − ML) score difference, descending."
  @spec disagreements([map()], non_neg_integer()) :: [map()]
  def disagreements(samples, n) do
    samples
    |> Enum.map(fn s -> Map.put(s, :divergence, abs(s.algorithm_score - s.ml_score)) end)
    |> Enum.sort_by(& &1.divergence, :desc)
    |> Enum.take(n)
  end

  @doc """
  Drift detector. Two modes:

    * `:ml_drift` — passes the current alg-vs-ML RMSE; flags when it
      exceeds `:threshold` (default 8.0 score points).
    * `:algorithm_drift` — passes the current Spearman correlation
      between algorithm score and contact distance plus a `:history`
      list of past correlations; flags when current is more than
      `:margin` (default 0.10) below the history median.
  """
  @spec detect_drift(:ml_drift | :algorithm_drift, number() | nil, keyword()) ::
          :ok | {:drift, String.t()}
  def detect_drift(:ml_drift, current_rmse, opts) do
    threshold = Keyword.get(opts, :threshold, 8.0)

    cond do
      is_nil(current_rmse) -> :ok
      current_rmse > threshold ->
        {:drift, "ML model has drifted: alg-vs-ML RMSE = #{Float.round(current_rmse, 2)} > #{threshold}"}
      true -> :ok
    end
  end

  def detect_drift(:algorithm_drift, current_corr, opts) do
    history = Keyword.get(opts, :history, [])
    margin = Keyword.get(opts, :margin, 0.10)

    cond do
      is_nil(current_corr) -> :ok
      history == [] -> :ok
      current_corr < median(history) - margin ->
        {:drift,
         "Algorithm correlation dropped: current #{format(current_corr)} vs history median " <>
           "#{format(median(history))} (margin #{margin})"}

      true ->
        :ok
    end
  end

  ## ── helpers ──────────────────────────────────────────────────────

  defp ranks(values) do
    values
    |> Enum.with_index()
    |> Enum.sort_by(&elem(&1, 0))
    |> Enum.chunk_by(&elem(&1, 0))
    |> Enum.flat_map_reduce(1, fn group, start ->
      avg_rank = start + (length(group) - 1) / 2.0
      tagged = Enum.map(group, fn {_v, idx} -> {idx, avg_rank} end)
      {tagged, start + length(group)}
    end)
    |> elem(0)
    |> Enum.sort_by(&elem(&1, 0))
    |> Enum.map(&elem(&1, 1))
  end

  defp zero_variance?(list), do: list |> Enum.uniq() |> length() == 1

  defp mean([]), do: nil
  defp mean(list), do: Enum.sum(list) / length(list)

  defp median([]), do: nil

  defp median(list) do
    sorted = Enum.sort(list)
    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 format(n) when is_number(n), do: :erlang.float_to_binary(n * 1.0, decimals: 3)
end