prop/lib/microwaveprop/propagation/scorer.ex

defmodule Microwaveprop.Propagation.Scorer do
  @moduledoc """
  Scoring functions for microwave propagation conditions.

  Each of the 10 factors produces a 0-100 score. The composite score is
  a weighted sum using weights from `BandConfig`. All thresholds and
  parameters are read from `BandConfig` — no hardcoded magic numbers here.
  """

  alias Microwaveprop.Propagation.BandConfig

  # ── Temperature conversion helpers ────────────────────────────────
  alias Microwaveprop.Propagation.Region

  # Compile-time inverse of the Marshall-Palmer b exponent (1/1.6).
  # Pre-computed to avoid one float division per rain pixel in the
  # dBZ → rain-rate hot path.
  @mp_inv_b 1 / 1.6

  @doc "Converts Fahrenheit to Celsius. Returns nil for nil input."
  @spec f_to_c(number() | nil) :: float() | nil
  def f_to_c(nil), do: nil
  def f_to_c(f), do: (f - 32) * 5 / 9

  @doc "Converts Celsius to Fahrenheit. Returns nil for nil input."
  @spec c_to_f(number() | nil) :: float() | nil
  def c_to_f(nil), do: nil
  def c_to_f(c), do: c * 9 / 5 + 32

  # ── Derived value helpers ─────────────────────────────────────────

  @doc """
  Computes absolute humidity in g/m3 from temperature and dewpoint (both Celsius).

  Formula: 217 * (6.112 * exp(17.67 * Td / (Td + 243.5))) / (T + 273.15)
  """
  @spec absolute_humidity(number(), number()) :: float()
  def absolute_humidity(temp_c, dewpoint_c) do
    e_sat = 6.112 * :math.exp(17.67 * dewpoint_c / (dewpoint_c + 243.5))
    217.0 * e_sat / (temp_c + 273.15)
  end

  @doc "Computes wind speed in knots from u and v components in m/s. Returns nil if either is nil."
  @spec wind_speed_kts(number() | nil, number() | nil) :: float() | nil
  def wind_speed_kts(nil, _v), do: nil
  def wind_speed_kts(_u, nil), do: nil

  def wind_speed_kts(u_ms, v_ms) do
    :math.sqrt(u_ms * u_ms + v_ms * v_ms) * 1.94384
  end

  @doc "Converts precipitation accumulation (mm) to rate (mm/hr). Returns 0.0 for nil, zero, or negative."
  @spec precip_to_rate_mmhr(number() | nil) :: float()
  def precip_to_rate_mmhr(nil), do: 0.0
  def precip_to_rate_mmhr(mm) when mm > 0, do: mm / 1
  def precip_to_rate_mmhr(_mm), do: 0.0

  @doc """
  Converts NEXRAD composite reflectivity (dBZ) to rain rate (mm/hr) via the
  Marshall-Palmer Z-R relationship `Z = 200 * R^1.6`, so `R = (Z/200)^(1/1.6)`
  with `Z = 10^(dBZ/10)`.

  Clipping rules:
  * Below 5 dBZ — not rain (ground clutter / clear air). Returns 0.0.
  * Above 150 mm/hr — hail contamination (55+ dBZ can return >100 mm/hr
    from pure Marshall-Palmer). The ITU-R P.838 rain table tops out at
    ~150 mm/hr so we clip there to stay inside the calibrated regime.
  """
  @spec dbz_to_rain_rate_mmhr(number() | nil) :: float()
  def dbz_to_rain_rate_mmhr(nil), do: 0.0
  def dbz_to_rain_rate_mmhr(dbz) when dbz < 5.0, do: 0.0

  def dbz_to_rain_rate_mmhr(dbz) do
    z = :math.pow(10, dbz / 10)
    r = :math.pow(z / 200, @mp_inv_b)
    min(r, 150.0)
  end

  # ── Factor 1: Humidity ────────────────────────────────────────────

  @doc """
  Scores absolute humidity (g/m3) for the given band.

  Beneficial bands (10 GHz): higher humidity improves ducting.
  Harmful bands (24 GHz+): humidity causes absorption loss.
  """
  @spec score_humidity(number(), map()) :: integer()
  def score_humidity(abs_humidity, %{humidity_effect: :beneficial}) do
    thresholds = BandConfig.humidity_beneficial_thresholds()

    case Enum.find(thresholds, fn {max_h, _score} -> abs_humidity < max_h end) do
      {_max, score} -> score
      nil -> BandConfig.humidity_beneficial_default()
    end
  end

  def score_humidity(abs_humidity, %{humidity_effect: :harmful, humidity_penalty: penalty}) do
    r = abs_humidity * penalty

    cond do
      r <= 6 -> 100
      r <= 9 -> round(95 - (r - 6) / 3 * 20)
      r <= 13 -> round(75 - (r - 9) / 4 * 30)
      r <= 18 -> round(45 - (r - 13) / 5 * 35)
      true -> max(0, round(10 - (r - 18) * 2))
    end
  end

  # ── Factor 2: Time of day ─────────────────────────────────────────

  @doc """
  Scores time of day for propagation conditions.

  Returns {score, label} where score is 0-100 and label describes the period.
  Uses longitude-based solar time offset (longitude / 15) so each grid point
  gets its own local time rather than a fixed timezone offset.
  """
  @spec score_time_of_day(integer(), integer(), integer(), number()) :: {integer(), String.t()}
  def score_time_of_day(utc_hour, utc_minute, month, longitude) do
    offset = longitude / 15
    local = :math.fmod(utc_hour + utc_minute / 60 + offset + 24, 24)
    sunrise = Enum.at(BandConfig.sunrise_table(), month - 1)
    d = local - sunrise

    classify_time_period(d, local)
  end

  defp classify_time_period(d, _local) when d >= -1.5 and d <= 1.5, do: {100, "Peak — inversion maximum"}

  defp classify_time_period(d, _local) when d > 1.5 and d <= 3.0, do: {78, "Good — inversion eroding"}

  defp classify_time_period(d, _local) when d > -3.0 and d < -1.5, do: {82, "Pre-dawn — inversion building"}

  defp classify_time_period(d, _local) when d > 3.0 and d <= 6.0, do: {38, "Marginal — boundary layer mixing"}

  defp classify_time_period(_d, local) when local >= 20 or local <= 1, do: {72, "Evening — cooling, inversion reforming"}

  defp classify_time_period(d, _local) when d > 6, do: {18, "Afternoon — full convective mixing"}

  defp classify_time_period(_d, _local), do: {55, "Night — gradual cooling"}

  # ── Factor 3: Temp-dewpoint depression ────────────────────────────

  @doc """
  Scores the temperature-dewpoint depression (both in Fahrenheit).

  Beneficial bands: moisture helps (tight depression = moist = better ducting,
  but too tight means fog which is bad).
  Harmful bands: dry air is better (wide depression = less absorption).
  """
  @spec score_td_depression(number(), number(), map()) :: integer()
  def score_td_depression(temp_f, dewpoint_f, %{humidity_effect: :beneficial}) do
    dep = temp_f - dewpoint_f

    cond do
      dep < 3 -> 40
      dep < 8 -> 75
      dep < 14 -> 85
      dep < 22 -> 70
      true -> 55
    end
  end

  def score_td_depression(temp_f, dewpoint_f, %{humidity_effect: :harmful}) do
    dep = temp_f - dewpoint_f

    cond do
      dep > 22 -> 96
      dep > 14 -> 80
      dep > 8 -> 60
      dep > 4 -> 38
      true -> 18
    end
  end

  # ── Factor 4: Refractivity ───────────────────────────────────────

  @doc """
  Scores minimum refractivity gradient and boundary layer depth.

  Nil gradient returns 50 (unknown). Otherwise walks thresholds from
  BandConfig and then applies an HPBL multiplier: shallow boundary layers
  amplify trapping (surface inversion → steep gradient → ducting), deep
  boundary layers indicate convective mixing that disrupts ducts.

  As of the 2026-04-25 algo revisions the HPBL boundary-layer multiplier
  is retired — `bl_depth_m` is accepted (and the schema still stores it
  for diagnostics) but it does not modify the score. See
  `docs/algo-reports/2026-04-25-algo-revisions.md` Recommendation 2: on
  the n=47,418 10 GHz HRRR-matched corpus rho_hpbl was +0.004, well below
  the 0.05 noise floor; the bin tables agreed within ±5 % across
  200–2,000 m. The previously reported "shallow-BL → longer-distance"
  effect was a small-sample artefact and disappeared once the matched
  corpus grew past ~5,000 contacts.
  """
  @spec score_refractivity(number() | nil, number() | nil, map()) :: integer()
  def score_refractivity(min_gradient, bl_depth_m, band_config) do
    score_refractivity(min_gradient, bl_depth_m, nil, band_config)
  end

  @doc """
  Scores refractivity, ignoring `best_duct_band_ghz`.

  Accepts the native-profile duct band so existing call sites
  (`Scorer.score_grid_point/2`, `Recalibrator`, the Rust comparator)
  continue to compile, but the previous 1.15× boost was retired in the
  2026-04-25 revisions: at 10 GHz on n=52,341 matched contacts, cells
  whose native duct supported the band ran *shorter* (198 km, n=173)
  than no-duct cells (211 km, n=44,658). Refractivity is now the
  gradient-only score.
  """
  @spec score_refractivity(number() | nil, number() | nil, number() | nil, map()) :: integer()
  def score_refractivity(min_gradient, bl_depth_m, best_duct_band_ghz, band_config) do
    score_refractivity(min_gradient, bl_depth_m, best_duct_band_ghz, nil, band_config)
  end

  @doc """
  Scores refractivity, ignoring `best_duct_band_ghz` and `bulk_richardson`.

  The Bulk Richardson gate existed only to suppress the native-duct
  boost during mechanically-mixed conditions. With the boost retired
  (see 4-arity doc), Richardson is irrelevant — every input collapses
  onto the gradient-only base score. Arity preserved so prod call
  sites don't need an emergency rewrite.
  """
  @spec score_refractivity(number() | nil, number() | nil, number() | nil, number() | nil, map()) ::
          integer()
  def score_refractivity(nil, _bl_depth_m, _best_duct_band_ghz, _bulk_richardson, _band_config), do: 50

  def score_refractivity(min_gradient, _bl_depth_m, _best_duct_band_ghz, _bulk_richardson, %{humidity_effect: effect}) do
    thresholds = BandConfig.refractivity_thresholds()

    case find_refractivity_threshold(min_gradient, thresholds, effect) do
      {:ok, score} ->
        score

      :none ->
        {beneficial_default, harmful_default} = BandConfig.refractivity_default()
        if effect == :beneficial, do: beneficial_default, else: harmful_default
    end
  end

  @doc """
  Inverse-sensor boost from commercial LOS link degradation.

  Nearby short-path commercial links (11/24/68 GHz in the DFW area) act as a
  physical refractivity anomaly sensor: when their rx_power drops 5+ dB below
  the 7-day baseline *without* an equipment cause, the same multipath fading
  that hurts their short Fresnel zone is what carries the beyond-LOS paths.
  This is a terminal boost applied once per composite score, not per-band —
  the underlying physics is band-agnostic at the resolution we can measure.

  Takes the incoming composite (or any 0-100) score and a degradation map
  from `Commercial.link_degradation_at/3`; returns a boosted integer score.

  * nil or empty result (`n_links == 0`) — no-op
  * <3 dB — noise floor, no-op
  * 3–8 dB — mild boost (+0 to +10)
  * ≥8 dB — strong boost (+10 to +25, clamped ≤ 100)
  """
  @spec commercial_link_boost(number(), map() | nil) :: integer()
  def commercial_link_boost(score, nil), do: clamp_score(score)

  def commercial_link_boost(score, %{n_links: 0}), do: clamp_score(score)

  def commercial_link_boost(score, %{degradation_db: db}) when db < 3.0, do: clamp_score(score)

  def commercial_link_boost(score, %{degradation_db: db}) when db < 8.0 do
    # 3 dB → +2; 7.9 dB → ~+10 (linear interpolation)
    bonus = 2 + (db - 3.0) * 8.0 / 5.0
    clamp_score(score + bonus)
  end

  def commercial_link_boost(score, %{degradation_db: db}) do
    # 8 dB → +10; 16 dB → +25 (linear)
    bonus = min(25.0, 10 + (db - 8.0) * 15.0 / 8.0)
    clamp_score(score + bonus)
  end

  defp clamp_score(value) do
    value |> round() |> max(0) |> min(100)
  end

  defp find_refractivity_threshold(gradient, thresholds, effect) do
    case Enum.find(thresholds, fn {max_grad, _ben, _harm} -> gradient < max_grad end) do
      {_max, beneficial_score, harmful_score} ->
        score =
          case effect do
            :beneficial -> beneficial_score
            :harmful -> harmful_score
          end

        {:ok, score}

      nil ->
        :none
    end
  end

  # ── Factor 5: Sky cover ──────────────────────────────────────────

  @doc "Scores sky cover percentage (0 = clear, 100 = overcast). Nil returns 50."
  @spec score_sky(number() | nil) :: integer()
  def score_sky(nil), do: 50

  def score_sky(pct) do
    cond do
      pct <= 6 -> 100
      pct <= 25 -> 88
      pct <= 50 -> 60
      pct <= 87 -> 25
      true -> 5
    end
  end

  # ── Factor 6: Season ─────────────────────────────────────────────

  @doc """
  Scores the seasonal effect for the given month and band config.

  When lat/lon are provided, applies a regional adjustment multiplier
  from `Propagation.Region` on top of the band's `seasonal_base` +
  `seasonal_adj`. This corrects for the meteorologist's observation
  that Gulf coast August is better than the uniform base suggests,
  while Corn Belt August is worse.
  """
  @spec score_season(integer(), float | nil, float | nil, map()) :: integer()
  def score_season(month, lat, lon, %{seasonal_base: base_map, seasonal_adj: adj_map}) do
    base = Map.get(base_map, month, 50)
    adj = Map.get(adj_map, month, 0)

    region_mult =
      if is_number(lat) and is_number(lon) do
        region = Region.for_point(lat, lon)
        Region.seasonal_adjustment(region, month)
      else
        1.0
      end

    round(min(100, max(0, (base + adj) * region_mult)))
  end

  # ── Factor 7: Wind ───────────────────────────────────────────────

  @doc "Scores wind speed in knots. Nil returns 50."
  @spec score_wind(number() | nil) :: integer()
  def score_wind(nil), do: 50

  def score_wind(speed_kts) do
    cond do
      speed_kts < 5 -> 100
      speed_kts < 10 -> 90
      speed_kts < 15 -> 75
      speed_kts < 20 -> 55
      speed_kts < 25 -> 35
      true -> 15
    end
  end

  # ── Factor 8: Rain attenuation ───────────────────────────────────

  @doc """
  Scores rain attenuation for the given rate (mm/hr) and band config.

  Uses ITU-R P.838 specific attenuation: gamma = k * R^alpha (dB/km).
  """
  @spec score_rain(number() | nil, map()) :: integer()
  def score_rain(nil, _band_config), do: 100
  def score_rain(rate, _band_config) when rate == 0, do: 100

  def score_rain(rain_rate_mmhr, %{rain_k: k, rain_alpha: alpha}) do
    gamma = k * :math.pow(rain_rate_mmhr, alpha)

    cond do
      gamma < 0.1 -> 95
      gamma < 0.5 -> 75
      gamma < 1.0 -> 50
      gamma < 2.0 -> 25
      gamma < 5.0 -> 10
      true -> 0
    end
  end

  # ── Factor 9: Precipitable water ──────────────────────────────────

  @doc """
  Scores precipitable water (mm).

  Beneficial bands (10 GHz): moderate PWAT is optimal for refractivity.
  Harmful bands (24+ GHz): lower PWAT is better (less water vapor absorption).
  """
  @spec score_pwat(number() | nil, map()) :: integer()
  def score_pwat(nil, _band_config), do: 60

  def score_pwat(pwat_mm, %{humidity_effect: :beneficial}) do
    cond do
      pwat_mm < 10 -> 55
      pwat_mm < 20 -> 75
      pwat_mm < 30 -> 90
      pwat_mm < 40 -> 70
      true -> 50
    end
  end

  def score_pwat(pwat_mm, %{humidity_effect: :harmful}) do
    cond do
      pwat_mm < 10 -> 95
      pwat_mm < 20 -> 80
      pwat_mm < 30 -> 60
      pwat_mm < 40 -> 35
      true -> 15
    end
  end

  # ── Factor 10: Pressure trend ────────────────────────────────────

  @doc """
  Scores barometric pressure and trend.

  Without previous reading, scores based on absolute pressure.
  With previous reading, scores based on pressure change (delta).
  """
  @spec score_pressure(number() | nil, number() | nil) :: integer()
  def score_pressure(nil, _previous_mb), do: 50

  def score_pressure(current_mb, nil) do
    cond do
      current_mb < 980 -> 88
      current_mb < 990 -> 82
      current_mb < 1000 -> 70
      current_mb < 1010 -> 55
      current_mb < 1020 -> 40
      true -> 30
    end
  end

  def score_pressure(current_mb, previous_mb) do
    delta = current_mb - previous_mb

    cond do
      delta > 2.5 -> 80
      delta > 0.8 -> 70
      delta > -0.5 -> 60
      delta > -2.0 -> 65
      true -> 45
    end
  end

  # ── Composite score ──────────────────────────────────────────────

  @doc """
  Precompute the four band-invariant factors so the grid scorer can
  reuse them across all 17 band iterations for a single point. Saves
  ~30% of the scoring loop by hoisting shared arithmetic out of the
  per-band call.
  """
  @spec precompute_band_invariants(map()) :: %{
          tod_score: integer(),
          sky_score: integer(),
          wind_score: integer(),
          pressure_score: integer()
        }
  def precompute_band_invariants(conditions) do
    {tod_score, _label} =
      score_time_of_day(conditions.utc_hour, conditions.utc_minute, conditions.month, conditions.longitude)

    %{
      tod_score: tod_score,
      sky_score: score_sky(conditions.sky_cover_pct),
      wind_score: score_wind(conditions.wind_speed_kts),
      pressure_score: score_pressure(conditions.pressure_mb, conditions.prev_pressure_mb)
    }
  end

  @doc """
  Computes the weighted composite propagation score.

  Takes a conditions map and band config, returns %{score: 0-100, factors: %{...}}.
  """
  @spec composite_score(map(), map()) :: %{score: integer(), factors: map()}
  def composite_score(conditions, band_config) do
    # Contract: band invariants are either all precomputed (via
    # `precompute_band_invariants/1`, typically merged in by the grid
    # scorer for per-point reuse across bands) or all absent. We check
    # one key and derive the rest from the same branch so we never
    # silently mix cached and freshly-computed values.
    %{
      tod_score: tod_score,
      sky_score: sky_score,
      wind_score: wind_score,
      pressure_score: pressure_score
    } =
      if Map.has_key?(conditions, :tod_score) do
        conditions
      else
        precompute_band_invariants(conditions)
      end

    factors = %{
      humidity: score_humidity(conditions.abs_humidity, band_config),
      time_of_day: tod_score,
      td_depression: score_td_depression(conditions.temp_f, conditions.dewpoint_f, band_config),
      refractivity:
        score_refractivity(
          conditions.min_refractivity_gradient,
          conditions.bl_depth_m,
          conditions[:best_duct_band_ghz],
          conditions[:bulk_richardson],
          band_config
        ),
      sky: sky_score,
      season: score_season(conditions.month, conditions[:latitude], conditions[:longitude], band_config),
      wind: wind_score,
      rain: score_rain(conditions.rain_rate_mmhr, band_config),
      pwat: score_pwat(conditions[:pwat_mm], band_config),
      pressure: pressure_score
    }

    weights = BandConfig.weights(band_config)

    weighted_sum =
      Enum.reduce(factors, 0.0, fn {factor, score}, acc ->
        acc + score * Map.fetch!(weights, factor)
      end)

    %{score: round(weighted_sum), factors: factors}
  end

  @doc """
  Merges multiple HRRR profiles along a path into a single conditions map.

  Strategy:
  - Beneficial factors (refractivity, pressure): use BEST (most favorable) along path
  - Harmful factors (rain, wind): use WORST (least favorable) along path
  - Other factors (temp, dewpoint, PWAT, BL depth): use path AVERAGE
  - Time/season/sky: taken from first profile (same for entire path)
  """
  @spec path_integrated_conditions([map()], map()) :: map() | nil
  def path_integrated_conditions(profiles, contact) do
    extracted = extract_profile_fields(profiles)
    build_path_conditions(extracted, contact)
  end

  defp extract_profile_fields(profiles) do
    Enum.reduce(profiles, {[], [], [], [], [], []}, fn p, {t, d, pr, g, b, pw} ->
      {
        if(is_nil(p.surface_temp_c), do: t, else: [p.surface_temp_c | t]),
        if(is_nil(p.surface_dewpoint_c), do: d, else: [p.surface_dewpoint_c | d]),
        if(is_nil(p.surface_pressure_mb), do: pr, else: [p.surface_pressure_mb | pr]),
        if(is_nil(p.min_refractivity_gradient), do: g, else: [p.min_refractivity_gradient | g]),
        if(is_nil(p.hpbl_m), do: b, else: [p.hpbl_m | b]),
        if(is_nil(p.pwat_mm), do: pw, else: [p.pwat_mm | pw])
      }
    end)
  end

  defp build_path_conditions({temps, _dewpoints, _pressures, _gradients, _bl_depths, _pwats}, _contact) when temps == [],
    do: nil

  defp build_path_conditions({_temps, dewpoints, _pressures, _gradients, _bl_depths, _pwats}, _contact)
       when dewpoints == [], do: nil

  defp build_path_conditions({temps, dewpoints, pressures, gradients, bl_depths, pwats}, contact) do
    lon = contact.pos1["lon"] || -97.0
    avg_temp_c = Enum.sum(temps) / length(temps)
    avg_dewpoint_c = Enum.sum(dewpoints) / length(dewpoints)

    %{
      abs_humidity: absolute_humidity(avg_temp_c, avg_dewpoint_c),
      temp_f: c_to_f(avg_temp_c),
      dewpoint_f: c_to_f(avg_dewpoint_c),
      wind_speed_kts: nil,
      sky_cover_pct: nil,
      utc_hour: contact.qso_timestamp.hour,
      utc_minute: contact.qso_timestamp.minute,
      month: contact.qso_timestamp.month,
      longitude: lon,
      pressure_mb: safe_min(pressures),
      prev_pressure_mb: nil,
      rain_rate_mmhr: 0.0,
      min_refractivity_gradient: safe_min(gradients),
      bl_depth_m: safe_avg(bl_depths),
      pwat_mm: safe_avg(pwats)
    }
  end

  defp safe_min([]), do: nil
  defp safe_min(list), do: Enum.min(list)

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