prop/lib/microwaveprop/pskr/calibration_sampler.ex

defmodule Microwaveprop.Pskr.CalibrationSampler do
  @moduledoc """
  Builds `Pskr.CalibrationSample` rows by joining PSKR spot density
  with HRRR atmospheric state and SWPC geomagnetic state at a given
  hour.

  Workflow:

  1. Pull every `pskr_spots_hourly` row for the target hour.
  2. Snap each row's midpoint to the propagation grid (0.125°) and
     bucket by `(band, snapped_lat, snapped_lon)`.
  3. For each bucket, aggregate the spot stats (sum count, count
     distinct paths, dedup modes, median distance).
  4. For each bucket, look up the nearest `hrrr_profiles` row to the
     midpoint at the target hour and copy the surface/refractivity
     fields onto the sample. HRRR is pulled in one window query and
     joined in-memory to keep the per-bucket cost a hashtable lookup.
  5. Look up SWPC Kp once per run.
  6. Bulk-upsert the result with `(hour_utc, band, midpoint_lat,
     midpoint_lon)` as the conflict target.

  The sampler is intentionally idempotent: re-running for the same
  hour is a no-op except for the `updated_at` bump, so the worker
  can safely overlap with manual reruns and the cron retry on
  failure.
  """

  import Ecto.Query

  alias Microwaveprop.Pskr.CalibrationSample
  alias Microwaveprop.Pskr.SpotHourly
  alias Microwaveprop.Repo
  alias Microwaveprop.SpaceWeather
  alias Microwaveprop.Weather.HrrrProfile

  require Logger

  @grid_step_deg 0.125
  # HRRR domain is ~3 km native, so a ±0.07° (~7 km) lat/lon window
  # plus ±60 min reliably finds at least one profile for any CONUS
  # midpoint at most hours.
  @hrrr_lookahead_deg 0.07
  @hrrr_lookahead_seconds 3600

  @doc """
  Build calibration samples for one hour. Returns the upsert count.
  Returns `0` when no PSKR spots were captured for that hour — for
  example, during PSKR client downtime.
  """
  @spec build_for_hour(DateTime.t()) :: non_neg_integer()
  def build_for_hour(%DateTime{} = hour_utc) do
    hour_utc = align_to_hour(hour_utc)
    spots = fetch_spots(hour_utc)

    if spots == [] do
      Logger.info("Pskr.CalibrationSampler: no spots for #{hour_utc} — skipping")
      0
    else
      hrrr_index = build_hrrr_index(hour_utc)
      kp = fetch_kp(hour_utc)

      samples =
        spots
        |> bucket_by_cell()
        |> Enum.map(&build_sample(&1, hour_utc, hrrr_index, kp))

      upsert(samples)
    end
  end

  # ── Spot fetch + cell bucketing ─────────────────────────────────

  defp fetch_spots(hour_utc) do
    Repo.all(
      from(s in SpotHourly,
        where: s.hour_utc == ^hour_utc and not is_nil(s.midpoint_lat) and not is_nil(s.midpoint_lon),
        select: %{
          band: s.band,
          midpoint_lat: s.midpoint_lat,
          midpoint_lon: s.midpoint_lon,
          spot_count: s.spot_count,
          distance_km: s.distance_km,
          modes: s.modes
        }
      )
    )
  end

  defp bucket_by_cell(spots) do
    Enum.group_by(spots, fn s -> {s.band, snap(s.midpoint_lat), snap(s.midpoint_lon)} end)
  end

  defp snap(deg), do: Float.round(deg / @grid_step_deg) * @grid_step_deg

  # ── HRRR feature lookup ─────────────────────────────────────────

  defp build_hrrr_index(hour_utc) do
    time_start = DateTime.add(hour_utc, -@hrrr_lookahead_seconds, :second)
    time_end = DateTime.add(hour_utc, @hrrr_lookahead_seconds, :second)

    Repo.all(
      from(h in HrrrProfile,
        where: h.valid_time >= ^time_start and h.valid_time <= ^time_end,
        select: %{
          lat: h.lat,
          lon: h.lon,
          surface_temp_c: h.surface_temp_c,
          surface_dewpoint_c: h.surface_dewpoint_c,
          pwat_mm: h.pwat_mm,
          surface_pressure_mb: h.surface_pressure_mb,
          min_refractivity_gradient: h.min_refractivity_gradient,
          hpbl_m: h.hpbl_m,
          ducting_detected: h.ducting_detected
        }
      )
    )
  end

  # Pick the HRRR row whose lat/lon is closest to (lat, lon). Linear
  # scan over the in-memory window — at the typical ~5k profile
  # window size and ~1k cells per hour, this is faster than 1k DB
  # round trips. Returns `nil` if no profile is in the window.
  defp nearest_hrrr(_lat, _lon, []), do: nil

  defp nearest_hrrr(lat, lon, profiles) do
    candidates =
      Enum.filter(profiles, fn p ->
        abs(p.lat - lat) <= @hrrr_lookahead_deg and abs(p.lon - lon) <= @hrrr_lookahead_deg
      end)

    case candidates do
      [] -> nil
      list -> Enum.min_by(list, fn p -> abs(p.lat - lat) + abs(p.lon - lon) end)
    end
  end

  # ── Kp lookup ───────────────────────────────────────────────────

  # Single Kp value applies cluster-wide for the hour. We prefer the
  # integer `kp_index` (the official 3-hour Kp) and fall back to the
  # truncated `estimated_kp` when SWPC hasn't published the official
  # value yet — same precedence as `Propagation.current_kp_index/0`.
  defp fetch_kp(_hour_utc) do
    case SpaceWeather.latest_kp() do
      %{kp_index: kp} when is_integer(kp) -> kp
      %{estimated_kp: kp} when is_number(kp) -> trunc(kp)
      _ -> nil
    end
  end

  # ── Per-bucket sample assembly ──────────────────────────────────

  defp build_sample({{band, lat, lon}, cell_spots}, hour_utc, hrrr_index, kp) do
    hrrr = nearest_hrrr(lat, lon, hrrr_index)
    now = DateTime.truncate(DateTime.utc_now(), :second)

    %{
      id: Ecto.UUID.generate(),
      hour_utc: hour_utc,
      band: band,
      midpoint_lat: lat,
      midpoint_lon: lon,
      spot_count: Enum.sum(Enum.map(cell_spots, & &1.spot_count)),
      distinct_paths: length(cell_spots),
      median_distance_km: median_distance(cell_spots),
      modes: dedup_modes(cell_spots),
      surface_temp_c: hrrr && hrrr.surface_temp_c,
      surface_dewpoint_c: hrrr && hrrr.surface_dewpoint_c,
      pwat_mm: hrrr && hrrr.pwat_mm,
      surface_pressure_mb: hrrr && hrrr.surface_pressure_mb,
      min_refractivity_gradient: hrrr && hrrr.min_refractivity_gradient,
      hpbl_m: hrrr && hrrr.hpbl_m,
      hrrr_ducting_detected: hrrr && hrrr.ducting_detected,
      kp_index: kp,
      inserted_at: now,
      updated_at: now
    }
  end

  defp median_distance(cell_spots) do
    cell_spots
    |> Enum.map(& &1.distance_km)
    |> Enum.reject(&is_nil/1)
    |> case do
      [] ->
        nil

      list ->
        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, else: Enum.at(sorted, mid)
    end
  end

  defp dedup_modes(cell_spots) do
    cell_spots |> Enum.flat_map(& &1.modes) |> Enum.uniq()
  end

  # ── Upsert ──────────────────────────────────────────────────────

  defp upsert([]), do: 0

  defp upsert(samples) do
    {count, _} =
      Repo.insert_all(CalibrationSample, samples,
        on_conflict:
          from(c in CalibrationSample,
            update: [
              set: [
                spot_count: fragment("EXCLUDED.spot_count"),
                distinct_paths: fragment("EXCLUDED.distinct_paths"),
                median_distance_km: fragment("EXCLUDED.median_distance_km"),
                modes: fragment("EXCLUDED.modes"),
                surface_temp_c: fragment("EXCLUDED.surface_temp_c"),
                surface_dewpoint_c: fragment("EXCLUDED.surface_dewpoint_c"),
                pwat_mm: fragment("EXCLUDED.pwat_mm"),
                surface_pressure_mb: fragment("EXCLUDED.surface_pressure_mb"),
                min_refractivity_gradient: fragment("EXCLUDED.min_refractivity_gradient"),
                hpbl_m: fragment("EXCLUDED.hpbl_m"),
                hrrr_ducting_detected: fragment("EXCLUDED.hrrr_ducting_detected"),
                kp_index: fragment("EXCLUDED.kp_index"),
                updated_at: fragment("EXCLUDED.updated_at")
              ]
            ]
          ),
        conflict_target: [:hour_utc, :band, :midpoint_lat, :midpoint_lon]
      )

    count
  end

  defp align_to_hour(%DateTime{} = dt) do
    %{dt | minute: 0, second: 0, microsecond: {0, 0}}
  end
end