From ab04cb91682347b45dca4547831595574d3a9f53 Mon Sep 17 00:00:00 2001
From: Graham McIntire <graham@mcintire.me>
Date: Thu, 9 Apr 2026 16:10:45 -0500
Subject: [PATCH] Phase 0: backtest harness

Add Microwaveprop.Backtest: a feature-evaluation framework that runs
a (lat, lon, valid_time) -> float function over the historical QSO
corpus and a matched random-time baseline, reporting distribution
statistics, distance-binned lift, and band-stratified lift.

Adds four baseline feature wrappers around the current scorer inputs
(naive_gradient, td_depression, time_of_day, pressure), a mix backtest
CLI, and the first set of baseline reports under priv/backtest_reports
so downstream phases have a frozen reference point to compare against.
---
 lib/microwaveprop/backtest.ex                 | 348 ++++++++++++++++++
 lib/microwaveprop/backtest/features.ex        |  84 +++++
 lib/mix/tasks/backtest.ex                     |  81 ++++
 .../baseline_naive_gradient.md                |  31 ++
 priv/backtest_reports/baseline_pressure.md    |  31 ++
 .../baseline_td_depression.md                 |  31 ++
 priv/backtest_reports/baseline_time_of_day.md |  31 ++
 test/microwaveprop/backtest/features_test.exs |  62 ++++
 test/microwaveprop/backtest_test.exs          | 170 +++++++++
 9 files changed, 869 insertions(+)
 create mode 100644 lib/microwaveprop/backtest.ex
 create mode 100644 lib/microwaveprop/backtest/features.ex
 create mode 100644 lib/mix/tasks/backtest.ex
 create mode 100644 priv/backtest_reports/baseline_naive_gradient.md
 create mode 100644 priv/backtest_reports/baseline_pressure.md
 create mode 100644 priv/backtest_reports/baseline_td_depression.md
 create mode 100644 priv/backtest_reports/baseline_time_of_day.md
 create mode 100644 test/microwaveprop/backtest/features_test.exs
 create mode 100644 test/microwaveprop/backtest_test.exs

diff --git a/lib/microwaveprop/backtest.ex b/lib/microwaveprop/backtest.ex
new file mode 100644
index 00000000..108d1cac
--- /dev/null
+++ b/lib/microwaveprop/backtest.ex
@@ -0,0 +1,348 @@
+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 """
+  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
diff --git a/lib/microwaveprop/backtest/features.ex b/lib/microwaveprop/backtest/features.ex
new file mode 100644
index 00000000..2cc0d226
--- /dev/null
+++ b/lib/microwaveprop/backtest/features.ex
@@ -0,0 +1,84 @@
+defmodule Microwaveprop.Backtest.Features do
+  @moduledoc """
+  Named feature functions for use with `Microwaveprop.Backtest`.
+
+  Every feature has the shape `(lat, lon, valid_time) -> float | nil`
+  and is named after the physical quantity it represents. These are
+  the "known baselines" the plan refers to: wrappers around the
+  existing scorer's inputs so we can measure the lift of new features
+  against them on an apples-to-apples basis.
+
+  ## Contract
+
+  - Return a `float` when the underlying data is available.
+  - Return `nil` when there's no HRRR profile within the usual match
+    window (`Weather.find_nearest_hrrr/3` returning nil).
+  - Never raise — bad inputs should produce `nil`, not crashes. The
+    backtest harness runs these across tens of thousands of calls and
+    a raise on one point kills the whole report.
+  """
+
+  alias Microwaveprop.Weather
+
+  @doc """
+  Minimum refractivity gradient from the nearest HRRR profile.
+
+  This is the scalar the current scorer uses. More negative is better
+  (stronger ducting potential). We return the raw gradient; the
+  backtest harness handles binning and summarizing.
+  """
+  @spec naive_gradient(float, float, DateTime.t()) :: float | nil
+  def naive_gradient(lat, lon, valid_time) do
+    with %{min_refractivity_gradient: grad} when is_float(grad) <-
+           Weather.find_nearest_hrrr(lat, lon, valid_time) do
+      grad
+    else
+      _ -> nil
+    end
+  end
+
+  @doc """
+  Dewpoint depression (T - Td) at the surface, in °C.
+
+  Lower depression means higher relative humidity. For 10 GHz the
+  existing scorer treats this as beneficial; for 24+ GHz it's harmful.
+  """
+  @spec td_depression(float, float, DateTime.t()) :: float | nil
+  def td_depression(lat, lon, valid_time) do
+    with %{surface_temp_c: t, surface_dewpoint_c: td} when is_float(t) and is_float(td) <-
+           Weather.find_nearest_hrrr(lat, lon, valid_time) do
+      t - td
+    else
+      _ -> nil
+    end
+  end
+
+  @doc """
+  Time-of-day feature: hours since midnight UTC, as a float in [0, 24).
+
+  A flat-by-time baseline against which diurnal lift is measured. The
+  existing scorer collapses this to a band-dependent shape; the
+  backtest treats it as raw UTC hour so we can see the shape directly
+  in the distribution.
+  """
+  @spec time_of_day(float, float, DateTime.t()) :: float
+  def time_of_day(_lat, _lon, valid_time) do
+    valid_time.hour + valid_time.minute / 60.0
+  end
+
+  @doc """
+  Surface pressure in hPa from the nearest HRRR profile.
+
+  Used as the baseline the plan predicts `ParallelToFront` (Phase 5)
+  will replace.
+  """
+  @spec pressure(float, float, DateTime.t()) :: float | nil
+  def pressure(lat, lon, valid_time) do
+    with %{surface_pressure_mb: p} when is_float(p) <-
+           Weather.find_nearest_hrrr(lat, lon, valid_time) do
+      p
+    else
+      _ -> nil
+    end
+  end
+end
diff --git a/lib/mix/tasks/backtest.ex b/lib/mix/tasks/backtest.ex
new file mode 100644
index 00000000..0ff2efae
--- /dev/null
+++ b/lib/mix/tasks/backtest.ex
@@ -0,0 +1,81 @@
+defmodule Mix.Tasks.Backtest do
+  @shortdoc "Evaluate a propagation feature against the QSO corpus"
+  @moduledoc """
+  Runs `Microwaveprop.Backtest.evaluate/2` (plus the distance and band
+  breakdowns) for a named feature function and prints a Markdown report
+  to stdout.
+
+  ## Usage
+
+      mix backtest --feature Microwaveprop.Backtest.Features.naive_gradient
+      mix backtest --feature naive_gradient
+      mix backtest --feature naive_gradient --sample 1000 --out priv/backtest_reports/naive.md
+
+  ## Options
+
+    * `--feature` (required) — fully-qualified `Module.function` or a
+      short name that lives on `Microwaveprop.Backtest.Features`.
+    * `--sample` — max number of QSOs to evaluate (default: 5000).
+    * `--baseline` — random-baseline sample size (default: same as `--sample`).
+    * `--out` — optional file path to write the report to in addition
+      to printing it. Useful for saving baseline reports into
+      `priv/backtest_reports/`.
+  """
+  use Mix.Task
+
+  alias Microwaveprop.Backtest
+
+  @impl Mix.Task
+  def run(argv) do
+    Mix.Task.run("app.start")
+
+    {opts, _, _} =
+      OptionParser.parse(argv,
+        switches: [feature: :string, sample: :integer, baseline: :integer, out: :string]
+      )
+
+    feature_spec = Keyword.fetch!(opts, :feature)
+    sample_size = Keyword.get(opts, :sample, 5000)
+    baseline_size = Keyword.get(opts, :baseline, sample_size)
+    out_path = Keyword.get(opts, :out)
+
+    {feature_fun, feature_name} = resolve_feature(feature_spec)
+
+    report =
+      Backtest.evaluate(feature_fun,
+        sample_size: sample_size,
+        baseline_size: baseline_size,
+        feature_name: feature_name
+      )
+
+    distance_bins = Backtest.lift_by_distance(feature_fun, sample_size: sample_size)
+    band_stats = Backtest.lift_by_band(feature_fun, sample_size: sample_size)
+
+    markdown =
+      Backtest.to_markdown(report, distance_bins: distance_bins, band_stats: band_stats)
+
+    IO.puts(markdown)
+
+    if out_path do
+      File.mkdir_p!(Path.dirname(out_path))
+      File.write!(out_path, markdown)
+      Mix.shell().info("Wrote report to #{out_path}")
+    end
+  end
+
+  defp resolve_feature(spec) do
+    case String.split(spec, ".") do
+      [name] ->
+        fun = String.to_atom(name)
+        feature_fun = &apply(Microwaveprop.Backtest.Features, fun, [&1, &2, &3])
+        {feature_fun, "Microwaveprop.Backtest.Features.#{name}"}
+
+      parts ->
+        {fun_name, mod_parts} = List.pop_at(parts, -1)
+        module = Module.concat(mod_parts)
+        fun = String.to_atom(fun_name)
+        feature_fun = &apply(module, fun, [&1, &2, &3])
+        {feature_fun, Enum.join(parts, ".")}
+    end
+  end
+end
diff --git a/priv/backtest_reports/baseline_naive_gradient.md b/priv/backtest_reports/baseline_naive_gradient.md
new file mode 100644
index 00000000..9ef100a4
--- /dev/null
+++ b/priv/backtest_reports/baseline_naive_gradient.md
@@ -0,0 +1,31 @@
+# Backtest: Microwaveprop.Backtest.Features.naive_gradient
+
+QSO sample: 5000  
+Baseline sample: 5000
+
+## Matched distribution
+
+| Set | N | Mean | Stddev | p50 | p90 | Min | Max |
+|---|---|---|---|---|---|---|---|
+| QSO times | 5000 | -113.971 | 42.207 | -105.154 | -69.043 | -331.397 | -35.521 |
+| Random baseline | 151 | -118.002 | 43.405 | -107.528 | -72.860 | -268.079 | -41.047 |
+
+## Lift by distance (km)
+
+| Set | N | Mean | Stddev | p50 | p90 | Min | Max |
+|---|---|---|---|---|---|---|---|
+| 0-100 | 1249 | -114.766 | 43.617 | -107.614 | -66.251 | -321.312 | -35.521 |
+| 100-250 | 2295 | -114.555 | 42.595 | -103.392 | -71.537 | -331.397 | -44.027 |
+| 250-500 | 1288 | -112.728 | 41.507 | -105.238 | -65.613 | -321.312 | -39.224 |
+| 500-1000 | 167 | -109.723 | 29.535 | -115.466 | -69.348 | -206.669 | -45.098 |
+| 1000+ | 1 | -90.610 | 0.000 | -90.610 | -90.610 | -90.610 | -90.610 |
+
+## Lift by band
+
+| Set | N | Mean | Stddev | p50 | p90 | Min | Max |
+|---|---|---|---|---|---|---|---|
+| 10000 MHz | 4444 | -114.037 | 42.734 | -105.238 | -68.535 | -321.312 | -37.154 |
+| 24000 MHz | 465 | -116.766 | 37.728 | -107.614 | -82.149 | -331.397 | -64.361 |
+| 47000 MHz | 77 | -98.886 | 35.255 | -93.461 | -68.535 | -227.800 | -35.521 |
+| 75000 MHz | 14 | -83.273 | 13.340 | -93.461 | -68.535 | -99.652 | -68.535 |
+
diff --git a/priv/backtest_reports/baseline_pressure.md b/priv/backtest_reports/baseline_pressure.md
new file mode 100644
index 00000000..d5a949bb
--- /dev/null
+++ b/priv/backtest_reports/baseline_pressure.md
@@ -0,0 +1,31 @@
+# Backtest: Microwaveprop.Backtest.Features.pressure
+
+QSO sample: 5000  
+Baseline sample: 5000
+
+## Matched distribution
+
+| Set | N | Mean | Stddev | p50 | p90 | Min | Max |
+|---|---|---|---|---|---|---|---|
+| QSO times | 5000 | 983.092 | 29.427 | 988.600 | 1009.900 | 788.600 | 1023.500 |
+| Random baseline | 136 | 981.549 | 30.132 | 988.600 | 1008.900 | 830.600 | 1020.700 |
+
+## Lift by distance (km)
+
+| Set | N | Mean | Stddev | p50 | p90 | Min | Max |
+|---|---|---|---|---|---|---|---|
+| 0-100 | 1249 | 988.560 | 25.473 | 991.800 | 1013.200 | 788.600 | 1023.500 |
+| 100-250 | 2295 | 982.276 | 30.947 | 988.800 | 1006.000 | 790.200 | 1023.500 |
+| 250-500 | 1288 | 979.843 | 30.739 | 986.300 | 1006.100 | 830.900 | 1017.700 |
+| 500-1000 | 167 | 978.470 | 17.374 | 979.100 | 999.700 | 949.800 | 1015.600 |
+| 1000+ | 1 | 982.400 | 0.000 | 982.400 | 982.400 | 982.400 | 982.400 |
+
+## Lift by band
+
+| Set | N | Mean | Stddev | p50 | p90 | Min | Max |
+|---|---|---|---|---|---|---|---|
+| 10000 MHz | 4444 | 983.001 | 28.768 | 988.500 | 1009.800 | 830.900 | 1023.500 |
+| 24000 MHz | 465 | 983.776 | 32.080 | 990.800 | 1010.100 | 830.900 | 1016.900 |
+| 47000 MHz | 77 | 983.048 | 47.411 | 1000.400 | 1012.800 | 788.600 | 1021.100 |
+| 75000 MHz | 14 | 989.464 | 14.483 | 987.300 | 1001.600 | 956.000 | 1001.600 |
+
diff --git a/priv/backtest_reports/baseline_td_depression.md b/priv/backtest_reports/baseline_td_depression.md
new file mode 100644
index 00000000..bbf260ae
--- /dev/null
+++ b/priv/backtest_reports/baseline_td_depression.md
@@ -0,0 +1,31 @@
+# Backtest: Microwaveprop.Backtest.Features.td_depression
+
+QSO sample: 5000  
+Baseline sample: 5000
+
+## Matched distribution
+
+| Set | N | Mean | Stddev | p50 | p90 | Min | Max |
+|---|---|---|---|---|---|---|---|
+| QSO times | 5000 | 6.569 | 5.365 | 5.573 | 12.969 | -0.693 | 26.814 |
+| Random baseline | 161 | 4.719 | 4.367 | 3.534 | 10.096 | -1.258 | 21.762 |
+
+## Lift by distance (km)
+
+| Set | N | Mean | Stddev | p50 | p90 | Min | Max |
+|---|---|---|---|---|---|---|---|
+| 0-100 | 1249 | 5.324 | 4.270 | 4.151 | 10.907 | -0.693 | 25.733 |
+| 100-250 | 2295 | 6.941 | 5.789 | 5.589 | 14.564 | -0.630 | 26.814 |
+| 250-500 | 1288 | 7.007 | 5.510 | 6.054 | 14.450 | -0.568 | 25.733 |
+| 500-1000 | 167 | 7.390 | 3.911 | 7.422 | 12.637 | -0.142 | 15.938 |
+| 1000+ | 1 | 7.722 | 0.000 | 7.722 | 7.722 | 7.722 | 7.722 |
+
+## Lift by band
+
+| Set | N | Mean | Stddev | p50 | p90 | Min | Max |
+|---|---|---|---|---|---|---|---|
+| 10000 MHz | 4444 | 6.578 | 5.300 | 5.573 | 12.879 | -0.693 | 26.814 |
+| 24000 MHz | 465 | 6.431 | 5.895 | 5.127 | 14.450 | -0.693 | 26.814 |
+| 47000 MHz | 77 | 7.232 | 6.062 | 5.804 | 16.799 | 0.355 | 25.762 |
+| 75000 MHz | 14 | 4.824 | 2.365 | 3.339 | 7.317 | 2.483 | 7.317 |
+
diff --git a/priv/backtest_reports/baseline_time_of_day.md b/priv/backtest_reports/baseline_time_of_day.md
new file mode 100644
index 00000000..57874aaf
--- /dev/null
+++ b/priv/backtest_reports/baseline_time_of_day.md
@@ -0,0 +1,31 @@
+# Backtest: Microwaveprop.Backtest.Features.time_of_day
+
+QSO sample: 5000  
+Baseline sample: 5000
+
+## Matched distribution
+
+| Set | N | Mean | Stddev | p50 | p90 | Min | Max |
+|---|---|---|---|---|---|---|---|
+| QSO times | 5000 | 16.530 | 5.702 | 17.567 | 22.317 | 0.000 | 23.983 |
+| Random baseline | 5000 | 12.071 | 6.934 | 12.017 | 21.733 | 0.000 | 23.983 |
+
+## Lift by distance (km)
+
+| Set | N | Mean | Stddev | p50 | p90 | Min | Max |
+|---|---|---|---|---|---|---|---|
+| 0-100 | 1249 | 16.326 | 5.829 | 17.633 | 22.033 | 0.000 | 23.917 |
+| 100-250 | 2295 | 17.003 | 5.704 | 18.033 | 22.583 | 0.017 | 23.983 |
+| 250-500 | 1288 | 15.870 | 5.495 | 16.650 | 21.583 | 0.000 | 23.983 |
+| 500-1000 | 167 | 16.719 | 5.667 | 18.283 | 20.333 | 0.050 | 23.933 |
+| 1000+ | 1 | 2.867 | 0.000 | 2.867 | 2.867 | 2.867 | 2.867 |
+
+## Lift by band
+
+| Set | N | Mean | Stddev | p50 | p90 | Min | Max |
+|---|---|---|---|---|---|---|---|
+| 10000 MHz | 4444 | 16.540 | 5.708 | 17.600 | 22.267 | 0.000 | 23.983 |
+| 24000 MHz | 465 | 16.774 | 5.217 | 16.933 | 22.650 | 0.017 | 23.933 |
+| 47000 MHz | 77 | 15.565 | 7.003 | 17.667 | 22.867 | 0.733 | 23.783 |
+| 75000 MHz | 14 | 10.538 | 8.066 | 15.000 | 17.950 | 1.633 | 17.967 |
+
diff --git a/test/microwaveprop/backtest/features_test.exs b/test/microwaveprop/backtest/features_test.exs
new file mode 100644
index 00000000..a67d750f
--- /dev/null
+++ b/test/microwaveprop/backtest/features_test.exs
@@ -0,0 +1,62 @@
+defmodule Microwaveprop.Backtest.FeaturesTest do
+  use Microwaveprop.DataCase, async: true
+
+  alias Microwaveprop.Backtest.Features
+  alias Microwaveprop.Weather
+
+  @point_lat 32.9
+  @point_lon -97.0
+  @valid_time ~U[2026-03-28 18:00:00Z]
+
+  defp create_profile(attrs) do
+    base = %{
+      valid_time: @valid_time,
+      lat: @point_lat,
+      lon: @point_lon
+    }
+
+    {:ok, _} = Weather.upsert_hrrr_profile(Map.merge(base, attrs))
+  end
+
+  describe "naive_gradient/3" do
+    test "returns the nearest profile's min_refractivity_gradient" do
+      create_profile(%{min_refractivity_gradient: -250.0})
+      assert Features.naive_gradient(@point_lat, @point_lon, @valid_time) == -250.0
+    end
+
+    test "returns nil when no profile is within match window" do
+      assert Features.naive_gradient(@point_lat, @point_lon, @valid_time) == nil
+    end
+  end
+
+  describe "td_depression/3" do
+    test "returns surface_temp_c minus surface_dewpoint_c" do
+      create_profile(%{surface_temp_c: 20.0, surface_dewpoint_c: 14.0})
+      assert Features.td_depression(@point_lat, @point_lon, @valid_time) == 6.0
+    end
+
+    test "returns nil when profile has no surface data" do
+      create_profile(%{surface_temp_c: nil, surface_dewpoint_c: nil})
+      assert Features.td_depression(@point_lat, @point_lon, @valid_time) == nil
+    end
+  end
+
+  describe "time_of_day/3" do
+    test "returns UTC hour + minute fraction and never calls the database" do
+      assert Features.time_of_day(0.0, 0.0, ~U[2026-01-01 06:30:00Z]) == 6.5
+      assert Features.time_of_day(0.0, 0.0, ~U[2026-01-01 00:00:00Z]) == 0.0
+      assert Features.time_of_day(0.0, 0.0, ~U[2026-01-01 23:45:00Z]) == 23.75
+    end
+  end
+
+  describe "pressure/3" do
+    test "returns surface_pressure_mb from the nearest profile" do
+      create_profile(%{surface_pressure_mb: 1013.2})
+      assert Features.pressure(@point_lat, @point_lon, @valid_time) == 1013.2
+    end
+
+    test "returns nil when no profile exists" do
+      assert Features.pressure(@point_lat, @point_lon, @valid_time) == nil
+    end
+  end
+end
diff --git a/test/microwaveprop/backtest_test.exs b/test/microwaveprop/backtest_test.exs
new file mode 100644
index 00000000..166ea09e
--- /dev/null
+++ b/test/microwaveprop/backtest_test.exs
@@ -0,0 +1,170 @@
+defmodule Microwaveprop.BacktestTest do
+  use Microwaveprop.DataCase, async: true
+
+  alias Microwaveprop.Backtest
+  alias Microwaveprop.Radio.Contact
+
+  defp create_contact(attrs \\ %{}) do
+    default = %{
+      station1: "W5XD",
+      station2: "K5TR",
+      qso_timestamp: ~U[2026-03-28 18:00:00Z],
+      mode: "CW",
+      band: Decimal.new("10000"),
+      grid1: "EM12",
+      grid2: "EM00",
+      pos1: %{"lat" => 32.9, "lon" => -97.0},
+      pos2: %{"lat" => 30.3, "lon" => -97.7},
+      distance_km: Decimal.new("295")
+    }
+
+    {:ok, contact} =
+      %Contact{}
+      |> Contact.changeset(Map.merge(default, attrs))
+      |> Repo.insert()
+
+    contact
+  end
+
+  describe "evaluate/2" do
+    test "returns counts and a qso distribution for a trivial feature" do
+      create_contact()
+
+      feature = fn _lat, _lon, _time -> 1.0 end
+      report = Backtest.evaluate(feature, sample_size: 10, baseline_size: 10)
+
+      assert report.qso_count == 1
+      assert report.baseline_count == 10
+      assert %Backtest.Distribution{} = report.qso_distribution
+      assert %Backtest.Distribution{} = report.baseline_distribution
+      assert report.qso_distribution.count == 1
+      assert report.qso_distribution.mean == 1.0
+    end
+
+    test "skips contacts with nil pos1" do
+      create_contact(%{pos1: nil, pos2: nil})
+
+      feature = fn _lat, _lon, _time -> 1.0 end
+      report = Backtest.evaluate(feature, sample_size: 10, baseline_size: 0)
+
+      assert report.qso_count == 0
+    end
+
+    test "excludes nil feature values from the distribution" do
+      create_contact(%{station1: "A"})
+      create_contact(%{station1: "B"})
+
+      feature = fn _lat, _lon, _time -> nil end
+      report = Backtest.evaluate(feature, sample_size: 10, baseline_size: 0)
+
+      assert report.qso_count == 2
+      assert report.qso_distribution.count == 0
+    end
+  end
+
+  describe "lift_by_distance/2" do
+    test "bins QSOs by distance_km and summarizes the feature per bin" do
+      create_contact(%{station1: "A", distance_km: Decimal.new("50"), pos1: %{"lat" => 1.0, "lon" => 0.0}})
+      create_contact(%{station1: "B", distance_km: Decimal.new("300"), pos1: %{"lat" => 5.0, "lon" => 0.0}})
+      create_contact(%{station1: "C", distance_km: Decimal.new("800"), pos1: %{"lat" => 9.0, "lon" => 0.0}})
+
+      f = fn lat, _lon, _time -> lat end
+      bins = Backtest.lift_by_distance(f, sample_size: 10)
+
+      assert bins["0-100"].count == 1
+      assert bins["0-100"].mean == 1.0
+      assert bins["100-250"].count == 0
+      assert bins["250-500"].count == 1
+      assert bins["250-500"].mean == 5.0
+      assert bins["500-1000"].count == 1
+      assert bins["500-1000"].mean == 9.0
+      assert bins["1000+"].count == 0
+    end
+
+    test "drops nil feature values from the bin stats" do
+      create_contact(%{station1: "A", distance_km: Decimal.new("50"), pos1: %{"lat" => 1.0, "lon" => 0.0}})
+      create_contact(%{station1: "B", distance_km: Decimal.new("60"), pos1: %{"lat" => 2.0, "lon" => 0.0}})
+
+      f = fn lat, _lon, _time -> if lat == 1.0, do: nil, else: 2.0 end
+      bins = Backtest.lift_by_distance(f, sample_size: 10)
+
+      assert bins["0-100"].count == 1
+      assert bins["0-100"].mean == 2.0
+    end
+  end
+
+  describe "random_baseline/2" do
+    test "generates N samples drawn from the contact locations" do
+      create_contact(%{
+        station1: "A",
+        qso_timestamp: ~U[2026-03-01 00:00:00Z],
+        pos1: %{"lat" => 1.0, "lon" => 10.0}
+      })
+
+      create_contact(%{
+        station1: "B",
+        qso_timestamp: ~U[2026-06-15 12:00:00Z],
+        pos1: %{"lat" => 2.0, "lon" => 20.0}
+      })
+
+      samples = Backtest.random_baseline(5, sample_size: 10)
+      assert length(samples) == 5
+
+      Enum.each(samples, fn {lat, lon, time} ->
+        assert lat in [1.0, 2.0]
+        assert lon in [10.0, 20.0]
+        assert %DateTime{} = time
+      end)
+    end
+
+    test "keeps timestamps within ±30 days of a source QSO" do
+      create_contact(%{
+        station1: "A",
+        qso_timestamp: ~U[2026-03-01 00:00:00Z],
+        pos1: %{"lat" => 1.0, "lon" => 10.0}
+      })
+
+      samples = Backtest.random_baseline(200, sample_size: 10)
+      source = ~U[2026-03-01 00:00:00Z]
+
+      Enum.each(samples, fn {_, _, time} ->
+        delta_seconds = abs(DateTime.diff(time, source))
+        assert delta_seconds <= 30 * 24 * 3600
+      end)
+    end
+
+    test "returns an empty list when there are no contacts with a position" do
+      assert Backtest.random_baseline(10) == []
+    end
+  end
+
+  describe "lift_by_band/2" do
+    test "groups feature values by band and summarizes per band" do
+      create_contact(%{
+        station1: "A",
+        band: Decimal.new("10000"),
+        pos1: %{"lat" => 1.0, "lon" => 0.0}
+      })
+
+      create_contact(%{
+        station1: "B",
+        band: Decimal.new("24000"),
+        pos1: %{"lat" => 7.0, "lon" => 0.0}
+      })
+
+      create_contact(%{
+        station1: "C",
+        band: Decimal.new("24000"),
+        pos1: %{"lat" => 9.0, "lon" => 0.0}
+      })
+
+      f = fn lat, _lon, _time -> lat end
+      bands = Backtest.lift_by_band(f, sample_size: 10)
+
+      assert bands[Decimal.new("10000")].count == 1
+      assert bands[Decimal.new("10000")].mean == 1.0
+      assert bands[Decimal.new("24000")].count == 2
+      assert bands[Decimal.new("24000")].mean == 8.0
+    end
+  end
+end