defmodule Localiser.Localisation.Calibration do
  @moduledoc false

  # Returns [{rssi, is_outlier}] using Tukey IQR fences on the full sample list.
  def classify_outliers([]), do: []

  def classify_outliers(samples) do
    {lo, hi} = iqr_fences(samples)
    Enum.map(samples, fn rssi -> {rssi, rssi < lo or rssi > hi} end)
  end

  # Live hint: is this single reading an outlier relative to the samples collected so far?
  # Returns false when fewer than 5 existing samples (not enough signal).
  def outlier?(_reading, existing) when length(existing) < 5, do: false

  def outlier?(reading, existing) do
    {lo, hi} = iqr_fences(existing)
    reading < lo or reading > hi
  end

  # OLS regression for RSSI = A - 10n * log10(d).
  # stages :: [%{distance: float, mean_rssi: float}]
  # Requires at least 2 stages with distinct distances.
  # Returns {:ok, {rssi_ref :: integer, path_loss_exp :: float}} or {:error, :insufficient_data}.
  def least_squares(stages) when length(stages) < 2, do: {:error, :insufficient_data}

  def least_squares(stages) do
    points = Enum.map(stages, fn %{distance: d, mean_rssi: rssi} ->
      {:math.log10(d), rssi}
    end)

    xs = Enum.map(points, &elem(&1, 0))
    ys = Enum.map(points, &elem(&1, 1))

    x_bar = mean(xs)
    y_bar = mean(ys)

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

    var_x = Enum.reduce(xs, 0.0, fn x, acc -> acc + (x - x_bar) * (x - x_bar) end)

    if var_x == 0.0 do
      {:error, :insufficient_data}
    else
      beta = cov_xy / var_x
      a = y_bar - beta * x_bar

      path_loss_exp = -beta / 10.0
      rssi_ref = round(a)

      {:ok, {rssi_ref, path_loss_exp}}
    end
  end

  # --- Private ---

  defp iqr_fences(samples) do
    sorted = Enum.sort(samples)
    n = length(sorted)

    q1 = percentile(sorted, n, 0.25)
    q3 = percentile(sorted, n, 0.75)
    iqr = q3 - q1

    {q1 - 1.5 * iqr, q3 + 1.5 * iqr}
  end

  defp percentile(sorted, n, p) do
    idx = p * (n - 1)
    lo = floor(idx)
    hi = ceil(idx)

    if lo == hi do
      Enum.at(sorted, lo) * 1.0
    else
      frac = idx - lo
      Enum.at(sorted, lo) * (1 - frac) + Enum.at(sorted, hi) * frac
    end
  end

  defp mean(list) do
    Enum.sum(list) / length(list)
  end
end