prop/lib/microwaveprop/backtest/features.ex

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.
  """

  import Ecto.Query

  alias Microwaveprop.Repo
  alias Microwaveprop.Weather
  alias Microwaveprop.Weather.HrrrNativeProfile

  @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

  @doc """
  Surface refractivity from the lowest level of the nearest native
  hybrid-sigma profile.

  Phase 1 sanity check: this should produce numbers comparable to the
  `naive_gradient` baseline but computed from native-level data. If
  the sign/magnitude look wrong, the native ingestion pipeline has a
  bug.

  Uses the ITU-R P.453-14 formula: N = 77.6*P/T + 3.73e5*e/T²
  where e (water vapor pressure) is derived from specific humidity.
  """
  @spec native_surface_refractivity(float, float, DateTime.t()) :: float | nil
  def native_surface_refractivity(lat, lon, valid_time) do
    with %HrrrNativeProfile{} = profile <- find_nearest_native(lat, lon, valid_time),
         t when is_float(t) <- profile.surface_temp_k,
         p when is_float(p) <- profile.surface_pressure_pa,
         q when is_float(q) <- profile.surface_spfh do
      # Water vapor pressure from specific humidity: e = q*P / (0.622 + 0.378*q)
      e = q * p / (0.622 + 0.378 * q)
      # N-units
      77.6 * p / (t * 100) + 3.73e5 * e / (t * t * 100)
    else
      _ -> nil
    end
  end

  @doc """
  Bulk Richardson number at the inversion top from the nearest native profile.

  Ri < 0.25 → turbulent (bad for propagation)
  Ri > 1.0 → laminar (good for propagation)

  Returns nil if no inversion detected or no native profile available.
  """
  @spec bulk_richardson(float, float, DateTime.t()) :: float | nil
  def bulk_richardson(lat, lon, valid_time) do
    with %HrrrNativeProfile{bulk_richardson: ri} when is_float(ri) <-
           find_nearest_native(lat, lon, valid_time) do
      ri
    else
      _ -> nil
    end
  end

  @doc """
  θₑ jump (K) across the inversion from the nearest native profile.

  Larger positive values indicate stronger thermodynamic decoupling.
  """
  @spec theta_e_jump(float, float, DateTime.t()) :: float | nil
  def theta_e_jump(lat, lon, valid_time) do
    with %HrrrNativeProfile{theta_e_jump_k: jump} when is_float(jump) <-
           find_nearest_native(lat, lon, valid_time) do
      jump
    else
      _ -> nil
    end
  end

  @doc """
  Wind shear magnitude (m/s) at the inversion top from the nearest native profile.
  """
  @spec shear_at_top(float, float, DateTime.t()) :: float | nil
  def shear_at_top(lat, lon, valid_time) do
    with %HrrrNativeProfile{shear_at_top_ms: shear} when is_float(shear) <-
           find_nearest_native(lat, lon, valid_time) do
      shear
    else
      _ -> nil
    end
  end

  defp find_nearest_native(lat, lon, valid_time) do
    dlat = 0.07
    dlon = 0.07
    time_start = DateTime.add(valid_time, -3600, :second)
    time_end = DateTime.add(valid_time, 3600, :second)

    HrrrNativeProfile
    |> where(
      [p],
      p.lat >= ^(lat - dlat) and p.lat <= ^(lat + dlat) and
        p.lon >= ^(lon - dlon) and p.lon <= ^(lon + dlon) and
        p.valid_time >= ^time_start and p.valid_time <= ^time_end
    )
    |> order_by([p],
      asc:
        fragment(
          "ABS(? - ?) + ABS(? - ?) + ABS(EXTRACT(EPOCH FROM ? - ?))",
          p.lat,
          ^lat,
          p.lon,
          ^lon,
          p.valid_time,
          ^valid_time
        )
    )
    |> limit(1)
    |> Repo.one()
  end
end