From b80878056d110fe19098775bdaca5de0c3f294eb Mon Sep 17 00:00:00 2001 From: Graham McIntire Date: Sun, 19 Apr 2026 15:42:40 -0500 Subject: [PATCH] feat(grid-rs): Rust worker for HRRR f01..f18 propagation chain MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Extracts the memory-hostile HRRR fetch → decode → score pipeline for forecast hours 1–18 into a separate Rust service (`prop-grid-rs`). Elixir retains f00 with its native-duct + NEXRAD + commercial-link enrichment; Rust handles the other 18 steps per hourly chain. Hand-off via a new `grid_tasks` table (`FOR UPDATE SKIP LOCKED` claim from Rust, Postgrex `NOTIFY propagation_ready` back to Elixir). Rust writes the same on-disk score-grid format Elixir already uses, to the same `/data/scores` NFS tree. Phase 1 ships in shadow mode with PROP_SCORES_DIR=/data/scores_shadow. Rust crate layout at `rust/prop_grid_rs/` (1:1 module parity with the Elixir source it ports): - grid, region, band_config, scores_file, sounding_params - scorer: all 10 factors + composite. Matches the Elixir scorer byte-for-byte across 115 golden-fixture samples (5 scenarios × 23 bands). - decoder: wgrib2 subprocess + Fortran-record lola-binary parser - fetcher: HRRR URL derivation + idx cache (1h TTL) + 8-way parallel byte-range downloads with 429/5xx retry - pipeline: end-to-end chain step, f00 rejected at the boundary - db: sqlx grid_tasks claim/complete + propagation_ready NOTIFY - bin/worker: tokio main loop, JSON logs, SIGTERM-safe 91 Rust tests + clippy -D warnings clean. 2,158 Elixir tests green. Elixir additions: - `GridTaskEnqueuer` seeds fh=1..18 rows from `PropagationGridWorker.seed_chain/0` - `NotifyListener` LISTEN → warm `ScoreCache` → PubSub fan-out - `ShadowComparator` diffs prod vs shadow `.ntms` bodies - `mix rust.golden` task writes the Rust-side golden fixture k8s: new `deployment-grid-rs.yaml` pinned to talos5 (32 GB NUC) via `prop-grid-rs=primary` nodeSelector + `workload=grid-rs:NoSchedule` toleration, 512 Mi limit, sharing the existing NFS `/data` mount. The plan document is at `plans/vivid-hatching-quail.md` (local to my workstation); phases 2 (cutover) and 3 (talos5 concurrency tuning) follow after 72h of shadow-mode parity. --- k8s/deployment-grid-rs.yaml | 94 + k8s/kustomization.yaml | 1 + .../propagation/grid_task_enqueuer.ex | 60 + .../propagation/notify_listener.ex | 100 + .../propagation/shadow_comparator.ex | 119 + .../workers/propagation_grid_worker.ex | 10 + lib/mix/tasks/rust.golden.ex | 214 ++ .../20260419151243_create_grid_tasks.exs | 48 + priv/rust_golden/scores.bincode | Bin 0 -> 15069 bytes rust/prop_grid_rs/.gitignore | 1 + rust/prop_grid_rs/Cargo.lock | 2952 +++++++++++++++++ rust/prop_grid_rs/Cargo.toml | 36 + rust/prop_grid_rs/Dockerfile | 33 + rust/prop_grid_rs/src/band_config.rs | 483 +++ rust/prop_grid_rs/src/bin/worker.rs | 118 + rust/prop_grid_rs/src/db.rs | 210 ++ rust/prop_grid_rs/src/decoder.rs | 355 ++ rust/prop_grid_rs/src/fetcher.rs | 445 +++ rust/prop_grid_rs/src/grid.rs | 82 + rust/prop_grid_rs/src/lib.rs | 30 + rust/prop_grid_rs/src/pipeline.rs | 307 ++ rust/prop_grid_rs/src/region.rs | 105 + rust/prop_grid_rs/src/scorer.rs | 949 ++++++ rust/prop_grid_rs/src/scores_file.rs | 284 ++ rust/prop_grid_rs/src/sounding_params.rs | 125 + rust/prop_grid_rs/tests/scorer_golden.rs | 84 + rust/prop_grid_rs/tests/scores.bincode | Bin 0 -> 15069 bytes .../propagation/grid_task_enqueuer_test.exs | 53 + 28 files changed, 7298 insertions(+) create mode 100644 k8s/deployment-grid-rs.yaml create mode 100644 lib/microwaveprop/propagation/grid_task_enqueuer.ex create mode 100644 lib/microwaveprop/propagation/notify_listener.ex create mode 100644 lib/microwaveprop/propagation/shadow_comparator.ex create mode 100644 lib/mix/tasks/rust.golden.ex create mode 100644 priv/repo/migrations/20260419151243_create_grid_tasks.exs create mode 100644 priv/rust_golden/scores.bincode create mode 100644 rust/prop_grid_rs/.gitignore create mode 100644 rust/prop_grid_rs/Cargo.lock create mode 100644 rust/prop_grid_rs/Cargo.toml create mode 100644 rust/prop_grid_rs/Dockerfile create mode 100644 rust/prop_grid_rs/src/band_config.rs create mode 100644 rust/prop_grid_rs/src/bin/worker.rs create mode 100644 rust/prop_grid_rs/src/db.rs create mode 100644 rust/prop_grid_rs/src/decoder.rs create mode 100644 rust/prop_grid_rs/src/fetcher.rs create mode 100644 rust/prop_grid_rs/src/grid.rs create mode 100644 rust/prop_grid_rs/src/lib.rs create mode 100644 rust/prop_grid_rs/src/pipeline.rs create mode 100644 rust/prop_grid_rs/src/region.rs create mode 100644 rust/prop_grid_rs/src/scorer.rs create mode 100644 rust/prop_grid_rs/src/scores_file.rs create mode 100644 rust/prop_grid_rs/src/sounding_params.rs create mode 100644 rust/prop_grid_rs/tests/scorer_golden.rs create mode 100644 rust/prop_grid_rs/tests/scores.bincode create mode 100644 test/microwaveprop/propagation/grid_task_enqueuer_test.exs diff --git a/k8s/deployment-grid-rs.yaml b/k8s/deployment-grid-rs.yaml new file mode 100644 index 00000000..5f262537 --- /dev/null +++ b/k8s/deployment-grid-rs.yaml @@ -0,0 +1,94 @@ +apiVersion: apps/v1 +kind: Deployment +metadata: + name: prop-grid-rs + namespace: prop + # Phase 1 (shadow) default: write to /data/scores_shadow so Elixir keeps + # owning the live /data/scores tree. Flip PROP_SCORES_DIR to /data/scores + # at cutover after ShadowComparator shows < 0.5 mean delta for 72h. +spec: + replicas: 1 + minReadySeconds: 5 + strategy: + type: RollingUpdate + rollingUpdate: + maxSurge: 0 + maxUnavailable: 1 + selector: + matchLabels: + app: prop-grid-rs + template: + metadata: + labels: + app: prop-grid-rs + tier: grid-rs + spec: + # Pin to talos5 — the 32 GB NUC with NVMe. The worker peaks at ~500 Mi + # steady-state but benefits from local NVMe for the NFS scores-dir + # writes (faster fsync) and from never sharing the node with a hot + # Elixir pod. Label talos5 with `prop-grid-rs=primary` and taint it + # `workload=grid-rs:NoSchedule` to enforce this. + nodeSelector: + prop-grid-rs: "primary" + tolerations: + - key: workload + operator: Equal + value: grid-rs + effect: NoSchedule + imagePullSecrets: + - name: forgejo-registry + securityContext: + runAsUser: 65534 + runAsNonRoot: true + fsGroup: 65534 + seccompProfile: + type: RuntimeDefault + containers: + - name: prop-grid-rs + # Image built by a new CI pipeline from rust/prop_grid_rs/. + # Multi-arch (amd64 + arm64) — reuses the repo's existing + # buildx wiring. + image: git.mcintire.me/graham/prop-grid-rs:main + imagePullPolicy: IfNotPresent + env: + # Elixir secrets bundle carries DATABASE_URL already. Rust + # reads it directly (sqlx connects with the same URL). + - name: HRRR_BASE_URL + value: "http://skippy.w5isp.com:8080" + - name: PROP_SCORES_DIR + value: "/data/scores_shadow" + - name: RUST_LOG + value: "info" + - name: PROP_GRID_RS_PG_CONNS + value: "4" + envFrom: + - secretRef: + name: prop-secrets + # No HTTP surface. Liveness is "the process is still running"; + # let the kubelet restart us if we crash. Claiming work from + # grid_tasks is the implicit readiness signal — an idle pod is + # still healthy (the queue is simply empty). + resources: + requests: + cpu: 100m + memory: 128Mi + limits: + cpu: "4" + memory: 512Mi + securityContext: + allowPrivilegeEscalation: false + runAsNonRoot: true + runAsUser: 65534 + capabilities: + drop: + - ALL + seccompProfile: + type: RuntimeDefault + volumeMounts: + - name: data + mountPath: /data + volumes: + - name: data + nfs: + server: 10.0.15.103 + path: /data diff --git a/k8s/kustomization.yaml b/k8s/kustomization.yaml index 875d6b43..0d59a70f 100644 --- a/k8s/kustomization.yaml +++ b/k8s/kustomization.yaml @@ -8,5 +8,6 @@ resources: - rbac.yaml - deployment.yaml - deployment-backfill.yaml + - deployment-grid-rs.yaml - service.yaml - metrics-service.yaml diff --git a/lib/microwaveprop/propagation/grid_task_enqueuer.ex b/lib/microwaveprop/propagation/grid_task_enqueuer.ex new file mode 100644 index 00000000..b264923c --- /dev/null +++ b/lib/microwaveprop/propagation/grid_task_enqueuer.ex @@ -0,0 +1,60 @@ +defmodule Microwaveprop.Propagation.GridTaskEnqueuer do + @moduledoc """ + Seeds `grid_tasks` rows for the Rust `prop-grid-rs` worker. + + Called from `PropagationGridWorker.seed_chain/0` alongside the existing + Elixir f00..f18 Oban fan-out. Rust only claims rows with + `forecast_hour > 0`; Elixir still owns the f00 analysis-hour chain + because of native-duct + NEXRAD + commercial-link enrichment. + + Inserts are idempotent via the `(run_time, forecast_hour)` unique + index — re-seeding the same cycle is a no-op. + """ + + alias Microwaveprop.Repo + + require Logger + + @max_forecast_hour 18 + + @spec seed(DateTime.t()) :: {:ok, non_neg_integer()} | {:error, term()} + def seed(%DateTime{} = run_time) do + run_time = DateTime.truncate(run_time, :second) + now = DateTime.truncate(DateTime.utc_now(), :microsecond) + + rows = + for fh <- 1..@max_forecast_hour do + %{ + id: Ecto.UUID.bingenerate(), + run_time: run_time, + forecast_hour: fh, + valid_time: DateTime.add(run_time, fh * 3600, :second), + status: "queued", + attempt: 0, + claimed_at: nil, + completed_at: nil, + error: nil, + inserted_at: now, + updated_at: now + } + end + + {count, _} = + Repo.insert_all("grid_tasks", rows, + on_conflict: :nothing, + conflict_target: [:run_time, :forecast_hour] + ) + + if count > 0 do + Logger.info("GridTaskEnqueuer: seeded #{count} grid_tasks for run_time=#{run_time}") + else + Logger.info("GridTaskEnqueuer: run_time=#{run_time} already seeded") + end + + {:ok, count} + rescue + e -> + Logger.error("GridTaskEnqueuer: seed failed: #{inspect(e)}") + {:error, e} + end +end diff --git a/lib/microwaveprop/propagation/notify_listener.ex b/lib/microwaveprop/propagation/notify_listener.ex new file mode 100644 index 00000000..9be71687 --- /dev/null +++ b/lib/microwaveprop/propagation/notify_listener.ex @@ -0,0 +1,100 @@ +defmodule Microwaveprop.Propagation.NotifyListener do + @moduledoc """ + Listens for `NOTIFY propagation_ready ''` from the Rust + `prop-grid-rs` worker. When a notification arrives, warms + `ScoreCache` for every band at that `valid_time` (reading from the + shared NFS `/data/scores//.ntms` tree that Rust just + wrote) and fans out the standard `"propagation:updated"` PubSub so + LiveView clients refresh without polling. + + No-op under `testing: :inline` Oban / when the Postgrex.Notifications + supervisor isn't started — the listener is additive; shadow mode and + unit tests don't need it running. + """ + + use GenServer + + alias Microwaveprop.Propagation + alias Microwaveprop.Propagation.BandConfig + alias Microwaveprop.Propagation.ScoreCache + alias Microwaveprop.Repo + + require Logger + + @channel "propagation_ready" + @pubsub Microwaveprop.PubSub + @topic "propagation:updated" + + @spec start_link(keyword()) :: GenServer.on_start() + def start_link(opts) do + GenServer.start_link(__MODULE__, opts, name: __MODULE__) + end + + @impl true + def init(_opts) do + send(self(), :subscribe) + {:ok, %{ref: nil, pid: nil}} + end + + @impl true + def handle_info(:subscribe, state) do + case start_notifications_conn() do + {:ok, pid, ref} -> + {:noreply, %{state | pid: pid, ref: ref}} + + {:error, reason} -> + Logger.warning("NotifyListener: failed to subscribe to #{@channel}: #{inspect(reason)}") + # Retry in 30s. An unreachable DB at boot shouldn't wedge the + # supervisor tree. + Process.send_after(self(), :subscribe, 30_000) + {:noreply, state} + end + end + + @impl true + def handle_info({:notification, _pid, _ref, @channel, payload}, state) do + case DateTime.from_iso8601(payload) do + {:ok, valid_time, _} -> + handle_propagation_ready(valid_time) + + _ -> + Logger.warning("NotifyListener: malformed #{@channel} payload: #{inspect(payload)}") + end + + {:noreply, state} + end + + def handle_info(_msg, state), do: {:noreply, state} + + defp start_notifications_conn do + config = Application.fetch_env!(:microwaveprop, Repo) + + case Postgrex.Notifications.start_link(config) do + {:ok, pid} -> + case Postgrex.Notifications.listen(pid, @channel) do + {:ok, ref} -> + {:ok, pid, ref} + + other -> + {:error, other} + end + + other -> + other + end + end + + defp handle_propagation_ready(valid_time) do + Enum.each(BandConfig.all_bands(), fn band -> + Propagation.warm_cache_and_broadcast(band.freq_mhz, valid_time) + end) + + # Prune anything older than 2 hours so the cache doesn't grow + # unbounded when Rust covers many forecast hours back-to-back. + ScoreCache.prune_older_than(DateTime.add(DateTime.utc_now(), -2, :hour)) + + Phoenix.PubSub.broadcast(@pubsub, @topic, {:propagation_updated, [valid_time]}) + + Logger.info("NotifyListener: warmed ScoreCache + broadcast for valid_time=#{valid_time}") + end +end diff --git a/lib/microwaveprop/propagation/shadow_comparator.ex b/lib/microwaveprop/propagation/shadow_comparator.ex new file mode 100644 index 00000000..39ee2d0d --- /dev/null +++ b/lib/microwaveprop/propagation/shadow_comparator.ex @@ -0,0 +1,119 @@ +defmodule Microwaveprop.Propagation.ShadowComparator do + @moduledoc """ + Shadow-mode validator for the Rust `prop-grid-rs` port. + + During Phase 1 shadow mode the Rust worker writes score-grid files to + `/data/scores_shadow//.ntms` while Elixir keeps writing to + `/data/scores/`. This module compares the two files and logs the + per-band delta distribution so drift is visible before cutover. + + Callable from IEx on demand or from a one-shot Mix task; not started + in the supervision tree (this is temporary instrumentation). + """ + + alias Microwaveprop.Propagation.ScoresFile + + require Logger + + @default_prod_dir "/data/scores" + @default_shadow_dir "/data/scores_shadow" + + @type diff :: %{ + band_mhz: non_neg_integer(), + valid_time: DateTime.t(), + total_cells: non_neg_integer(), + matching_cells: non_neg_integer(), + missing_shadow_cells: non_neg_integer(), + max_abs_delta: non_neg_integer(), + mean_abs_delta: float() + } + + @spec compare(non_neg_integer(), DateTime.t(), keyword()) :: + {:ok, diff()} | {:error, term()} + def compare(band_mhz, %DateTime{} = valid_time, opts \\ []) do + prod_dir = Keyword.get(opts, :prod_dir, @default_prod_dir) + shadow_dir = Keyword.get(opts, :shadow_dir, @default_shadow_dir) + + with {:ok, prod} <- read_file(prod_dir, band_mhz, valid_time), + {:ok, shadow} <- read_file(shadow_dir, band_mhz, valid_time) do + {:ok, diff_payloads(band_mhz, valid_time, prod, shadow)} + end + end + + defp read_file(dir, band_mhz, valid_time) do + path = band_mhz |> ScoresFile.path_for(valid_time) |> swap_base(dir) + + case File.read(path) do + {:ok, binary} -> decode_scores_file(binary) + {:error, reason} -> {:error, {:read, path, reason}} + end + end + + defp swap_base(path, new_base) do + suffix = Path.relative_to(path, ScoresFile.base_dir()) + Path.join(new_base, suffix) + end + + # Inline parse of the v1 scores-file on-disk format so this module + # doesn't reach into `ScoresFile`'s private `decode/1`. The format is + # documented in `ScoresFile`'s moduledoc; only the body size and + # dimensions matter for a diff, so we intentionally keep this minimal. + # The literal "NTMS" magic is the file-header marker already on disk. + defp decode_scores_file( + <<"NTMS", 1::unsigned-8, _band_mhz::unsigned-little-32, _valid_time_unix::signed-little-64, + _lat_min::float-little-32, _lon_min::float-little-32, _step::float-little-32, n_rows::unsigned-little-16, + n_cols::unsigned-little-16, body::binary>> + ) do + expected = n_rows * n_cols + + if byte_size(body) == expected do + {:ok, %{n_rows: n_rows, n_cols: n_cols, scores: body}} + else + {:error, :invalid_format} + end + end + + defp decode_scores_file(_), do: {:error, :invalid_format} + + defp diff_payloads(band_mhz, valid_time, prod, shadow) do + %{scores: p_body, n_cols: n_cols, n_rows: n_rows} = prod + %{scores: s_body} = shadow + + {total, matching, missing, max_abs, sum_abs} = + for row <- 0..(n_rows - 1), + col <- 0..(n_cols - 1), + reduce: {0, 0, 0, 0, 0} do + {t, m, miss, max_d, sum_d} -> + idx = row * n_cols + col + p = :binary.at(p_body, idx) + s = :binary.at(s_body, idx) + + cond do + p == 255 and s == 255 -> + {t, m, miss, max_d, sum_d} + + p != 255 and s == 255 -> + {t + 1, m, miss + 1, max_d, sum_d} + + p == s -> + {t + 1, m + 1, miss, max_d, sum_d} + + true -> + d = abs(p - s) + {t + 1, m, miss, max(max_d, d), sum_d + d} + end + end + + mean = if total > 0, do: sum_abs / total, else: 0.0 + + %{ + band_mhz: band_mhz, + valid_time: valid_time, + total_cells: total, + matching_cells: matching, + missing_shadow_cells: missing, + max_abs_delta: max_abs, + mean_abs_delta: mean + } + end +end diff --git a/lib/microwaveprop/workers/propagation_grid_worker.ex b/lib/microwaveprop/workers/propagation_grid_worker.ex index 943b75ee..0e318698 100644 --- a/lib/microwaveprop/workers/propagation_grid_worker.ex +++ b/lib/microwaveprop/workers/propagation_grid_worker.ex @@ -45,6 +45,7 @@ defmodule Microwaveprop.Workers.PropagationGridWorker do alias Microwaveprop.Propagation alias Microwaveprop.Propagation.BandConfig alias Microwaveprop.Propagation.Grid + alias Microwaveprop.Propagation.GridTaskEnqueuer alias Microwaveprop.Propagation.ProfilesFile alias Microwaveprop.Propagation.ScoreCache alias Microwaveprop.Weather @@ -96,6 +97,15 @@ defmodule Microwaveprop.Workers.PropagationGridWorker do end Oban.insert_all(jobs) + + # Seed the grid_tasks table in parallel so the Rust `prop-grid-rs` + # worker can pick up f01..f18 from a shared NFS mount. During Phase + # 1 shadow mode this runs alongside the Elixir chain; the Rust + # writes land in /data/scores_shadow and are diffed by + # `ShadowComparator`. At cutover the Elixir fan-out above is + # trimmed to just fh=0. + _ = GridTaskEnqueuer.seed(run_time) + :ok end diff --git a/lib/mix/tasks/rust.golden.ex b/lib/mix/tasks/rust.golden.ex new file mode 100644 index 00000000..05c59b57 --- /dev/null +++ b/lib/mix/tasks/rust.golden.ex @@ -0,0 +1,214 @@ +defmodule Mix.Tasks.Rust.Golden do + @shortdoc "Generate golden test fixtures for the Rust scorer port" + @moduledoc """ + Dump Elixir-scorer output for a sample of points to + `priv/rust_golden/scores.bincode`. The Rust crate loads this fixture + in its scorer tests to assert byte-for-byte parity with the Elixir + implementation. + + The fixture format is a little-endian binary: + + n_samples : u32 + per-sample: + abs_humidity : f64 + temp_f : f64 + dewpoint_f : f64 + wind_speed_kts : f64 (NaN = None) + sky_cover_pct : f64 (NaN = None) + utc_hour : u8 + utc_minute : u8 + month : u8 + longitude : f64 + latitude : f64 (NaN = None) + pressure_mb : f64 (NaN = None) + prev_pressure_mb : f64 (NaN = None) + rain_rate_mmhr : f64 (NaN = None) + min_refractivity_gradient : f64 (NaN = None) + bl_depth_m : f64 (NaN = None) + pwat_mm : f64 (NaN = None) + best_duct_band_ghz: f64 (NaN = None) + bulk_richardson : f64 (NaN = None) + band_mhz : u32 + expected_score : i32 + + No Rust dependency on bincode's stable format — we write fixed-width + fields in a fixed order and the Rust loader reads them in the same + order. Matches the Rust `Conditions` struct field-for-field. + + mix rust.golden # writes priv/rust_golden/scores.bincode + mix rust.golden --print 5 # pretty-print first 5 samples + + """ + use Mix.Task + + alias Microwaveprop.Propagation.BandConfig + alias Microwaveprop.Propagation.Scorer + + @impl true + def run(args) do + {opts, _, _} = OptionParser.parse(args, strict: [print: :integer]) + samples = build_samples() + + path = Path.join(["priv", "rust_golden", "scores.bincode"]) + File.mkdir_p!(Path.dirname(path)) + File.write!(path, encode(samples)) + Mix.shell().info("wrote #{length(samples)} samples to #{path}") + + if n = opts[:print] do + samples |> Enum.take(n) |> Enum.each(&IO.inspect/1) + end + end + + defp build_samples do + bands = BandConfig.all_bands() + + conditions_set = [ + summer_peak(), + winter_dry_cold(), + stormy_afternoon(), + dawn_duct(), + humid_gulf() + ] + + for {c, idx} <- Enum.with_index(conditions_set), band <- bands do + result = Scorer.composite_score(c, band) + {idx, c, band.freq_mhz, result.score} + end + end + + defp summer_peak do + %{ + abs_humidity: 10.0, + temp_f: 80, + dewpoint_f: 70, + wind_speed_kts: 5, + sky_cover_pct: 20, + utc_hour: 11, + utc_minute: 15, + month: 6, + longitude: -75.0, + pressure_mb: 1020, + prev_pressure_mb: 1018, + rain_rate_mmhr: 0, + min_refractivity_gradient: -350, + bl_depth_m: 400, + pwat_mm: 25.0 + } + end + + defp winter_dry_cold do + %{ + abs_humidity: 2.5, + temp_f: 28, + dewpoint_f: 14, + wind_speed_kts: 18, + sky_cover_pct: 10, + utc_hour: 14, + utc_minute: 0, + month: 1, + longitude: -97.0, + pressure_mb: 1028, + prev_pressure_mb: 1026, + rain_rate_mmhr: 0, + min_refractivity_gradient: -60, + bl_depth_m: 1800, + pwat_mm: 5.0 + } + end + + defp stormy_afternoon do + %{ + abs_humidity: 18.0, + temp_f: 85, + dewpoint_f: 74, + wind_speed_kts: 22, + sky_cover_pct: 95, + utc_hour: 21, + utc_minute: 0, + month: 7, + longitude: -97.0, + pressure_mb: 1002, + prev_pressure_mb: 1008, + rain_rate_mmhr: 8.0, + min_refractivity_gradient: -80, + bl_depth_m: 2200, + pwat_mm: 38.0 + } + end + + defp dawn_duct do + %{ + abs_humidity: 14.0, + temp_f: 72, + dewpoint_f: 70, + wind_speed_kts: 3, + sky_cover_pct: 5, + utc_hour: 11, + utc_minute: 30, + month: 5, + longitude: -95.0, + pressure_mb: 1018, + prev_pressure_mb: 1017, + rain_rate_mmhr: 0, + min_refractivity_gradient: -220, + bl_depth_m: 140, + pwat_mm: 30.0 + } + end + + defp humid_gulf do + %{ + abs_humidity: 22.0, + temp_f: 88, + dewpoint_f: 77, + wind_speed_kts: 8, + sky_cover_pct: 55, + utc_hour: 18, + utc_minute: 0, + month: 8, + longitude: -93.0, + latitude: 29.5, + pressure_mb: 1014, + prev_pressure_mb: 1014, + rain_rate_mmhr: 0, + min_refractivity_gradient: -110, + bl_depth_m: 700, + pwat_mm: 52.0 + } + end + + defp encode(samples) do + header = <> + body = for {_idx, c, band_mhz, score} <- samples, into: <<>>, do: encode_sample(c, band_mhz, score) + header <> body + end + + defp encode_sample(c, band_mhz, score) do + f(c.abs_humidity) <> + f(c.temp_f) <> + f(c.dewpoint_f) <> + opt(c[:wind_speed_kts]) <> + opt(c[:sky_cover_pct]) <> + <> <> + f(c.longitude) <> + opt(c[:latitude]) <> + opt(c[:pressure_mb]) <> + opt(c[:prev_pressure_mb]) <> + opt(c[:rain_rate_mmhr]) <> + opt(c[:min_refractivity_gradient]) <> + opt(c[:bl_depth_m]) <> + opt(c[:pwat_mm]) <> + opt(c[:best_duct_band_ghz]) <> + opt(c[:bulk_richardson]) <> + <> + end + + # IEEE 754 quiet-NaN bit pattern, little-endian. Used as the sentinel + # for `None` on the Rust side — `f64::is_nan` → `Option::None`. + @nan_bytes <<0x7FF8_0000_0000_0000::unsigned-little-64>> + + defp f(n) when is_float(n), do: <> + defp f(n) when is_integer(n), do: <> + defp opt(nil), do: @nan_bytes + defp opt(v), do: f(v) +end diff --git a/priv/repo/migrations/20260419151243_create_grid_tasks.exs b/priv/repo/migrations/20260419151243_create_grid_tasks.exs new file mode 100644 index 00000000..d166e630 --- /dev/null +++ b/priv/repo/migrations/20260419151243_create_grid_tasks.exs @@ -0,0 +1,48 @@ +defmodule Microwaveprop.Repo.Migrations.CreateGridTasks do + use Ecto.Migration + + def change do + # Hand-off queue between Elixir's PropagationGridWorker (seed) and the + # external Rust `prop-grid-rs` worker. Each row is one forecast hour's + # worth of fetch → decode → score → write-ntms work. Rust claims with + # `SELECT ... FOR UPDATE SKIP LOCKED` on (status='queued'), marks + # 'running' while working, and transitions to 'done' / 'failed' on + # completion. On 'done' it NOTIFYs `propagation_ready` so Elixir pods + # can warm ScoreCache without polling. + create table(:grid_tasks, primary_key: false) do + add :id, :binary_id, primary_key: true + + # HRRR cycle + step. (run_time, forecast_hour) is the natural key; an + # hourly cron re-enqueueing the same step is a no-op via the unique + # index rather than a conflict that has to be handled. + add :run_time, :utc_datetime, null: false + add :forecast_hour, :integer, null: false + add :valid_time, :utc_datetime, null: false + + # 'queued' → 'running' → 'done' | 'failed'. Not an enum so we can add + # more states later without a migration. + add :status, :string, null: false, default: "queued" + + # How many times a worker has claimed this row. Incremented on each + # claim so runaway retries are observable without scanning logs. + add :attempt, :integer, null: false, default: 0 + + add :claimed_at, :utc_datetime_usec + add :completed_at, :utc_datetime_usec + add :error, :text + + timestamps(type: :utc_datetime) + end + + # Uniqueness on the natural key lets the seeder run `INSERT ... ON CONFLICT + # DO NOTHING` and stay idempotent across retries of the seed job. + create unique_index(:grid_tasks, [:run_time, :forecast_hour]) + + # Rust's claim query filters on status first. A partial index on just the + # claimable subset keeps the common lookup cheap as finished rows + # accumulate — the table stays small (19 rows per hour, pruned on a cron) + # but the principle is cheap either way. + create index(:grid_tasks, [:status], where: "status = 'queued'", name: :grid_tasks_queued_idx) + create index(:grid_tasks, [:completed_at]) + end +end diff --git a/priv/rust_golden/scores.bincode b/priv/rust_golden/scores.bincode new file mode 100644 index 0000000000000000000000000000000000000000..1a8be7f1bb9f6e95246f7bffdbb66d8b7c72cffa GIT binary patch literal 15069 zcmc(lKWh|07>D0WV!RO4L)5}Q@rpx11mht=A&H~$OiUnNE|;7dEG&e8jUa+nZlkSN zT3F}@5YpLhX_*folGbmK()rEII-XzP*j?z_xbzVzF{Xe`aOa$^nd(>Mfk-d z?A6R3vsL&ojw3XV;A6q2@ENgk46j0Q21jUI#u3V=X3OlP{Pf*yjYaq}RkQmn!eOgs z4_JiZ&zh~W2rs9{=;63`WsDvRyK`gKVLM{?+qgLv!Jbdh2>ZV$>^574o0leOgg2dk z5hl{Gm``~%Wjkfd9XD@!eshn?kd?$YhAy1fQ^rQWUoXR9^Ro+S)p}r5;_t>rzef-x z@k=a%Bz~Snki^3XlK2I-3X*uZ3X*uZ3X=Fmwj(6*jzy5f!_Ngt{1#gUNqmn*ki^3X zl6V+F5)UIt;$Z|yJd7ZTrxCowciBlPiGRoN8=^`NqjfGfiJ2z@y&E=zYN#Ueof~~DMuy#PPj-r1Yzj^7)Fr9!w8c2yqZEf zLJ|*GK@#tLXT%Muf+Qa92uXa_SELG(_l2 zcwhEGXcq~Zb;VMYt^SdTwQcB|CDoEntu^@@h%UYx(B=LDQg&;`c9s696 z#OIYYQUytTURfgul6bfyB=PwZ3Zx2>cpAY=`~o{8B=InUBpybP#KQ=Z_>7R=xZ8G| zr*1RN*JjLO&7HJ*7&&fGhAb9h(_WJIDnp~+6+2}(Jb3Ry+G@&CiC>A0evcqX;$Z|y ze21-qBpybP#KQ=Zco;zvzs!yWNqqhY8JU!lco;zv4?h src/bin/worker.rs \ + && echo '' > src/lib.rs \ + && cargo build --release --bin worker \ + && rm -rf src + +COPY src ./src +RUN cargo build --release --bin worker + +# Runtime image: slim debian plus wgrib2. libexpat/libjasper pulled in +# transitively by wgrib2 for PNG/JPEG2000 GRIB2 decompression. +FROM debian:trixie-slim AS runtime +RUN apt-get update && apt-get install -y --no-install-recommends \ + wgrib2 \ + ca-certificates \ + tini \ + && rm -rf /var/lib/apt/lists/* + +COPY --from=builder /src/target/release/worker /usr/local/bin/prop-grid-rs + +USER 65534:65534 +ENTRYPOINT ["/usr/bin/tini", "--", "/usr/local/bin/prop-grid-rs"] diff --git a/rust/prop_grid_rs/src/band_config.rs b/rust/prop_grid_rs/src/band_config.rs new file mode 100644 index 00000000..928dcff7 --- /dev/null +++ b/rust/prop_grid_rs/src/band_config.rs @@ -0,0 +1,483 @@ +//! Band configuration tables. 1:1 port of +//! `lib/microwaveprop/propagation/band_config.ex`. +//! +//! Keep the data layout identical to the Elixir module; any numerical +//! drift here will show up immediately in golden-fixture scorer tests. + +use std::sync::OnceLock; + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum HumidityEffect { + Beneficial, + Harmful, +} + +#[derive(Debug, Clone, Copy)] +pub struct Weights { + pub humidity: f64, + pub time_of_day: f64, + pub td_depression: f64, + pub refractivity: f64, + pub sky: f64, + pub season: f64, + pub wind: f64, + pub rain: f64, + pub pressure: f64, + pub pwat: f64, +} + +impl Weights { + #[allow(clippy::too_many_arguments)] + pub const fn new( + humidity: f64, + time_of_day: f64, + td_depression: f64, + refractivity: f64, + sky: f64, + season: f64, + wind: f64, + rain: f64, + pressure: f64, + pwat: f64, + ) -> Self { + Self { + humidity, + time_of_day, + td_depression, + refractivity, + sky, + season, + wind, + rain, + pressure, + pwat, + } + } +} + +/// Default weight vector fit by gradient descent on 5,000 QSOs +/// (2026-04-11). Used as fallback for bands that don't carry a `weights` +/// override. +pub const DEFAULT_WEIGHTS: Weights = Weights::new( + 0.1243, // humidity + 0.0496, // time_of_day + 0.0978, // td_depression + 0.1049, // refractivity + 0.08, // sky + 0.1112, // season + 0.08, // wind + 0.1362, // rain + 0.1032, // pressure + 0.1128, // pwat +); + +/// Sunrise hour by month (1..=12) in local solar time. Index = month - 1. +pub const SUNRISE_TABLE: [f64; 12] = [ + 7.4, 7.3, 7.0, 6.7, 6.35, 6.25, 6.35, 6.65, 6.9, 7.1, 7.35, 7.45, +]; + +pub const HUMIDITY_BENEFICIAL_THRESHOLDS: &[(u8, i32)] = &[ + (4, 55), + (7, 70), + (10, 82), + (14, 90), + (18, 95), + (22, 88), +]; +pub const HUMIDITY_BENEFICIAL_DEFAULT: i32 = 75; + +/// `(gradient_max, score_beneficial, score_harmful)`. First entry whose +/// `gradient_max` exceeds the observed gradient wins. Calibrated for +/// HRRR-derived gradients. +pub const REFRACTIVITY_THRESHOLDS: &[(i32, i32, i32)] = &[ + (-200, 98, 85), + (-150, 92, 80), + (-100, 82, 72), + (-75, 68, 62), + (-55, 55, 55), + (-40, 48, 48), +]; +pub const REFRACTIVITY_DEFAULT: (i32, i32) = (42, 42); + +#[derive(Debug, Clone)] +pub struct BandConfig { + pub freq_mhz: u32, + pub label: &'static str, + pub o2_db_km: f64, + pub h2o_coeff: f64, + pub humidity_effect: HumidityEffect, + pub humidity_penalty: f64, + pub rain_k: f64, + pub rain_alpha: f64, + /// Per-month base score for seasonal factor, index = month - 1. + pub seasonal_base: [i32; 12], + /// Per-month additive adjustment, index = month - 1. + pub seasonal_adj: [i32; 12], + pub typical_range_km: u32, + pub extended_range_km: u32, + pub exceptional_range_km: u32, + /// `None` falls back to `DEFAULT_WEIGHTS`. + pub weights: Option, +} + +impl BandConfig { + pub fn weights(&self) -> Weights { + self.weights.unwrap_or(DEFAULT_WEIGHTS) + } +} + +/// Fixed-size literal helper; month values given as (month, score) pairs +/// in Elixir map form. Panics at construction if a month is out of range. +const fn months(v: [(u8, i32); 12]) -> [i32; 12] { + // Build by iterating; can't use const fn loop with indirect writes + // cleanly, so use a match per position. + let mut out = [0; 12]; + let mut i = 0; + while i < 12 { + let (m, s) = v[i]; + // Month values always 1..=12, stored at index month - 1 is the + // convention in Elixir — we rely on the literal being in month + // order here (it is in every band config). + assert!(m >= 1 && m <= 12); + out[(m - 1) as usize] = s; + i += 1; + } + out +} + +const fn adj_from(values: &[(u8, i32)]) -> [i32; 12] { + let mut out = [0; 12]; + let mut i = 0; + while i < values.len() { + let (m, s) = values[i]; + assert!(m >= 1 && m <= 12); + out[(m - 1) as usize] = s; + i += 1; + } + out +} + +// The seasonal_base tables are identical across many bands; pull out the +// common shape once so the 22 band definitions stay scannable. +const BASE_LOW_BAND: [i32; 12] = months([ + (1, 38), (2, 40), (3, 22), (4, 55), (5, 68), (6, 90), + (7, 95), (8, 75), (9, 78), (10, 88), (11, 78), (12, 25), +]); + +const BASE_50: [i32; 12] = months([ + (1, 30), (2, 32), (3, 25), (4, 50), (5, 70), (6, 95), + (7, 95), (8, 70), (9, 65), (10, 70), (11, 60), (12, 25), +]); + +const BASE_24G: [i32; 12] = months([ + (1, 88), (2, 84), (3, 68), (4, 62), (5, 51), (6, 34), + (7, 18), (8, 18), (9, 48), (10, 68), (11, 96), (12, 88), +]); + +const BASE_47G: [i32; 12] = months([ + (1, 90), (2, 88), (3, 78), (4, 68), (5, 55), (6, 38), + (7, 22), (8, 22), (9, 48), (10, 74), (11, 96), (12, 90), +]); + +const BASE_68G: [i32; 12] = months([ + (1, 90), (2, 88), (3, 78), (4, 65), (5, 50), (6, 32), + (7, 18), (8, 18), (9, 44), (10, 70), (11, 92), (12, 90), +]); + +const BASE_75G: [i32; 12] = months([ + (1, 90), (2, 90), (3, 80), (4, 68), (5, 55), (6, 38), + (7, 22), (8, 22), (9, 48), (10, 74), (11, 96), (12, 90), +]); + +const BASE_122G: [i32; 12] = months([ + (1, 92), (2, 90), (3, 78), (4, 62), (5, 45), (6, 28), + (7, 15), (8, 15), (9, 38), (10, 68), (11, 92), (12, 92), +]); + +const BASE_134G: [i32; 12] = months([ + (1, 92), (2, 90), (3, 78), (4, 65), (5, 48), (6, 30), + (7, 18), (8, 18), (9, 42), (10, 70), (11, 92), (12, 92), +]); + +const BASE_142G: [i32; 12] = months([ + (1, 92), (2, 90), (3, 78), (4, 65), (5, 47), (6, 28), + (7, 16), (8, 16), (9, 40), (10, 68), (11, 92), (12, 92), +]); + +const BASE_145G: [i32; 12] = months([ + (1, 92), (2, 90), (3, 78), (4, 64), (5, 46), (6, 27), + (7, 15), (8, 15), (9, 38), (10, 66), (11, 92), (12, 92), +]); + +const BASE_241G: [i32; 12] = months([ + (1, 95), (2, 92), (3, 75), (4, 55), (5, 35), (6, 15), + (7, 8), (8, 8), (9, 30), (10, 65), (11, 95), (12, 95), +]); + +const BASE_288G: [i32; 12] = months([ + (1, 95), (2, 92), (3, 75), (4, 55), (5, 35), (6, 14), + (7, 7), (8, 7), (9, 28), (10, 64), (11, 95), (12, 95), +]); + +const BASE_322G: [i32; 12] = months([ + (1, 96), (2, 92), (3, 74), (4, 52), (5, 32), (6, 12), + (7, 6), (8, 6), (9, 26), (10, 62), (11, 96), (12, 96), +]); + +const BASE_403G: [i32; 12] = months([ + (1, 96), (2, 92), (3, 74), (4, 52), (5, 32), (6, 12), + (7, 6), (8, 6), (9, 26), (10, 62), (11, 96), (12, 96), +]); + +const ADJ_NONE: [i32; 12] = [0; 12]; +const ADJ_24G: [i32; 12] = adj_from(&[(5, -4), (6, -8), (7, -10), (8, -10), (9, -4)]); + +fn bands() -> &'static [BandConfig] { + static CELL: OnceLock> = OnceLock::new(); + CELL.get_or_init(build_bands) +} + +#[allow(clippy::too_many_lines)] +fn build_bands() -> Vec { + use HumidityEffect::*; + vec![ + BandConfig { + freq_mhz: 50, label: "50 MHz", + o2_db_km: 0.0, h2o_coeff: 0.0, humidity_effect: Beneficial, humidity_penalty: 0.0, + rain_k: 0.0, rain_alpha: 1.0, + seasonal_base: BASE_50, seasonal_adj: ADJ_NONE, + typical_range_km: 400, extended_range_km: 1500, exceptional_range_km: 4000, + weights: None, + }, + BandConfig { + freq_mhz: 144, label: "144 MHz", + o2_db_km: 0.0, h2o_coeff: 0.0, humidity_effect: Beneficial, humidity_penalty: 0.0, + rain_k: 0.0, rain_alpha: 1.0, + seasonal_base: BASE_LOW_BAND, seasonal_adj: ADJ_NONE, + typical_range_km: 350, extended_range_km: 900, exceptional_range_km: 2500, + weights: None, + }, + BandConfig { + freq_mhz: 222, label: "222 MHz", + o2_db_km: 0.0, h2o_coeff: 0.0, humidity_effect: Beneficial, humidity_penalty: 0.0, + rain_k: 0.0, rain_alpha: 1.0, + seasonal_base: BASE_LOW_BAND, seasonal_adj: ADJ_NONE, + typical_range_km: 320, extended_range_km: 800, exceptional_range_km: 2000, + weights: Some(Weights::new(0.1593, 0.0350, 0.1250, 0.1401, 0.0706, 0.1276, 0.0706, 0.0120, 0.0681, 0.1916)), + }, + BandConfig { + freq_mhz: 432, label: "432 MHz", + o2_db_km: 0.0, h2o_coeff: 0.0, humidity_effect: Beneficial, humidity_penalty: 0.0, + rain_k: 0.0, rain_alpha: 1.0, + seasonal_base: BASE_LOW_BAND, seasonal_adj: ADJ_NONE, + typical_range_km: 300, extended_range_km: 700, exceptional_range_km: 1800, + weights: Some(Weights::new(0.2061, 0.0329, 0.1200, 0.1186, 0.0663, 0.1106, 0.0663, 0.0113, 0.0809, 0.1870)), + }, + BandConfig { + freq_mhz: 902, label: "902 MHz", + o2_db_km: 0.006, h2o_coeff: 0.0, humidity_effect: Beneficial, humidity_penalty: 0.0, + rain_k: 0.000, rain_alpha: 1.0, + seasonal_base: BASE_LOW_BAND, seasonal_adj: ADJ_NONE, + typical_range_km: 400, extended_range_km: 800, exceptional_range_km: 1500, + weights: Some(Weights::new(0.2201, 0.0437, 0.1102, 0.0888, 0.0783, 0.1197, 0.0783, 0.0133, 0.0839, 0.1635)), + }, + BandConfig { + freq_mhz: 1_296, label: "1296 MHz", + o2_db_km: 0.006, h2o_coeff: 0.0, humidity_effect: Beneficial, humidity_penalty: 0.0, + rain_k: 0.000, rain_alpha: 1.0, + seasonal_base: BASE_LOW_BAND, seasonal_adj: ADJ_NONE, + typical_range_km: 350, extended_range_km: 700, exceptional_range_km: 1200, + weights: Some(Weights::new(0.2131, 0.0494, 0.0898, 0.1392, 0.0797, 0.1108, 0.0797, 0.0136, 0.0568, 0.1678)), + }, + BandConfig { + freq_mhz: 2_304, label: "2304 MHz", + o2_db_km: 0.006, h2o_coeff: 0.0, humidity_effect: Beneficial, humidity_penalty: 0.0, + rain_k: 0.001, rain_alpha: 1.15, + seasonal_base: BASE_LOW_BAND, seasonal_adj: ADJ_NONE, + typical_range_km: 300, extended_range_km: 600, exceptional_range_km: 1000, + weights: Some(Weights::new(0.1941, 0.0387, 0.1385, 0.1638, 0.0625, 0.0868, 0.0625, 0.0336, 0.0432, 0.1762)), + }, + BandConfig { + freq_mhz: 3_400, label: "3400 MHz", + o2_db_km: 0.006, h2o_coeff: 0.0, humidity_effect: Beneficial, humidity_penalty: 0.0, + rain_k: 0.002, rain_alpha: 1.20, + seasonal_base: BASE_LOW_BAND, seasonal_adj: ADJ_NONE, + typical_range_km: 250, extended_range_km: 550, exceptional_range_km: 900, + weights: Some(Weights::new(0.1996, 0.0398, 0.1251, 0.1310, 0.0642, 0.0893, 0.0642, 0.0489, 0.0568, 0.1811)), + }, + BandConfig { + freq_mhz: 5_760, label: "5760 MHz", + o2_db_km: 0.007, h2o_coeff: 0.0, humidity_effect: Beneficial, humidity_penalty: 0.0, + rain_k: 0.005, rain_alpha: 1.25, + seasonal_base: BASE_LOW_BAND, seasonal_adj: ADJ_NONE, + typical_range_km: 220, extended_range_km: 500, exceptional_range_km: 1000, + weights: Some(Weights::new(0.1829, 0.0423, 0.1205, 0.0881, 0.0683, 0.0949, 0.0683, 0.0822, 0.0602, 0.1925)), + }, + BandConfig { + freq_mhz: 10_000, label: "10 GHz", + o2_db_km: 0.007, h2o_coeff: 0.0, humidity_effect: Beneficial, humidity_penalty: 0.0, + rain_k: 0.010, rain_alpha: 1.28, + seasonal_base: BASE_LOW_BAND, seasonal_adj: ADJ_NONE, + typical_range_km: 200, extended_range_km: 500, exceptional_range_km: 1000, + weights: None, + }, + BandConfig { + freq_mhz: 24_000, label: "24 GHz", + o2_db_km: 0.02, h2o_coeff: 0.002, humidity_effect: Harmful, humidity_penalty: 1.6, + rain_k: 0.070, rain_alpha: 1.07, + seasonal_base: BASE_24G, seasonal_adj: ADJ_24G, + typical_range_km: 100, extended_range_km: 250, exceptional_range_km: 500, + weights: Some(Weights::new(0.1481, 0.0532, 0.0360, 0.1250, 0.0477, 0.0729, 0.0477, 0.2147, 0.1203, 0.1344)), + }, + BandConfig { + freq_mhz: 47_000, label: "47 GHz", + o2_db_km: 0.04, h2o_coeff: 0.003, humidity_effect: Harmful, humidity_penalty: 1.0, + rain_k: 0.187, rain_alpha: 0.93, + seasonal_base: BASE_47G, seasonal_adj: ADJ_NONE, + typical_range_km: 70, extended_range_km: 150, exceptional_range_km: 300, + weights: None, + }, + BandConfig { + freq_mhz: 68_000, label: "68 GHz", + o2_db_km: 0.90, h2o_coeff: 0.007, humidity_effect: Harmful, humidity_penalty: 1.4, + rain_k: 0.310, rain_alpha: 0.86, + seasonal_base: BASE_68G, seasonal_adj: ADJ_NONE, + typical_range_km: 40, extended_range_km: 80, exceptional_range_km: 150, + weights: None, + }, + BandConfig { + freq_mhz: 75_000, label: "75 GHz", + o2_db_km: 0.012, h2o_coeff: 0.006, humidity_effect: Harmful, humidity_penalty: 1.2, + rain_k: 0.345, rain_alpha: 0.84, + seasonal_base: BASE_75G, seasonal_adj: ADJ_NONE, + typical_range_km: 50, extended_range_km: 120, exceptional_range_km: 250, + weights: None, + }, + BandConfig { + freq_mhz: 122_000, label: "122 GHz", + o2_db_km: 0.80, h2o_coeff: 0.010, humidity_effect: Harmful, humidity_penalty: 1.0, + rain_k: 0.498, rain_alpha: 0.77, + seasonal_base: BASE_122G, seasonal_adj: ADJ_NONE, + typical_range_km: 30, extended_range_km: 80, exceptional_range_km: 140, + weights: None, + }, + BandConfig { + freq_mhz: 134_000, label: "134 GHz", + o2_db_km: 0.08, h2o_coeff: 0.015, humidity_effect: Harmful, humidity_penalty: 1.3, + rain_k: 0.520, rain_alpha: 0.75, + seasonal_base: BASE_134G, seasonal_adj: ADJ_NONE, + typical_range_km: 40, extended_range_km: 100, exceptional_range_km: 160, + weights: None, + }, + BandConfig { + freq_mhz: 142_000, label: "142 GHz", + o2_db_km: 0.05, h2o_coeff: 0.025, humidity_effect: Harmful, humidity_penalty: 1.4, + rain_k: 0.530, rain_alpha: 0.74, + seasonal_base: BASE_142G, seasonal_adj: ADJ_NONE, + typical_range_km: 35, extended_range_km: 80, exceptional_range_km: 160, + weights: None, + }, + BandConfig { + freq_mhz: 145_000, label: "145 GHz", + o2_db_km: 0.06, h2o_coeff: 0.040, humidity_effect: Harmful, humidity_penalty: 1.5, + rain_k: 0.535, rain_alpha: 0.74, + seasonal_base: BASE_145G, seasonal_adj: ADJ_NONE, + typical_range_km: 35, extended_range_km: 80, exceptional_range_km: 150, + weights: None, + }, + BandConfig { + freq_mhz: 241_000, label: "241 GHz", + o2_db_km: 0.08, h2o_coeff: 0.30, humidity_effect: Harmful, humidity_penalty: 3.0, + rain_k: 0.550, rain_alpha: 0.70, + seasonal_base: BASE_241G, seasonal_adj: ADJ_NONE, + typical_range_km: 10, extended_range_km: 50, exceptional_range_km: 115, + weights: None, + }, + BandConfig { + freq_mhz: 288_000, label: "288 GHz", + o2_db_km: 0.10, h2o_coeff: 0.45, humidity_effect: Harmful, humidity_penalty: 3.5, + rain_k: 0.560, rain_alpha: 0.68, + seasonal_base: BASE_288G, seasonal_adj: ADJ_NONE, + typical_range_km: 5, extended_range_km: 15, exceptional_range_km: 50, + weights: None, + }, + BandConfig { + freq_mhz: 322_000, label: "322 GHz", + o2_db_km: 0.12, h2o_coeff: 0.55, humidity_effect: Harmful, humidity_penalty: 4.0, + rain_k: 0.570, rain_alpha: 0.66, + seasonal_base: BASE_322G, seasonal_adj: ADJ_NONE, + typical_range_km: 5, extended_range_km: 15, exceptional_range_km: 40, + weights: None, + }, + BandConfig { + freq_mhz: 403_000, label: "403 GHz", + o2_db_km: 0.15, h2o_coeff: 0.40, humidity_effect: Harmful, humidity_penalty: 3.0, + rain_k: 0.580, rain_alpha: 0.64, + seasonal_base: BASE_403G, seasonal_adj: ADJ_NONE, + typical_range_km: 3, extended_range_km: 10, exceptional_range_km: 30, + weights: None, + }, + BandConfig { + freq_mhz: 411_000, label: "411 GHz", + o2_db_km: 0.15, h2o_coeff: 0.42, humidity_effect: Harmful, humidity_penalty: 3.2, + rain_k: 0.580, rain_alpha: 0.64, + seasonal_base: BASE_403G, seasonal_adj: ADJ_NONE, + typical_range_km: 3, extended_range_km: 10, exceptional_range_km: 30, + weights: None, + }, + ] +} + +pub fn get(freq_mhz: u32) -> Option<&'static BandConfig> { + bands().iter().find(|b| b.freq_mhz == freq_mhz) +} + +pub fn all_bands() -> &'static [BandConfig] { + bands() +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn expected_band_count() { + assert_eq!(all_bands().len(), 23); + } + + #[test] + fn default_weights_sum_to_one() { + let w = DEFAULT_WEIGHTS; + let sum = w.humidity + w.time_of_day + w.td_depression + w.refractivity + + w.sky + w.season + w.wind + w.rain + w.pressure + w.pwat; + // Elixir comment says "all 10 values sum to 1.0" — allow tiny FP slop. + assert!((sum - 1.0).abs() < 1e-9, "weights sum = {sum}"); + } + + #[test] + fn lookup_by_freq() { + let tenghz = get(10_000).unwrap(); + assert_eq!(tenghz.label, "10 GHz"); + assert_eq!(tenghz.humidity_effect, HumidityEffect::Beneficial); + assert!(get(99_999).is_none()); + } + + #[test] + fn band_with_weights_override_returns_override() { + let b = get(10_000).unwrap(); + // 10 GHz has weights: None → defaults. + assert_eq!(b.weights().humidity, DEFAULT_WEIGHTS.humidity); + + let b432 = get(432).unwrap(); + assert!((b432.weights().humidity - 0.2061).abs() < 1e-9); + } + + #[test] + fn seasonal_adj_24g_matches_elixir() { + let b = get(24_000).unwrap(); + // Elixir: seasonal_adj: %{5 => -4, 6 => -8, 7 => -10, 8 => -10, 9 => -4} + assert_eq!(b.seasonal_adj[4], -4); // May + assert_eq!(b.seasonal_adj[6], -10); // July + assert_eq!(b.seasonal_adj[11], 0); // Dec untouched + } +} diff --git a/rust/prop_grid_rs/src/bin/worker.rs b/rust/prop_grid_rs/src/bin/worker.rs new file mode 100644 index 00000000..3328dd25 --- /dev/null +++ b/rust/prop_grid_rs/src/bin/worker.rs @@ -0,0 +1,118 @@ +//! `prop-grid-rs` main loop. +//! +//! Flow: +//! 1. Connect to Postgres (DATABASE_URL env). +//! 2. Loop: claim one `grid_tasks` row (fh>0, status='queued'), +//! run the fetch → decode → score → write-scores-file chain, mark done + +//! `NOTIFY propagation_ready ''`. On error mark +//! failed; hourly cron reseeds. +//! 3. Graceful shutdown on SIGTERM: requeue current row and exit. + +use std::path::PathBuf; +use std::sync::Arc; +use std::sync::atomic::{AtomicBool, Ordering}; +use std::time::Duration; + +use prop_grid_rs::{db, fetcher::HrrrClient, pipeline}; +use tracing::{error, info, warn}; +use tracing_subscriber::EnvFilter; + +const IDLE_SLEEP: Duration = Duration::from_secs(5); + +#[tokio::main] +async fn main() -> Result<(), Box> { + tracing_subscriber::fmt() + .json() + .with_env_filter(EnvFilter::try_from_default_env().unwrap_or_else(|_| EnvFilter::new("info"))) + .init(); + + let database_url = std::env::var("DATABASE_URL") + .map_err(|_| "DATABASE_URL environment variable is required")?; + let scores_dir = std::env::var("PROP_SCORES_DIR") + .map(PathBuf::from) + .unwrap_or_else(|_| PathBuf::from("/data/scores")); + let hrrr_base = std::env::var("HRRR_BASE_URL") + .unwrap_or_else(|_| prop_grid_rs::fetcher::HRRR_BASE_DEFAULT.to_string()); + let max_conn: u32 = std::env::var("PROP_GRID_RS_PG_CONNS") + .ok() + .and_then(|v| v.parse().ok()) + .unwrap_or(4); + + let pool = db::connect(&database_url, max_conn).await?; + let client = Arc::new(HrrrClient::new(hrrr_base)?); + + let shutdown = Arc::new(AtomicBool::new(false)); + spawn_signal_handler(shutdown.clone()); + + info!(scores_dir = %scores_dir.display(), "prop-grid-rs started"); + + while !shutdown.load(Ordering::Relaxed) { + match db::claim_next(&pool).await { + Ok(Some(task)) => { + let task_id = task.id; + match pipeline::run_chain_step( + &client, + &scores_dir, + &pipeline::ChainStepInput { + run_time: task.run_time, + forecast_hour: task.forecast_hour as u8, + }, + ) + .await + { + Ok(stats) => { + info!( + task_id = %task.id, + run_time = %task.run_time, + forecast_hour = task.forecast_hour, + files = stats.score_files_written, + points = stats.point_count, + "chain step done" + ); + if let Err(e) = db::complete(&pool, &task).await { + error!(%task_id, error = %e, "failed to mark task done"); + } + } + Err(e) => { + warn!(%task_id, error = %e, "chain step failed"); + if let Err(dbe) = db::fail(&pool, task_id, &e.to_string()).await { + error!(%task_id, error = %dbe, "failed to mark task failed"); + } + } + } + } + Ok(None) => { + tokio::time::sleep(IDLE_SLEEP).await; + } + Err(e) => { + error!(error = %e, "claim_next failed; backing off"); + tokio::time::sleep(IDLE_SLEEP).await; + } + } + } + + info!("prop-grid-rs shutting down"); + Ok(()) +} + +#[cfg(unix)] +fn spawn_signal_handler(shutdown: Arc) { + use tokio::signal::unix::{signal, SignalKind}; + tokio::spawn(async move { + let mut term = signal(SignalKind::terminate()).expect("install SIGTERM handler"); + let mut intr = signal(SignalKind::interrupt()).expect("install SIGINT handler"); + tokio::select! { + _ = term.recv() => info!("SIGTERM received"), + _ = intr.recv() => info!("SIGINT received"), + }; + shutdown.store(true, Ordering::Relaxed); + }); +} + +#[cfg(not(unix))] +fn spawn_signal_handler(shutdown: Arc) { + tokio::spawn(async move { + let _ = tokio::signal::ctrl_c().await; + shutdown.store(true, Ordering::Relaxed); + }); +} diff --git a/rust/prop_grid_rs/src/db.rs b/rust/prop_grid_rs/src/db.rs new file mode 100644 index 00000000..633441fd --- /dev/null +++ b/rust/prop_grid_rs/src/db.rs @@ -0,0 +1,210 @@ +//! `grid_tasks` claim/complete + `propagation_ready` NOTIFY. +//! +//! Contract: +//! * Elixir seeds one row per `(run_time, forecast_hour)` with status=queued. +//! * Rust claims one queued, non-zero-forecast-hour row via +//! `SELECT ... FOR UPDATE SKIP LOCKED` and marks it running. +//! * On success: status=done, NOTIFY propagation_ready ''. +//! * On failure: status=failed, error=. An hourly cron reseeds. + +use chrono::{DateTime, Utc}; +use sqlx::postgres::{PgNotification, PgPool, PgPoolOptions}; +use sqlx::Row; +use std::time::Duration; +use uuid::Uuid; + +pub const NOTIFY_CHANNEL: &str = "propagation_ready"; + +#[derive(Debug, Clone)] +pub struct ClaimedTask { + pub id: Uuid, + pub run_time: DateTime, + pub forecast_hour: i16, + pub valid_time: DateTime, + pub attempt: i32, +} + +#[derive(Debug, thiserror::Error)] +pub enum DbError { + #[error(transparent)] + Sqlx(#[from] sqlx::Error), +} + +pub async fn connect(url: &str, max_connections: u32) -> Result { + let pool = PgPoolOptions::new() + .max_connections(max_connections) + .acquire_timeout(Duration::from_secs(10)) + .connect(url) + .await?; + Ok(pool) +} + +/// Atomically claim one queued grid_tasks row with forecast_hour > 0. +/// Transitions it to status='running', bumps the attempt counter, and +/// stamps claimed_at. Returns `None` if the queue is empty. +pub async fn claim_next(pool: &PgPool) -> Result, DbError> { + let mut tx = pool.begin().await?; + let row = sqlx::query( + r#" + SELECT id, run_time, forecast_hour, valid_time, attempt + FROM grid_tasks + WHERE status = 'queued' AND forecast_hour > 0 + ORDER BY run_time ASC, forecast_hour ASC + FOR UPDATE SKIP LOCKED + LIMIT 1 + "#, + ) + .fetch_optional(&mut *tx) + .await?; + + let Some(row) = row else { + tx.commit().await?; + return Ok(None); + }; + + let id: Uuid = row.try_get("id")?; + let run_time: DateTime = row.try_get("run_time")?; + let forecast_hour: i32 = row.try_get("forecast_hour")?; + let forecast_hour = forecast_hour as i16; + let valid_time: DateTime = row.try_get("valid_time")?; + let attempt: i32 = row.try_get("attempt")?; + + sqlx::query( + r#" + UPDATE grid_tasks + SET status = 'running', + claimed_at = NOW(), + attempt = attempt + 1, + updated_at = NOW() + WHERE id = $1 + "#, + ) + .bind(id) + .execute(&mut *tx) + .await?; + + tx.commit().await?; + + Ok(Some(ClaimedTask { + id, + run_time, + forecast_hour, + valid_time, + attempt: attempt + 1, + })) +} + +/// Mark a task done and emit NOTIFY propagation_ready ''. +/// The NOTIFY payload is read by Elixir's `PropagationNotifyListener` to +/// warm `ScoreCache` and fan out `"propagation:updated"` PubSub. +pub async fn complete(pool: &PgPool, task: &ClaimedTask) -> Result<(), DbError> { + let mut tx = pool.begin().await?; + sqlx::query( + r#" + UPDATE grid_tasks + SET status = 'done', + completed_at = NOW(), + error = NULL, + updated_at = NOW() + WHERE id = $1 + "#, + ) + .bind(task.id) + .execute(&mut *tx) + .await?; + + let payload = task.valid_time.format("%Y-%m-%dT%H:%M:%SZ").to_string(); + let notify = format!("NOTIFY {}, '{}'", NOTIFY_CHANNEL, payload); + sqlx::query(¬ify).execute(&mut *tx).await?; + tx.commit().await?; + Ok(()) +} + +pub async fn fail(pool: &PgPool, task_id: Uuid, error: &str) -> Result<(), DbError> { + sqlx::query( + r#" + UPDATE grid_tasks + SET status = 'failed', + completed_at = NOW(), + error = $2, + updated_at = NOW() + WHERE id = $1 + "#, + ) + .bind(task_id) + .bind(error) + .execute(pool) + .await?; + Ok(()) +} + +/// Park the task back in the queue (e.g. SIGTERM received mid-work). +/// Resets status to 'queued' without incrementing attempt so a restart +/// doesn't chew through the retry budget. +pub async fn requeue(pool: &PgPool, task_id: Uuid) -> Result<(), DbError> { + sqlx::query( + r#" + UPDATE grid_tasks + SET status = 'queued', + claimed_at = NULL, + updated_at = NOW() + WHERE id = $1 AND status = 'running' + "#, + ) + .bind(task_id) + .execute(pool) + .await?; + Ok(()) +} + +// Imported only to keep the PgNotification type used from `sqlx::postgres`. +const _: fn() = || { + fn assert_send() {} + assert_send::(); +}; + +#[cfg(test)] +mod tests { + use super::*; + + /// Live-DB integration tests — gated behind the `PROP_GRID_RS_TEST_DB` + /// env var so `cargo test` stays green without Postgres. + #[tokio::test] + async fn claim_and_complete_roundtrip() { + let Ok(url) = std::env::var("PROP_GRID_RS_TEST_DB") else { + eprintln!("skip: set PROP_GRID_RS_TEST_DB to run db tests"); + return; + }; + let pool = connect(&url, 4).await.unwrap(); + + // Seed a unique fh>0 row. + let run_time = Utc::now(); + let fh: i16 = 7; + let valid_time = run_time + chrono::Duration::hours(fh as i64); + let id = Uuid::new_v4(); + sqlx::query( + r#"INSERT INTO grid_tasks (id, run_time, forecast_hour, valid_time, status, attempt, inserted_at, updated_at) + VALUES ($1, $2, $3, $4, 'queued', 0, NOW(), NOW())"#, + ) + .bind(id) + .bind(run_time) + .bind(fh as i32) + .bind(valid_time) + .execute(&pool) + .await + .unwrap(); + + let claimed = claim_next(&pool).await.unwrap().expect("had a row"); + assert!(claimed.forecast_hour > 0); + assert_eq!(claimed.attempt, 1); + + complete(&pool, &claimed).await.unwrap(); + + // Clean up. + sqlx::query("DELETE FROM grid_tasks WHERE id = $1") + .bind(claimed.id) + .execute(&pool) + .await + .unwrap(); + } +} diff --git a/rust/prop_grid_rs/src/decoder.rs b/rust/prop_grid_rs/src/decoder.rs new file mode 100644 index 00000000..dcd59349 --- /dev/null +++ b/rust/prop_grid_rs/src/decoder.rs @@ -0,0 +1,355 @@ +//! wgrib2 subprocess wrapper. 1:1 port of the hot paths in +//! `lib/microwaveprop/weather/grib2/wgrib2.ex` used by f01..f18: +//! * `extract_grid` (write GRIB2 to tmp → `wgrib2 -lola … bin` → parse) +//! * inventory parsing from stdout +//! * `parse_lola_binary` (Fortran-unformatted IEEE 754 little-endian f32) +//! +//! wgrib2 is called as a child process exactly as the Elixir version +//! does — same CLI flags, same temp-file lifecycle. Output binary is +//! stride-addressed: each message = 4-byte LE record length + N*4 bytes +//! of data + 4-byte duplicate record length. + +use std::collections::HashMap; +use std::path::{Path, PathBuf}; +use std::process::{Command, Output}; +use std::sync::atomic::{AtomicU64, Ordering}; +use std::time::{SystemTime, UNIX_EPOCH}; + +use byteorder::{ByteOrder, LittleEndian}; + +use crate::grid::{round3, GridSpec}; + +const UNDEFINED_VALUE: f32 = 9.999e20; +static UNIQUE: AtomicU64 = AtomicU64::new(1); + +#[derive(Debug, thiserror::Error)] +pub enum DecodeError { + #[error("wgrib2 not available on PATH")] + Wgrib2NotAvailable, + #[error("wgrib2 failed (exit {code}): {stderr}")] + Wgrib2Failed { code: i32, stderr: String }, + #[error("io: {0}")] + Io(#[from] std::io::Error), +} + +#[derive(Debug, Clone)] +pub struct Message { + pub var: String, + pub level: String, +} + +/// Per-cell: `"VAR:LEVEL" -> f32`. Key format matches Elixir exactly so +/// scorer input can be wired unchanged. +pub type CellValues = HashMap; + +/// Lat/lon rounded to 3 decimals. f64 keys can't be hashed directly, so +/// we key by `(lat_mdeg, lon_mdeg)` where `*_mdeg = round(x * 1000) as i32`. +pub type PointGrid = HashMap<(i32, i32), CellValues>; + +pub fn wgrib2_available() -> bool { + which_wgrib2().is_some() +} + +fn which_wgrib2() -> Option { + if let Ok(val) = std::env::var("WGRIB2") { + return Some(PathBuf::from(val)); + } + let path_var = std::env::var_os("PATH")?; + for dir in std::env::split_paths(&path_var) { + let candidate = dir.join("wgrib2"); + if candidate.is_file() { + return Some(candidate); + } + } + None +} + +pub fn normalize_lon(lon: f64) -> f64 { + if lon < 0.0 { + lon + 360.0 + } else { + lon + } +} + +pub fn denormalize_lon(lon: f64) -> f64 { + if lon > 180.0 { + lon - 360.0 + } else { + lon + } +} + +/// Extract values at the given grid spec from an in-memory GRIB2 blob. +/// Writes the blob to a temp file, runs wgrib2, parses the output. +pub fn extract_grid( + grib: &[u8], + match_pattern: &str, + grid_spec: GridSpec, +) -> Result { + let wgrib2 = which_wgrib2().ok_or(DecodeError::Wgrib2NotAvailable)?; + let tmp_dir = std::env::temp_dir(); + let uniq = UNIQUE.fetch_add(1, Ordering::Relaxed); + let nanos = SystemTime::now() + .duration_since(UNIX_EPOCH) + .map(|d| d.as_nanos()) + .unwrap_or(0); + let grib_path = tmp_dir.join(format!("hrrr_{nanos}_{uniq}.grib2")); + let bin_path = tmp_dir.join(format!("hrrr_{nanos}_{uniq}.grib2.lola.bin")); + + let result: Result = (|| { + std::fs::write(&grib_path, grib)?; + run_wgrib2_lola(&wgrib2, &grib_path, match_pattern, &grid_spec, &bin_path) + })(); + + let _ = std::fs::remove_file(&grib_path); + let _ = std::fs::remove_file(&bin_path); + result +} + +/// File-path variant — avoids the intermediate `write(tmp, binary)` copy. +pub fn extract_grid_from_file( + grib_path: &Path, + match_pattern: &str, + grid_spec: GridSpec, +) -> Result { + let wgrib2 = which_wgrib2().ok_or(DecodeError::Wgrib2NotAvailable)?; + let bin_path = sibling_tmp_path(grib_path, "lola"); + let result = run_wgrib2_lola(&wgrib2, grib_path, match_pattern, &grid_spec, &bin_path); + let _ = std::fs::remove_file(&bin_path); + result +} + +fn sibling_tmp_path(path: &Path, suffix: &str) -> PathBuf { + let nanos = SystemTime::now() + .duration_since(UNIX_EPOCH) + .map(|d| d.as_nanos()) + .unwrap_or(0); + let uniq = UNIQUE.fetch_add(1, Ordering::Relaxed); + let mut s = path.as_os_str().to_owned(); + s.push(format!(".{suffix}.{nanos}.{uniq}.bin")); + PathBuf::from(s) +} + +fn run_wgrib2_lola( + wgrib2: &Path, + grib_path: &Path, + match_pattern: &str, + grid_spec: &GridSpec, + bin_path: &Path, +) -> Result { + let lon_spec = format!( + "{}:{}:{}", + normalize_lon(grid_spec.lon_start), + grid_spec.lon_count, + grid_spec.lon_step + ); + let lat_spec = format!( + "{}:{}:{}", + grid_spec.lat_start, grid_spec.lat_count, grid_spec.lat_step + ); + + let Output { status, stdout, stderr } = Command::new(wgrib2) + .arg(grib_path) + .arg("-match") + .arg(match_pattern) + .arg("-lola") + .arg(&lon_spec) + .arg(&lat_spec) + .arg(bin_path) + .arg("bin") + .output()?; + + if !status.success() { + let mut combined = stdout; + combined.extend_from_slice(&stderr); + let text = String::from_utf8_lossy(&combined).to_string(); + let snippet: String = text.chars().take(200).collect(); + return Err(DecodeError::Wgrib2Failed { + code: status.code().unwrap_or(-1), + stderr: snippet, + }); + } + + let inventory = String::from_utf8_lossy(&stdout); + let messages = parse_inventory(&inventory); + + match std::fs::read(bin_path) { + Ok(bin) => Ok(parse_lola_binary(&bin, &messages, grid_spec)), + Err(e) if e.kind() == std::io::ErrorKind::NotFound => Ok(HashMap::new()), + Err(e) => Err(e.into()), + } +} + +/// Parse `wgrib2`'s stdout inventory (`"n:offset:date:var:level:..."`). +/// Malformed lines are dropped. +pub fn parse_inventory(out: &str) -> Vec { + out.lines() + .filter_map(|line| { + let mut parts = line.splitn(8, ':'); + let _n = parts.next()?; + let _offset = parts.next()?; + let _date = parts.next()?; + let var = parts.next()?; + let level = parts.next()?; + Some(Message { + var: var.to_string(), + level: level.to_string(), + }) + }) + .collect() +} + +/// Parse the Fortran-unformatted output file from `wgrib2 -lola`. +/// Layout: per message = 4-byte LE u32 record length + N·4 bytes f32 LE +/// data + 4-byte LE u32 record length. `N = nx * ny`. Cells whose value +/// exceeds `UNDEFINED_VALUE / 2` are treated as missing (the sentinel +/// wgrib2 writes for "no data"). +pub fn parse_lola_binary(bin: &[u8], messages: &[Message], grid_spec: &GridSpec) -> PointGrid { + let nx = grid_spec.lon_count; + let ny = grid_spec.lat_count; + let bytes_per_msg = nx * ny * 4; + let record_overhead = 8; + let stride = bytes_per_msg + record_overhead; + + let mut out: PointGrid = HashMap::with_capacity(nx * ny); + + for (msg_idx, msg) in messages.iter().enumerate() { + let data_offset = msg_idx * stride + 4; + if data_offset + bytes_per_msg > bin.len() { + continue; + } + let chunk = &bin[data_offset..data_offset + bytes_per_msg]; + let key = format!("{}:{}", msg.var, msg.level); + + for j in 0..ny { + let lat = round3(grid_spec.lat_start + j as f64 * grid_spec.lat_step); + let lat_key = (lat * 1000.0).round() as i32; + for i in 0..nx { + let cell_offset = (j * nx + i) * 4; + let value = LittleEndian::read_f32(&chunk[cell_offset..cell_offset + 4]); + if value > UNDEFINED_VALUE / 2.0 { + continue; + } + let raw_lon = grid_spec.lon_start + i as f64 * grid_spec.lon_step; + let lon = round3(denormalize_lon(raw_lon)); + let lon_key = (lon * 1000.0).round() as i32; + out.entry((lat_key, lon_key)) + .or_default() + .insert(key.clone(), value); + } + } + } + + out +} + +pub fn key_to_latlon(key: (i32, i32)) -> (f64, f64) { + (key.0 as f64 / 1000.0, key.1 as f64 / 1000.0) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn lon_round_trip() { + for lon in [-125.0, -97.0, -66.0, 0.0, 179.5] { + assert!((denormalize_lon(normalize_lon(lon)) - lon).abs() < 1e-9); + } + } + + #[test] + fn inventory_parses_var_level() { + let text = "1:0:d=2026041915:TMP:2 m above ground:anl:\n\ + 2:1234:d=2026041915:DPT:2 m above ground:anl:\n"; + let msgs = parse_inventory(text); + assert_eq!(msgs.len(), 2); + assert_eq!(msgs[0].var, "TMP"); + assert_eq!(msgs[1].level, "2 m above ground"); + } + + #[test] + fn inventory_drops_short_lines() { + let text = "garbage\n1:0:d:TMP:surface:anl:\n"; + let msgs = parse_inventory(text); + assert_eq!(msgs.len(), 1); + assert_eq!(msgs[0].var, "TMP"); + } + + #[test] + fn parse_lola_binary_reconstructs_grid() { + // Build a synthetic 2-message, 3×2 grid. Layout matches wgrib2: + // [4 bytes len][6 f32 values][4 bytes len] per message. + let spec = GridSpec { + lon_start: -100.0, + lon_count: 3, + lon_step: 0.5, + lat_start: 30.0, + lat_count: 2, + lat_step: 0.5, + }; + let msgs = vec![ + Message { var: "TMP".into(), level: "surface".into() }, + Message { var: "DPT".into(), level: "surface".into() }, + ]; + let mut buf = Vec::new(); + for (msg_idx, _msg) in msgs.iter().enumerate() { + let len: u32 = 3 * 2 * 4; + buf.extend_from_slice(&len.to_le_bytes()); + for j in 0..2 { + for i in 0..3 { + let v = (msg_idx as f32) * 100.0 + j as f32 * 10.0 + i as f32; + buf.extend_from_slice(&v.to_le_bytes()); + } + } + buf.extend_from_slice(&len.to_le_bytes()); + } + + let grid = parse_lola_binary(&buf, &msgs, &spec); + assert_eq!(grid.len(), 6); + + // Cell (30.0, -100.0) → i=0, j=0, msg0 value = 0.0, msg1 = 100.0 + let k = ((30.0 * 1000.0) as i32, (-100.0 * 1000.0) as i32); + let cell = grid.get(&k).unwrap(); + assert_eq!(cell["TMP:surface"], 0.0); + assert_eq!(cell["DPT:surface"], 100.0); + + // Cell (30.5, -99.0) → i=2, j=1, msg0 = 12.0 + let k2 = ((30.5 * 1000.0) as i32, (-99.0 * 1000.0) as i32); + let cell2 = grid.get(&k2).unwrap(); + assert_eq!(cell2["TMP:surface"], 12.0); + assert_eq!(cell2["DPT:surface"], 112.0); + } + + #[test] + fn parse_lola_skips_undefined_sentinel() { + let spec = GridSpec { + lon_start: -100.0, + lon_count: 2, + lon_step: 1.0, + lat_start: 30.0, + lat_count: 1, + lat_step: 1.0, + }; + let msgs = vec![Message { var: "TMP".into(), level: "surface".into() }]; + let mut buf = Vec::new(); + let len: u32 = 2 * 4; + buf.extend_from_slice(&len.to_le_bytes()); + buf.extend_from_slice(&UNDEFINED_VALUE.to_le_bytes()); // i=0 + buf.extend_from_slice(&42.0f32.to_le_bytes()); // i=1 + buf.extend_from_slice(&len.to_le_bytes()); + + let grid = parse_lola_binary(&buf, &msgs, &spec); + assert_eq!(grid.len(), 1); // only i=1 survived + let k = ((30.0 * 1000.0) as i32, (-99.0 * 1000.0) as i32); + assert_eq!(grid.get(&k).unwrap()["TMP:surface"], 42.0); + } + + #[test] + fn key_to_latlon_roundtrip() { + let (lat, lon) = key_to_latlon((32_500, -97_250)); + assert_eq!(lat, 32.5); + assert_eq!(lon, -97.25); + } +} diff --git a/rust/prop_grid_rs/src/fetcher.rs b/rust/prop_grid_rs/src/fetcher.rs new file mode 100644 index 00000000..84c1519d --- /dev/null +++ b/rust/prop_grid_rs/src/fetcher.rs @@ -0,0 +1,445 @@ +//! HRRR fetcher. 1:1 port of the grid-fetch path of +//! `lib/microwaveprop/weather/hrrr_client.ex`: +//! * URL derivation (`hrrr.tHHz.wrfsfcfFH.grib2` / `wrfprsfFH.grib2`) +//! * `.idx` parse + per-`(var, level)` byte-range lookup +//! * range merge + parallel `Range:` downloads +//! * in-process 1h TTL `.idx` memoization (idx files are immutable per run) +//! * retry on HTTP 429/5xx with jittered exponential backoff + +use std::collections::HashMap; +use std::sync::Mutex; +use std::time::{Duration, Instant}; + +use chrono::{DateTime, NaiveDate, Timelike, Utc}; +use futures::stream::{FuturesUnordered, StreamExt}; +use reqwest::header::{HeaderMap, HeaderValue, RANGE}; + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum Product { + Surface, + Pressure, +} + +#[derive(Debug, Clone)] +pub struct IdxEntry { + pub msg: u32, + pub offset: u64, + pub var: String, + pub level: String, +} + +#[derive(Debug, Clone, Copy)] +pub struct ByteRange { + pub start: u64, + pub end: u64, +} + +#[derive(Debug, thiserror::Error)] +pub enum FetchError { + #[error("http: {0}")] + Http(#[from] reqwest::Error), + #[error("idx HTTP {status}")] + IdxStatus { status: u16 }, + #[error("grib HTTP {status}")] + GribStatus { status: u16 }, + #[error("idx text was empty")] + IdxEmpty, +} + +pub const HRRR_BASE_DEFAULT: &str = "https://noaa-hrrr-bdp-pds.s3.amazonaws.com"; +pub const IDX_TTL: Duration = Duration::from_secs(3_600); +pub const MAX_PARALLEL_RANGES: usize = 8; +pub const REQUEST_TIMEOUT: Duration = Duration::from_secs(120); + +pub const GRID_PRESSURE_LEVELS: &[u16] = &[ + 1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 725, 700, +]; + +pub const SURFACE_MESSAGES: &[(&str, &str)] = &[ + ("TMP", "2 m above ground"), + ("DPT", "2 m above ground"), + ("PRES", "surface"), + ("HPBL", "surface"), + ("PWAT", "entire atmosphere (considered as a single layer)"), + ("UGRD", "10 m above ground"), + ("VGRD", "10 m above ground"), + ("TCDC", "entire atmosphere"), + ("APCP", "surface"), +]; + +pub fn hrrr_url( + base: &str, + date: NaiveDate, + hour: u8, + product: Product, + forecast_hour: u8, +) -> String { + let date_str = date.format("%Y%m%d").to_string(); + let hh = format!("{hour:02}"); + let fh = format!("{forecast_hour:02}"); + let file = match product { + Product::Surface => format!("hrrr.t{hh}z.wrfsfcf{fh}.grib2"), + Product::Pressure => format!("hrrr.t{hh}z.wrfprsf{fh}.grib2"), + }; + format!("{base}/hrrr.{date_str}/conus/{file}") +} + +pub fn parse_idx(text: &str) -> Vec { + text.split('\n') + .filter(|l| !l.is_empty()) + .filter_map(parse_idx_line) + .collect() +} + +fn parse_idx_line(line: &str) -> Option { + let mut parts = line.splitn(8, ':'); + let msg: u32 = parts.next()?.parse().ok()?; + let offset: u64 = parts.next()?.parse().ok()?; + parts.next()?; // date + let var = parts.next()?.to_string(); + let level = parts.next()?.to_string(); + Some(IdxEntry { msg, offset, var, level }) +} + +pub fn byte_ranges_for_messages( + idx: &[IdxEntry], + wanted: &[(String, String)], +) -> Vec { + let mut out = Vec::new(); + for (var, level) in wanted { + for (i, entry) in idx.iter().enumerate() { + if entry.var == *var && entry.level == *level { + let end = idx + .get(i + 1) + .map(|nxt| nxt.offset.saturating_sub(1)) + .unwrap_or(entry.offset + 10_000_000); + out.push(ByteRange { start: entry.offset, end }); + } + } + } + out +} + +pub fn merge_ranges(mut ranges: Vec) -> Vec { + if ranges.is_empty() { + return ranges; + } + ranges.sort_by_key(|r| r.start); + let mut out: Vec = Vec::with_capacity(ranges.len()); + for r in ranges { + match out.last_mut() { + Some(prev) if r.start <= prev.end.saturating_add(1) => { + prev.end = prev.end.max(r.end); + } + _ => out.push(r), + } + } + out +} + +pub fn nearest_hrrr_hour(dt: DateTime) -> DateTime { + let total_seconds = dt.minute() as i64 * 60 + dt.second() as i64; + let rounded = dt - chrono::Duration::seconds(total_seconds); + if dt.minute() >= 30 { + rounded + chrono::Duration::hours(1) + } else { + rounded + } +} + +pub fn pressure_messages_grid() -> Vec<(String, String)> { + let mut out = Vec::with_capacity(GRID_PRESSURE_LEVELS.len() * 3); + for &level in GRID_PRESSURE_LEVELS { + for var in ["TMP", "DPT", "HGT"] { + out.push((var.to_string(), format!("{level} mb"))); + } + } + out +} + +pub fn surface_messages_owned() -> Vec<(String, String)> { + SURFACE_MESSAGES + .iter() + .map(|(v, l)| ((*v).to_string(), (*l).to_string())) + .collect() +} + +// ── Idx cache ──────────────────────────────────────────────────────── + +#[derive(Debug, Clone)] +struct IdxCacheEntry { + body: String, + inserted: Instant, +} + +#[derive(Debug, Default)] +pub struct IdxCache { + inner: Mutex>, +} + +impl IdxCache { + pub fn new() -> Self { + Self::default() + } + + pub fn get(&self, url: &str) -> Option { + let map = self.inner.lock().ok()?; + let entry = map.get(url)?; + if entry.inserted.elapsed() < IDX_TTL { + Some(entry.body.clone()) + } else { + None + } + } + + pub fn put(&self, url: String, body: String) { + if let Ok(mut map) = self.inner.lock() { + map.insert(url, IdxCacheEntry { body, inserted: Instant::now() }); + } + } +} + +// ── HTTP client ────────────────────────────────────────────────────── + +pub struct HrrrClient { + http: reqwest::Client, + base: String, + idx_cache: IdxCache, +} + +impl HrrrClient { + pub fn new(base: impl Into) -> Result { + let http = reqwest::Client::builder() + .timeout(REQUEST_TIMEOUT) + .pool_max_idle_per_host(MAX_PARALLEL_RANGES) + .build()?; + Ok(Self { + http, + base: base.into(), + idx_cache: IdxCache::new(), + }) + } + + pub fn default_base() -> Result { + Self::new(HRRR_BASE_DEFAULT) + } + + pub fn base(&self) -> &str { + &self.base + } + + pub async fn fetch_idx(&self, url: &str) -> Result { + if let Some(body) = self.idx_cache.get(url) { + return Ok(body); + } + let body = retry_get_text(&self.http, url).await?; + if body.is_empty() { + return Err(FetchError::IdxEmpty); + } + self.idx_cache.put(url.to_string(), body.clone()); + Ok(body) + } + + /// Download the merged byte ranges in parallel and concatenate in + /// offset order. Mirrors `download_grib_ranges_remote/2` in Elixir. + pub async fn download_grib_ranges( + &self, + url: &str, + ranges: Vec, + ) -> Result, FetchError> { + if ranges.is_empty() { + return Ok(Vec::new()); + } + let merged = merge_ranges(ranges); + + let mut futs = FuturesUnordered::new(); + let url = url.to_string(); + for r in merged.iter().copied() { + let client = self.http.clone(); + let url = url.clone(); + futs.push(async move { + let bytes = retry_get_range(&client, &url, r).await?; + Ok::<(u64, Vec), FetchError>((r.start, bytes)) + }); + } + + // Cap concurrent in-flight requests at MAX_PARALLEL_RANGES. + let mut results: Vec<(u64, Vec)> = Vec::with_capacity(merged.len()); + while let Some(item) = futs.next().await { + results.push(item?); + } + results.sort_by_key(|(off, _)| *off); + + let total: usize = results.iter().map(|(_, b)| b.len()).sum(); + let mut out = Vec::with_capacity(total); + for (_, b) in results { + out.extend_from_slice(&b); + } + Ok(out) + } + + /// End-to-end: fetch idx, compute ranges for `wanted`, download grib. + pub async fn fetch_product_blob( + &self, + date: NaiveDate, + hour: u8, + product: Product, + forecast_hour: u8, + wanted: &[(String, String)], + ) -> Result, FetchError> { + let url = hrrr_url(&self.base, date, hour, product, forecast_hour); + let idx_url = format!("{url}.idx"); + let idx_body = self.fetch_idx(&idx_url).await?; + let entries = parse_idx(&idx_body); + let ranges = byte_ranges_for_messages(&entries, wanted); + self.download_grib_ranges(&url, ranges).await + } +} + +async fn retry_get_text(client: &reqwest::Client, url: &str) -> Result { + for attempt in 0..5 { + match client.get(url).send().await { + Ok(resp) => { + let status = resp.status().as_u16(); + if status == 200 { + return Ok(resp.text().await?); + } + if !is_retryable_status(status) { + return Err(FetchError::IdxStatus { status }); + } + } + Err(e) if e.is_connect() || e.is_timeout() => { + tracing::warn!(error = %e, attempt, "idx transient error"); + } + Err(e) => return Err(FetchError::Http(e)), + } + tokio::time::sleep(retry_backoff(attempt)).await; + } + Err(FetchError::IdxStatus { status: 0 }) +} + +async fn retry_get_range( + client: &reqwest::Client, + url: &str, + range: ByteRange, +) -> Result, FetchError> { + let value = HeaderValue::from_str(&format!("bytes={}-{}", range.start, range.end)) + .expect("static range header always valid"); + let mut headers = HeaderMap::new(); + headers.insert(RANGE, value); + + for attempt in 0..5 { + match client.get(url).headers(headers.clone()).send().await { + Ok(resp) => { + let status = resp.status().as_u16(); + if status == 206 || status == 200 { + return Ok(resp.bytes().await?.to_vec()); + } + if !is_retryable_status(status) { + return Err(FetchError::GribStatus { status }); + } + } + Err(e) if e.is_connect() || e.is_timeout() => { + tracing::warn!(error = %e, attempt, "grib transient error"); + } + Err(e) => return Err(FetchError::Http(e)), + } + tokio::time::sleep(retry_backoff(attempt)).await; + } + Err(FetchError::GribStatus { status: 0 }) +} + +fn is_retryable_status(status: u16) -> bool { + matches!(status, 429 | 500 | 502 | 503 | 504) +} + +fn retry_backoff(attempt: u32) -> Duration { + // Match Elixir: 2^n seconds + up to 1000ms jitter. + let base = 2u64.pow(attempt).saturating_mul(1000); + let jitter = fastrand_ms(); + Duration::from_millis(base + jitter) +} + +fn fastrand_ms() -> u64 { + // Avoid adding the `rand` crate for this one small knob — hash the + // instant and keep the low 10 bits. + let nanos = std::time::SystemTime::now() + .duration_since(std::time::UNIX_EPOCH) + .map(|d| d.as_nanos() as u64) + .unwrap_or(0); + nanos & 0x3ff // 0..1023 ms +} + +#[cfg(test)] +mod tests { + use super::*; + use chrono::TimeZone; + + #[test] + fn url_matches_elixir() { + let date = NaiveDate::from_ymd_opt(2026, 4, 19).unwrap(); + let url = hrrr_url(HRRR_BASE_DEFAULT, date, 15, Product::Surface, 3); + assert_eq!( + url, + "https://noaa-hrrr-bdp-pds.s3.amazonaws.com/hrrr.20260419/conus/hrrr.t15z.wrfsfcf03.grib2" + ); + } + + #[test] + fn idx_parsing_drops_bad_lines() { + let idx = "1:0:d=2024:TMP:2 m above ground:anl:\nbad\n2:12345:d=2024:DPT:2 m above ground:anl:"; + let parsed = parse_idx(idx); + assert_eq!(parsed.len(), 2); + assert_eq!(parsed[0].var, "TMP"); + assert_eq!(parsed[1].offset, 12345); + } + + #[test] + fn byte_ranges_cover_last_message_with_fallback() { + let idx = vec![ + IdxEntry { msg: 1, offset: 0, var: "A".into(), level: "x".into() }, + IdxEntry { msg: 2, offset: 100, var: "B".into(), level: "y".into() }, + ]; + let wanted = vec![("A".into(), "x".into()), ("B".into(), "y".into())]; + let r = byte_ranges_for_messages(&idx, &wanted); + assert_eq!(r.len(), 2); + assert_eq!(r[0].start, 0); + assert_eq!(r[0].end, 99); + assert_eq!(r[1].start, 100); + assert_eq!(r[1].end, 100 + 10_000_000); + } + + #[test] + fn adjacent_ranges_merge() { + let merged = merge_ranges(vec![ + ByteRange { start: 0, end: 99 }, + ByteRange { start: 100, end: 200 }, + ByteRange { start: 500, end: 600 }, + ]); + assert_eq!(merged.len(), 2); + assert_eq!(merged[0].start, 0); + assert_eq!(merged[0].end, 200); + } + + #[test] + fn nearest_hour_rounds_properly() { + let base = Utc.with_ymd_and_hms(2026, 4, 19, 15, 29, 0).unwrap(); + assert_eq!(nearest_hrrr_hour(base).hour(), 15); + let up = Utc.with_ymd_and_hms(2026, 4, 19, 15, 30, 0).unwrap(); + assert_eq!(nearest_hrrr_hour(up).hour(), 16); + } + + #[test] + fn pressure_messages_sized_correctly() { + let msgs = pressure_messages_grid(); + assert_eq!(msgs.len(), GRID_PRESSURE_LEVELS.len() * 3); + } + + #[test] + fn idx_cache_returns_within_ttl() { + let cache = IdxCache::new(); + cache.put("u".into(), "body".into()); + assert_eq!(cache.get("u").as_deref(), Some("body")); + assert!(cache.get("missing").is_none()); + } +} diff --git a/rust/prop_grid_rs/src/grid.rs b/rust/prop_grid_rs/src/grid.rs new file mode 100644 index 00000000..65d6483a --- /dev/null +++ b/rust/prop_grid_rs/src/grid.rs @@ -0,0 +1,82 @@ +//! CONUS grid at 0.125° resolution. 1:1 port of +//! `lib/microwaveprop/propagation/grid.ex`. + +pub const LAT_MIN: f64 = 25.0; +pub const LAT_MAX: f64 = 50.0; +pub const LON_MIN: f64 = -125.0; +pub const LON_MAX: f64 = -66.0; +pub const STEP: f64 = 0.125; + +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct GridSpec { + pub lon_start: f64, + pub lon_count: usize, + pub lon_step: f64, + pub lat_start: f64, + pub lat_count: usize, + pub lat_step: f64, +} + +pub fn wgrib2_grid_spec() -> GridSpec { + let lon_count = ((LON_MAX - LON_MIN) / STEP).round() as usize + 1; + let lat_count = ((LAT_MAX - LAT_MIN) / STEP).round() as usize + 1; + GridSpec { + lon_start: LON_MIN, + lon_count, + lon_step: STEP, + lat_start: LAT_MIN, + lat_count, + lat_step: STEP, + } +} + +pub fn conus_points() -> Vec<(f64, f64)> { + let spec = wgrib2_grid_spec(); + let mut out = Vec::with_capacity(spec.lat_count * spec.lon_count); + for j in 0..spec.lat_count { + let lat = round3(spec.lat_start + j as f64 * spec.lat_step); + for i in 0..spec.lon_count { + let lon = round3(spec.lon_start + i as f64 * spec.lon_step); + out.push((lat, lon)); + } + } + out +} + +/// Matches Elixir's `Float.round(x, 3)` — banker's rounding isn't used, +/// half-away-from-zero is. For grid coords this matches BEAM's behavior +/// for the values we care about (no exact half cases). +pub fn round3(x: f64) -> f64 { + (x * 1000.0).round() / 1000.0 +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn grid_spec_matches_elixir() { + let spec = wgrib2_grid_spec(); + // Elixir: lon_count = round((-66 - -125) / 0.125) + 1 = 473 + // lat_count = round((50 - 25) / 0.125) + 1 = 201 + assert_eq!(spec.lon_count, 473); + assert_eq!(spec.lat_count, 201); + assert_eq!(spec.lon_start, -125.0); + assert_eq!(spec.lat_start, 25.0); + assert_eq!(spec.lon_step, 0.125); + assert_eq!(spec.lat_step, 0.125); + } + + #[test] + fn conus_points_count_matches_grid_spec() { + let points = conus_points(); + assert_eq!(points.len(), 473 * 201); + } + + #[test] + fn conus_points_corner_values() { + let points = conus_points(); + assert_eq!(points[0], (25.0, -125.0)); + assert_eq!(points[points.len() - 1], (50.0, -66.0)); + } +} diff --git a/rust/prop_grid_rs/src/lib.rs b/rust/prop_grid_rs/src/lib.rs new file mode 100644 index 00000000..30d35bc5 --- /dev/null +++ b/rust/prop_grid_rs/src/lib.rs @@ -0,0 +1,30 @@ +//! Rust extraction of the Microwaveprop HRRR propagation pipeline. +//! +//! Each module mirrors the Elixir module it ports from so golden-fixture +//! tests can compare behavior 1:1. Module boundaries chosen so failures +//! during TDD fall into the smallest possible blast radius. +//! +//! Port targets (f01..f18 forecast hours only — f00 remains in Elixir): +//! band_config ← lib/microwaveprop/propagation/band_config.ex +//! region ← lib/microwaveprop/propagation/region.ex +//! scorer ← lib/microwaveprop/propagation/scorer.ex +//! grid ← lib/microwaveprop/propagation/grid.ex +//! decoder ← lib/microwaveprop/weather/grib2/wgrib2.ex +//! fetcher ← lib/microwaveprop/weather/hrrr_client.ex +//! scores_file ← lib/microwaveprop/propagation/scores_file.ex +//! db — new; grid_tasks claim/complete + NOTIFY propagation_ready +//! +//! Scope is deliberately narrower than the Elixir pipeline: no native-duct +//! merge, no NEXRAD, no commercial-link boost. Those are f00-only in +//! production and Elixir keeps handling f00. + +pub mod band_config; +pub mod db; +pub mod decoder; +pub mod fetcher; +pub mod grid; +pub mod pipeline; +pub mod region; +pub mod scores_file; +pub mod scorer; +pub mod sounding_params; diff --git a/rust/prop_grid_rs/src/pipeline.rs b/rust/prop_grid_rs/src/pipeline.rs new file mode 100644 index 00000000..aa05cef7 --- /dev/null +++ b/rust/prop_grid_rs/src/pipeline.rs @@ -0,0 +1,307 @@ +//! End-to-end f01..f18 chain step: fetch surface + pressure GRIBs, decode, +//! build per-cell `Conditions`, score every band, write score-grid files +//! (`/.ntms`) to the shared scores directory. +//! +//! This is the Rust equivalent of `PropagationGridWorker.process_forecast_hour/4` +//! minus the f00-only enrichment (native duct, NEXRAD, commercial link +//! boost, ProfilesFile write) — Elixir retains f00. + +use std::collections::HashMap; +use std::path::Path; + +use chrono::{DateTime, Datelike, Timelike, Utc}; + +use crate::band_config; +use crate::decoder::{self, CellValues, PointGrid}; +use crate::fetcher::{self, HrrrClient, Product}; +use crate::grid::wgrib2_grid_spec; +use crate::scores_file::{self, ScorePoint}; +use crate::scorer::{self, Conditions}; +use crate::sounding_params::{self, Level}; + +#[derive(Debug, thiserror::Error)] +pub enum PipelineError { + #[error("fetch: {0}")] + Fetch(#[from] fetcher::FetchError), + #[error("decode: {0}")] + Decode(#[from] decoder::DecodeError), + #[error("write: {0}")] + Write(#[from] scores_file::WriteError), + #[error("forecast hour 0 is reserved for Elixir")] + F00Reserved, + #[error("surface grib was empty")] + EmptyGrib, +} + +#[derive(Debug, Clone)] +pub struct ChainStepInput { + pub run_time: DateTime, + pub forecast_hour: u8, +} + +impl ChainStepInput { + pub fn valid_time(&self) -> DateTime { + self.run_time + chrono::Duration::hours(self.forecast_hour as i64) + } +} + +#[derive(Debug, Clone)] +pub struct ChainStepStats { + pub score_files_written: u32, + pub point_count: u32, + pub band_count: u32, +} + +pub async fn run_chain_step( + client: &HrrrClient, + scores_dir: &Path, + step: &ChainStepInput, +) -> Result { + if step.forecast_hour == 0 { + return Err(PipelineError::F00Reserved); + } + let valid_time = step.valid_time(); + let date = step.run_time.date_naive(); + let hour = step.run_time.hour() as u8; + + let sfc_wanted = fetcher::surface_messages_owned(); + let prs_wanted = fetcher::pressure_messages_grid(); + let sfc_fut = client.fetch_product_blob( + date, + hour, + Product::Surface, + step.forecast_hour, + &sfc_wanted, + ); + let prs_fut = client.fetch_product_blob( + date, + hour, + Product::Pressure, + step.forecast_hour, + &prs_wanted, + ); + let (sfc_blob, prs_blob) = tokio::try_join!(sfc_fut, prs_fut)?; + if sfc_blob.is_empty() { + return Err(PipelineError::EmptyGrib); + } + + let grid_spec = wgrib2_grid_spec(); + let sfc_pattern = match_pattern(&sfc_wanted); + let prs_pattern = match_pattern(&prs_wanted); + + // Decoder runs wgrib2 subprocess — block so we don't hog the tokio + // reactor. `spawn_blocking` lifts it onto the dedicated pool. + let grid_spec_sfc = grid_spec; + let sfc_grid = tokio::task::spawn_blocking(move || { + decoder::extract_grid(&sfc_blob, &sfc_pattern, grid_spec_sfc) + }) + .await + .expect("blocking join")?; + + let grid_spec_prs = grid_spec; + let prs_grid = tokio::task::spawn_blocking(move || { + decoder::extract_grid(&prs_blob, &prs_pattern, grid_spec_prs) + }) + .await + .expect("blocking join")?; + + let merged = merge_grids(sfc_grid, prs_grid); + let point_count = merged.len() as u32; + + // Score every point × every band; one ScorePoint per (cell, band). + let bands = band_config::all_bands(); + let mut per_band: HashMap> = + HashMap::with_capacity(bands.len()); + + for (key, cell) in &merged { + let (lat, lon) = decoder::key_to_latlon(*key); + let conditions = match cell_to_conditions(cell, lat, lon, &valid_time) { + Some(c) => c, + None => continue, + }; + let invariants = scorer::precompute_band_invariants(&conditions); + for band in bands { + let r = scorer::composite_score_with(&conditions, band, Some(invariants)); + per_band + .entry(band.freq_mhz) + .or_default() + .push(ScorePoint { lat, lon, score: r.score }); + } + } + + let scores_dir = scores_dir.to_path_buf(); + let mut files_written = 0; + for (band_mhz, scores) in &per_band { + let dir = scores_dir.clone(); + let band_mhz = *band_mhz; + let scores = scores.clone(); + tokio::task::spawn_blocking(move || { + scores_file::write_atomic(&dir, band_mhz, valid_time, &scores) + }) + .await + .expect("blocking join")?; + files_written += 1; + } + + Ok(ChainStepStats { + score_files_written: files_written, + point_count, + band_count: bands.len() as u32, + }) +} + +/// wgrib2 `-match` regex: `":(A|B|C):"`. Var names only — levels are +/// post-filtered when we read the binary back. +fn match_pattern(wanted: &[(String, String)]) -> String { + let mut vars: Vec<&str> = wanted.iter().map(|(v, _)| v.as_str()).collect(); + vars.sort(); + vars.dedup(); + format!(":({}):", vars.join("|")) +} + +fn merge_grids(mut sfc: PointGrid, prs: PointGrid) -> PointGrid { + for (k, v) in prs { + sfc.entry(k).or_default().extend(v); + } + sfc +} + +fn cell_to_conditions( + cell: &CellValues, + lat: f64, + lon: f64, + valid_time: &DateTime, +) -> Option { + let tmp_k = cell.get("TMP:2 m above ground").copied()?; + let dpt_k = cell.get("DPT:2 m above ground").copied()?; + let temp_c = (tmp_k - 273.15) as f64; + let dewpoint_c = (dpt_k - 273.15) as f64; + let temp_f = scorer::c_to_f(temp_c); + let dewpoint_f = scorer::c_to_f(dewpoint_c); + let abs_humidity = scorer::absolute_humidity(temp_c, dewpoint_c); + + let wind_speed_kts = match ( + cell.get("UGRD:10 m above ground").copied(), + cell.get("VGRD:10 m above ground").copied(), + ) { + (Some(u), Some(v)) => Some(scorer::wind_speed_kts(u as f64, v as f64)), + _ => None, + }; + + let sky_cover_pct = cell.get("TCDC:entire atmosphere").copied().map(|v| v as f64); + let pwat_mm = cell + .get("PWAT:entire atmosphere (considered as a single layer)") + .copied() + .map(|v| v as f64); + let pressure_mb = cell + .get("PRES:surface") + .copied() + .map(|pa| pa as f64 / 100.0); + let rain_rate_mmhr = cell + .get("APCP:surface") + .copied() + .map(|mm| mm as f64) + .filter(|v| *v > 0.0); + let bl_depth_m = cell.get("HPBL:surface").copied().map(|v| v as f64); + + let levels: Vec = fetcher::GRID_PRESSURE_LEVELS + .iter() + .filter_map(|&p| { + let pres_mb = p as f64; + let t = cell.get(&format!("TMP:{p} mb")).copied()?; + let d = cell.get(&format!("DPT:{p} mb")).copied(); + let h = cell.get(&format!("HGT:{p} mb")).copied()?; + Some(Level { + pres_mb, + hght_m: h as f64, + tmpc: (t - 273.15) as f64, + dwpc: d.map(|v| (v - 273.15) as f64), + }) + }) + .collect(); + let min_refractivity_gradient = sounding_params::min_refractivity_gradient(levels); + + Some(Conditions { + abs_humidity, + temp_f, + dewpoint_f, + wind_speed_kts, + sky_cover_pct, + utc_hour: valid_time.hour() as u8, + utc_minute: valid_time.minute() as u8, + month: valid_time.month() as u8, + longitude: lon, + latitude: Some(lat), + pressure_mb, + prev_pressure_mb: None, + rain_rate_mmhr, + min_refractivity_gradient, + bl_depth_m, + pwat_mm, + best_duct_band_ghz: None, + bulk_richardson: None, + }) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn match_pattern_dedups_vars() { + let wanted = vec![ + ("TMP".into(), "2 m above ground".into()), + ("DPT".into(), "2 m above ground".into()), + ("TMP".into(), "1000 mb".into()), + ]; + assert_eq!(match_pattern(&wanted), ":(DPT|TMP):"); + } + + #[test] + fn cell_to_conditions_maps_fields() { + use chrono::TimeZone; + let mut cell = CellValues::new(); + cell.insert("TMP:2 m above ground".into(), 298.15); // 25 °C + cell.insert("DPT:2 m above ground".into(), 293.15); // 20 °C + cell.insert("UGRD:10 m above ground".into(), 3.0); + cell.insert("VGRD:10 m above ground".into(), 4.0); + cell.insert("TCDC:entire atmosphere".into(), 30.0); + cell.insert("PRES:surface".into(), 101_000.0); + cell.insert("HPBL:surface".into(), 400.0); + cell.insert("PWAT:entire atmosphere (considered as a single layer)".into(), 25.0); + + let vt = Utc.with_ymd_and_hms(2026, 6, 15, 18, 0, 0).unwrap(); + let c = cell_to_conditions(&cell, 32.0, -97.0, &vt).unwrap(); + assert!((c.temp_f - 77.0).abs() < 0.1); + assert!((c.dewpoint_f - 68.0).abs() < 0.1); + assert!((c.wind_speed_kts.unwrap() - 9.72).abs() < 0.1); + assert!((c.pressure_mb.unwrap() - 1010.0).abs() < 0.01); + assert_eq!(c.sky_cover_pct, Some(30.0)); + assert_eq!(c.bl_depth_m, Some(400.0)); + assert_eq!(c.pwat_mm, Some(25.0)); + assert_eq!(c.month, 6); + assert_eq!(c.longitude, -97.0); + } + + #[test] + fn cell_missing_surface_returns_none() { + use chrono::TimeZone; + let cell = CellValues::new(); + let vt = Utc.with_ymd_and_hms(2026, 6, 15, 18, 0, 0).unwrap(); + assert!(cell_to_conditions(&cell, 32.0, -97.0, &vt).is_none()); + } + + #[test] + fn rejects_forecast_hour_zero() { + use chrono::TimeZone; + let client = HrrrClient::default_base().unwrap(); + let dir = tempfile::tempdir().unwrap(); + let step = ChainStepInput { + run_time: Utc.with_ymd_and_hms(2026, 4, 19, 15, 0, 0).unwrap(), + forecast_hour: 0, + }; + let rt = tokio::runtime::Runtime::new().unwrap(); + let err = rt.block_on(run_chain_step(&client, dir.path(), &step)).unwrap_err(); + assert!(matches!(err, PipelineError::F00Reserved)); + } +} diff --git a/rust/prop_grid_rs/src/region.rs b/rust/prop_grid_rs/src/region.rs new file mode 100644 index 00000000..95025ba9 --- /dev/null +++ b/rust/prop_grid_rs/src/region.rs @@ -0,0 +1,105 @@ +//! Climatological region lookup. 1:1 port of +//! `lib/microwaveprop/propagation/region.ex`. + +#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] +pub enum Region { + GulfCoast, + Southeast, + SouthernPlains, + CornBelt, + Northeast, + DesertSouthwest, + PacificNorthwest, + MountainWest, + Other, +} + +/// Lat min (inclusive), lat max (exclusive), lon min (inclusive), lon max (exclusive). +struct BBox(Region, f64, f64, f64, f64); + +// Ordered to match Elixir's `Enum.find_value` short-circuit. +const REGIONS: &[BBox] = &[ + BBox(Region::GulfCoast, 25.0, 32.0, -100.0, -80.0), + BBox(Region::Southeast, 30.0, 37.0, -90.0, -75.0), + BBox(Region::SouthernPlains, 32.0, 38.0, -105.0, -93.0), + BBox(Region::CornBelt, 38.0, 48.0, -100.0, -82.0), + BBox(Region::Northeast, 38.0, 48.0, -82.0, -67.0), + BBox(Region::DesertSouthwest, 30.0, 38.0, -120.0, -105.0), + BBox(Region::PacificNorthwest, 42.0, 50.0, -125.0, -115.0), + BBox(Region::MountainWest, 38.0, 48.0, -115.0, -100.0), +]; + +pub fn for_point(lat: f64, lon: f64) -> Region { + for bbox in REGIONS { + let BBox(name, lat_min, lat_max, lon_min, lon_max) = *bbox; + if lat >= lat_min && lat < lat_max && lon >= lon_min && lon < lon_max { + return name; + } + } + Region::Other +} + +/// Multiplier applied to the band's `seasonal_base` score. Unknown +/// region/month combinations return 1.0. Range [0.7, 1.3]. +pub fn seasonal_adjustment(region: Region, month: u8) -> f64 { + match (region, month) { + (Region::GulfCoast, 6) => 0.95, + (Region::GulfCoast, 7) => 0.95, + (Region::GulfCoast, 8) => 1.15, + (Region::GulfCoast, 9) => 1.10, + (Region::CornBelt, 6) => 1.05, + (Region::CornBelt, 7) => 0.85, + (Region::CornBelt, 8) => 0.80, + (Region::CornBelt, 9) => 0.95, + (Region::Southeast, 7) => 0.97, + (Region::Southeast, 8) => 1.08, + (Region::SouthernPlains, 8) => 1.05, + (Region::DesertSouthwest, 7) => 1.10, + (Region::DesertSouthwest, 8) => 1.10, + (Region::DesertSouthwest, 9) => 1.05, + (Region::PacificNorthwest, 6) => 1.15, + (Region::PacificNorthwest, 7) => 1.20, + (Region::PacificNorthwest, 8) => 1.20, + (Region::PacificNorthwest, 9) => 1.10, + _ => 1.0, + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn gulf_coast_points() { + // Houston ≈ (29.8, -95.4) — Gulf Coast bbox. + assert_eq!(for_point(29.8, -95.4), Region::GulfCoast); + } + + #[test] + fn corn_belt_points() { + // Des Moines ≈ (41.6, -93.6) — Corn Belt bbox. + assert_eq!(for_point(41.6, -93.6), Region::CornBelt); + } + + #[test] + fn non_conus_returns_other() { + // Alaska (61.2, -149.9) — outside every bbox. + assert_eq!(for_point(61.2, -149.9), Region::Other); + } + + #[test] + fn adjustments_match_elixir_table() { + assert_eq!(seasonal_adjustment(Region::GulfCoast, 8), 1.15); + assert_eq!(seasonal_adjustment(Region::CornBelt, 7), 0.85); + assert_eq!(seasonal_adjustment(Region::PacificNorthwest, 7), 1.20); + assert_eq!(seasonal_adjustment(Region::Other, 8), 1.0); + assert_eq!(seasonal_adjustment(Region::CornBelt, 1), 1.0); + } + + #[test] + fn precedence_matches_elixir_ordering() { + // (30.5, -90.0) is inside BOTH gulf_coast (25..32, -100..-80) and + // southeast (30..37, -90..-75). Elixir hits gulf_coast first. + assert_eq!(for_point(30.5, -90.0), Region::GulfCoast); + } +} diff --git a/rust/prop_grid_rs/src/scorer.rs b/rust/prop_grid_rs/src/scorer.rs new file mode 100644 index 00000000..afc364ac --- /dev/null +++ b/rust/prop_grid_rs/src/scorer.rs @@ -0,0 +1,949 @@ +//! Propagation scorer — 1:1 port of `lib/microwaveprop/propagation/scorer.ex`. +//! +//! Every factor returns a 0..=100 integer score. The composite is +//! `round(Σ factor · weight)` with weights from `band_config`. All +//! thresholds live in `band_config` / `region` — nothing magic here. +//! +//! Behavior parity tests mirror the Elixir `ScorerTest` suite line-by-line; +//! integer scores match exactly. +//! +//! Elixir `Kernel.round/1` and Rust `f64::round` both use round-half-away- +//! from-zero, so no custom rounding helper is needed. + +use crate::band_config::{ + BandConfig, HumidityEffect, Weights, DEFAULT_WEIGHTS, + HUMIDITY_BENEFICIAL_DEFAULT, HUMIDITY_BENEFICIAL_THRESHOLDS, + REFRACTIVITY_DEFAULT, REFRACTIVITY_THRESHOLDS, SUNRISE_TABLE, +}; +use crate::region; + +const BULK_RICHARDSON_STABLE_MAX: f64 = 25.0; + +/// Inputs to `composite_score`. All f64 everywhere to stay close to +/// Erlang's single-float world. +#[derive(Debug, Clone, Default)] +pub struct Conditions { + pub abs_humidity: f64, + pub temp_f: f64, + pub dewpoint_f: f64, + pub wind_speed_kts: Option, + pub sky_cover_pct: Option, + pub utc_hour: u8, + pub utc_minute: u8, + pub month: u8, + pub longitude: f64, + pub latitude: Option, + pub pressure_mb: Option, + pub prev_pressure_mb: Option, + pub rain_rate_mmhr: Option, + pub min_refractivity_gradient: Option, + pub bl_depth_m: Option, + pub pwat_mm: Option, + pub best_duct_band_ghz: Option, + pub bulk_richardson: Option, +} + +#[derive(Debug, Clone, Copy)] +pub struct BandInvariants { + pub tod: i32, + pub sky: i32, + pub wind: i32, + pub pressure: i32, +} + +#[derive(Debug, Clone)] +pub struct Factors { + pub humidity: i32, + pub time_of_day: i32, + pub td_depression: i32, + pub refractivity: i32, + pub sky: i32, + pub season: i32, + pub wind: i32, + pub rain: i32, + pub pwat: i32, + pub pressure: i32, +} + +#[derive(Debug, Clone)] +pub struct Composite { + pub score: i32, + pub factors: Factors, +} + +// ── Helpers ────────────────────────────────────────────────────────── + +pub fn f_to_c(f: f64) -> f64 { + (f - 32.0) * 5.0 / 9.0 +} + +pub fn c_to_f(c: f64) -> f64 { + c * 9.0 / 5.0 + 32.0 +} + +pub fn absolute_humidity(temp_c: f64, dewpoint_c: f64) -> f64 { + let e_sat = 6.112 * (17.67 * dewpoint_c / (dewpoint_c + 243.5)).exp(); + 217.0 * e_sat / (temp_c + 273.15) +} + +pub fn wind_speed_kts(u_ms: f64, v_ms: f64) -> f64 { + (u_ms * u_ms + v_ms * v_ms).sqrt() * 1.94384 +} + +pub fn precip_to_rate_mmhr(mm: Option) -> f64 { + match mm { + Some(v) if v > 0.0 => v, + _ => 0.0, + } +} + +/// Marshall-Palmer Z-R: `Z = 200 · R^1.6` → `R = (Z/200)^(1/1.6)`, `Z = 10^(dBZ/10)`. +/// Clips below 5 dBZ to 0 (ground clutter) and above 150 mm/hr (hail). +pub fn dbz_to_rain_rate_mmhr(dbz: Option) -> f64 { + match dbz { + None => 0.0, + Some(v) if v < 5.0 => 0.0, + Some(v) => { + let z = 10f64.powf(v / 10.0); + (z / 200.0).powf(1.0 / 1.6).min(150.0) + } + } +} + +fn clamp_score_f(v: f64) -> i32 { + (v.round() as i64).clamp(0, 100) as i32 +} + +fn clamp_score(v: i32) -> i32 { + v.clamp(0, 100) +} + +// ── Factor 1: humidity ─────────────────────────────────────────────── + +pub fn score_humidity(abs_humidity: f64, band: &BandConfig) -> i32 { + match band.humidity_effect { + HumidityEffect::Beneficial => HUMIDITY_BENEFICIAL_THRESHOLDS + .iter() + .find(|(max, _)| abs_humidity < *max as f64) + .map(|(_, s)| *s) + .unwrap_or(HUMIDITY_BENEFICIAL_DEFAULT), + HumidityEffect::Harmful => { + let r = abs_humidity * band.humidity_penalty; + if r <= 6.0 { + 100 + } else if r <= 9.0 { + (95.0 - (r - 6.0) / 3.0 * 20.0).round() as i32 + } else if r <= 13.0 { + (75.0 - (r - 9.0) / 4.0 * 30.0).round() as i32 + } else if r <= 18.0 { + (45.0 - (r - 13.0) / 5.0 * 35.0).round() as i32 + } else { + ((10.0 - (r - 18.0) * 2.0).round() as i32).max(0) + } + } + } +} + +// ── Factor 2: time of day ──────────────────────────────────────────── + +pub fn score_time_of_day( + utc_hour: u8, + utc_minute: u8, + month: u8, + longitude: f64, +) -> (i32, &'static str) { + let offset = longitude / 15.0; + let raw = utc_hour as f64 + utc_minute as f64 / 60.0 + offset + 24.0; + // Match Elixir :math.fmod/2 — same semantics as fmod. + let local = raw.rem_euclid(24.0); + let sunrise = SUNRISE_TABLE[(month - 1) as usize]; + let d = local - sunrise; + + if (-1.5..=1.5).contains(&d) { + (100, "Peak — inversion maximum") + } else if d > 1.5 && d <= 3.0 { + (78, "Good — inversion eroding") + } else if d > -3.0 && d < -1.5 { + (82, "Pre-dawn — inversion building") + } else if d > 3.0 && d <= 6.0 { + (38, "Marginal — boundary layer mixing") + } else if local >= 20.0 || local <= 1.0 { + (72, "Evening — cooling, inversion reforming") + } else if d > 6.0 { + (18, "Afternoon — full convective mixing") + } else { + (55, "Night — gradual cooling") + } +} + +// ── Factor 3: T-Td depression ──────────────────────────────────────── + +pub fn score_td_depression(temp_f: f64, dewpoint_f: f64, band: &BandConfig) -> i32 { + let dep = temp_f - dewpoint_f; + match band.humidity_effect { + HumidityEffect::Beneficial => { + if dep < 3.0 { + 40 + } else if dep < 8.0 { + 75 + } else if dep < 14.0 { + 85 + } else if dep < 22.0 { + 70 + } else { + 55 + } + } + HumidityEffect::Harmful => { + if dep > 22.0 { + 96 + } else if dep > 14.0 { + 80 + } else if dep > 8.0 { + 60 + } else if dep > 4.0 { + 38 + } else { + 18 + } + } + } +} + +// ── Factor 4: refractivity gradient (+ HPBL multiplier + native duct boost) ── + +pub fn score_refractivity( + min_gradient: Option, + bl_depth_m: Option, + best_duct_band_ghz: Option, + bulk_richardson: Option, + band: &BandConfig, +) -> i32 { + let Some(gradient) = min_gradient else { + return 50; + }; + + let base_match = REFRACTIVITY_THRESHOLDS + .iter() + .find(|(max, _, _)| gradient < *max as f64); + + let base = match base_match { + Some((_, b, h)) => match band.humidity_effect { + HumidityEffect::Beneficial => *b, + HumidityEffect::Harmful => *h, + }, + None => { + let (b, h) = REFRACTIVITY_DEFAULT; + if band.humidity_effect == HumidityEffect::Beneficial { + b + } else { + h + } + } + }; + + let with_hpbl = apply_hpbl_multiplier(base, bl_depth_m); + apply_native_duct_boost(with_hpbl, best_duct_band_ghz, bulk_richardson, band.freq_mhz) +} + +fn apply_hpbl_multiplier(base: i32, hpbl: Option) -> i32 { + let mul = match hpbl { + None => return base, + Some(v) if v < 200.0 => 1.10, + Some(v) if v < 500.0 => 1.05, + Some(v) if v < 1500.0 => 1.0, + Some(v) if v < 2000.0 => 0.92, + Some(_) => 0.78, + }; + clamp_score_f(base as f64 * mul) +} + +fn apply_native_duct_boost( + base: i32, + best_duct_band_ghz: Option, + bulk_richardson: Option, + freq_mhz: u32, +) -> i32 { + let Some(duct_band_ghz) = best_duct_band_ghz else { + return base; + }; + let target_ghz = freq_mhz as f64 / 1000.0; + let stable = match bulk_richardson { + None => true, + Some(r) => r < BULK_RICHARDSON_STABLE_MAX, + }; + if duct_band_ghz >= target_ghz && stable { + clamp_score_f(base as f64 * 1.15) + } else { + base + } +} + +// ── Factor 5: sky cover ────────────────────────────────────────────── + +pub fn score_sky(pct: Option) -> i32 { + match pct { + None => 50, + Some(p) if p <= 6.0 => 100, + Some(p) if p <= 25.0 => 88, + Some(p) if p <= 50.0 => 60, + Some(p) if p <= 87.0 => 25, + Some(_) => 5, + } +} + +// ── Factor 6: season ───────────────────────────────────────────────── + +pub fn score_season(month: u8, lat: Option, lon: Option, band: &BandConfig) -> i32 { + let idx = (month - 1) as usize; + let base = band.seasonal_base[idx]; + let adj = band.seasonal_adj[idx]; + + let region_mult = match (lat, lon) { + (Some(la), Some(lo)) => region::seasonal_adjustment(region::for_point(la, lo), month), + _ => 1.0, + }; + + let raw = (base + adj) as f64 * region_mult; + clamp_score_f(raw) +} + +// ── Factor 7: wind ─────────────────────────────────────────────────── + +pub fn score_wind(speed_kts: Option) -> i32 { + match speed_kts { + None => 50, + Some(s) if s < 5.0 => 100, + Some(s) if s < 10.0 => 90, + Some(s) if s < 15.0 => 75, + Some(s) if s < 20.0 => 55, + Some(s) if s < 25.0 => 35, + Some(_) => 15, + } +} + +// ── Factor 8: rain attenuation (ITU-R P.838) ──────────────────────── + +pub fn score_rain(rate_mmhr: Option, band: &BandConfig) -> i32 { + match rate_mmhr { + None => 100, + Some(0.0) => 100, + Some(v) => { + let gamma = band.rain_k * v.powf(band.rain_alpha); + if gamma < 0.1 { + 95 + } else if gamma < 0.5 { + 75 + } else if gamma < 1.0 { + 50 + } else if gamma < 2.0 { + 25 + } else if gamma < 5.0 { + 10 + } else { + 0 + } + } + } +} + +// ── Factor 9: precipitable water ──────────────────────────────────── + +pub fn score_pwat(pwat_mm: Option, band: &BandConfig) -> i32 { + let Some(v) = pwat_mm else { return 60 }; + match band.humidity_effect { + HumidityEffect::Beneficial => { + if v < 10.0 { + 55 + } else if v < 20.0 { + 75 + } else if v < 30.0 { + 90 + } else if v < 40.0 { + 70 + } else { + 50 + } + } + HumidityEffect::Harmful => { + if v < 10.0 { + 95 + } else if v < 20.0 { + 80 + } else if v < 30.0 { + 60 + } else if v < 40.0 { + 35 + } else { + 15 + } + } + } +} + +// ── Factor 10: pressure ───────────────────────────────────────────── + +pub fn score_pressure(current_mb: Option, previous_mb: Option) -> i32 { + let Some(current) = current_mb else { return 50 }; + match previous_mb { + None => { + if current < 980.0 { + 88 + } else if current < 990.0 { + 82 + } else if current < 1000.0 { + 70 + } else if current < 1010.0 { + 55 + } else if current < 1020.0 { + 40 + } else { + 30 + } + } + Some(prev) => { + let delta = current - prev; + if delta > 2.5 { + 80 + } else if delta > 0.8 { + 70 + } else if delta > -0.5 { + 60 + } else if delta > -2.0 { + 65 + } else { + 45 + } + } + } +} + +// ── Composite ─────────────────────────────────────────────────────── + +pub fn precompute_band_invariants(c: &Conditions) -> BandInvariants { + let (tod, _) = score_time_of_day(c.utc_hour, c.utc_minute, c.month, c.longitude); + BandInvariants { + tod, + sky: score_sky(c.sky_cover_pct), + wind: score_wind(c.wind_speed_kts), + pressure: score_pressure(c.pressure_mb, c.prev_pressure_mb), + } +} + +pub fn composite_score(c: &Conditions, band: &BandConfig) -> Composite { + composite_score_with(c, band, None) +} + +pub fn composite_score_with( + c: &Conditions, + band: &BandConfig, + precomputed: Option, +) -> Composite { + let inv = precomputed.unwrap_or_else(|| precompute_band_invariants(c)); + + let factors = Factors { + humidity: score_humidity(c.abs_humidity, band), + time_of_day: inv.tod, + td_depression: score_td_depression(c.temp_f, c.dewpoint_f, band), + refractivity: score_refractivity( + c.min_refractivity_gradient, + c.bl_depth_m, + c.best_duct_band_ghz, + c.bulk_richardson, + band, + ), + sky: inv.sky, + season: score_season(c.month, c.latitude, Some(c.longitude).filter(|_| c.latitude.is_some()), band), + wind: inv.wind, + rain: score_rain(c.rain_rate_mmhr, band), + pwat: score_pwat(c.pwat_mm, band), + pressure: inv.pressure, + }; + + let w = band.weights.unwrap_or(DEFAULT_WEIGHTS); + let weighted = weighted_sum(&factors, &w); + Composite { + score: clamp_score_f(weighted), + factors, + } +} + +fn weighted_sum(f: &Factors, w: &Weights) -> f64 { + f.humidity as f64 * w.humidity + + f.time_of_day as f64 * w.time_of_day + + f.td_depression as f64 * w.td_depression + + f.refractivity as f64 * w.refractivity + + f.sky as f64 * w.sky + + f.season as f64 * w.season + + f.wind as f64 * w.wind + + f.rain as f64 * w.rain + + f.pressure as f64 * w.pressure + + f.pwat as f64 * w.pwat +} + +/// Commercial-link inverse-sensor boost. Not used in the grid hot path +/// (f00-only — Elixir keeps that) but ported so unit-level correctness +/// stays testable from the Rust side. +pub fn commercial_link_boost(score: i32, n_links: u32, degradation_db: f64) -> i32 { + if n_links == 0 || degradation_db < 3.0 { + return clamp_score(score); + } + if degradation_db < 8.0 { + let bonus = 2.0 + (degradation_db - 3.0) * 8.0 / 5.0; + clamp_score_f(score as f64 + bonus) + } else { + let bonus = (10.0 + (degradation_db - 8.0) * 15.0 / 8.0).min(25.0); + clamp_score_f(score as f64 + bonus) + } +} + +// ── Tests (mirror Elixir ScorerTest) ──────────────────────────────── + +#[cfg(test)] +mod tests { + use super::*; + use crate::band_config; + + fn b10g() -> &'static BandConfig { + band_config::get(10_000).unwrap() + } + fn b24g() -> &'static BandConfig { + band_config::get(24_000).unwrap() + } + fn b75g() -> &'static BandConfig { + band_config::get(75_000).unwrap() + } + + const EAST_LON: f64 = -75.0; + + #[test] + fn f_to_c_boundaries() { + assert!((f_to_c(32.0) - 0.0).abs() < 0.01); + assert!((f_to_c(212.0) - 100.0).abs() < 0.01); + } + + #[test] + fn c_to_f_boundaries() { + assert!((c_to_f(0.0) - 32.0).abs() < 0.01); + assert!((c_to_f(100.0) - 212.0).abs() < 0.01); + } + + #[test] + fn absolute_humidity_typical_summer() { + let ah = absolute_humidity(25.0, 20.0); + assert!(ah > 15.0 && ah < 20.0, "{ah}"); + } + + #[test] + fn absolute_humidity_cold_dry() { + let ah = absolute_humidity(0.0, -5.0); + assert!(ah > 2.0 && ah < 5.0, "{ah}"); + } + + #[test] + fn wind_speed_3_4_kts() { + assert!((wind_speed_kts(3.0, 4.0) - 9.72).abs() < 0.1); + } + + #[test] + fn precip_to_rate_handles_nil_zero_negative() { + assert_eq!(precip_to_rate_mmhr(Some(2.5)), 2.5); + assert_eq!(precip_to_rate_mmhr(Some(0.0)), 0.0); + assert_eq!(precip_to_rate_mmhr(Some(-1.0)), 0.0); + assert_eq!(precip_to_rate_mmhr(None), 0.0); + } + + #[test] + fn dbz_to_rain_rate_clips_and_interpolates() { + assert_eq!(dbz_to_rain_rate_mmhr(None), 0.0); + assert_eq!(dbz_to_rain_rate_mmhr(Some(3.0)), 0.0); + assert!((dbz_to_rain_rate_mmhr(Some(20.0)) - 0.66).abs() < 0.05); + assert!((dbz_to_rain_rate_mmhr(Some(30.0)) - 2.7).abs() < 0.2); + assert!((dbz_to_rain_rate_mmhr(Some(45.0)) - 23.7).abs() < 1.0); + assert!(dbz_to_rain_rate_mmhr(Some(55.0)) <= 150.0); + assert!(dbz_to_rain_rate_mmhr(Some(65.0)) <= 150.0); + } + + // ── humidity ───────────────────────────────────────────────── + + #[test] + fn humidity_beneficial_table() { + assert_eq!(score_humidity(3.0, b10g()), 55); + assert_eq!(score_humidity(9.0, b10g()), 82); + assert_eq!(score_humidity(16.0, b10g()), 95); + assert_eq!(score_humidity(20.0, b10g()), 88); + assert_eq!(score_humidity(25.0, b10g()), 75); + } + + #[test] + fn humidity_harmful_piecewise() { + assert_eq!(score_humidity(2.0, b24g()), 100); + assert_eq!(score_humidity(5.0, b24g()), 82); + assert_eq!(score_humidity(10.0, b24g()), 24); + assert_eq!(score_humidity(20.0, b24g()), 0); + } + + // ── time of day ────────────────────────────────────────────── + + #[test] + fn time_of_day_dawn_peak() { + let (s, l) = score_time_of_day(11, 15, 6, EAST_LON); + assert_eq!(s, 100); + assert!(l.contains("Peak")); + } + + #[test] + fn time_of_day_afternoon() { + let (s, _) = score_time_of_day(22, 0, 6, EAST_LON); + assert_eq!(s, 18); + } + + #[test] + fn time_of_day_pre_dawn() { + let (s, l) = score_time_of_day(9, 0, 6, EAST_LON); + assert_eq!(s, 82); + assert!(l.contains("Pre-dawn")); + } + + #[test] + fn time_of_day_evening() { + let (s, l) = score_time_of_day(1, 0, 6, EAST_LON); + assert_eq!(s, 72); + assert!(l.contains("Evening")); + } + + #[test] + fn time_of_day_western_longitude_shifts_earlier() { + let (s, _) = score_time_of_day(13, 24, 1, -90.0); + assert_eq!(s, 100); + } + + #[test] + fn time_of_day_same_utc_different_longitudes() { + let (east, _) = score_time_of_day(18, 0, 6, -75.0); + let (west, _) = score_time_of_day(18, 0, 6, -120.0); + assert!(west > east); + } + + // ── td depression ──────────────────────────────────────────── + + #[test] + fn td_beneficial() { + assert_eq!(score_td_depression(72.0, 70.0, b10g()), 40); + assert_eq!(score_td_depression(80.0, 70.0, b10g()), 85); + assert_eq!(score_td_depression(100.0, 50.0, b10g()), 55); + } + + #[test] + fn td_harmful() { + assert_eq!(score_td_depression(72.0, 70.0, b24g()), 18); + assert_eq!(score_td_depression(80.0, 70.0, b24g()), 60); + assert_eq!(score_td_depression(100.0, 50.0, b24g()), 96); + } + + // ── refractivity ───────────────────────────────────────────── + + const BASELINE_BL: Option = Some(750.0); + + #[test] + fn refractivity_strong_gradient_beneficial() { + assert_eq!(score_refractivity(Some(-600.0), BASELINE_BL, None, None, b10g()), 98); + } + + #[test] + fn refractivity_strong_gradient_harmful() { + assert_eq!(score_refractivity(Some(-250.0), BASELINE_BL, None, None, b24g()), 85); + } + + #[test] + fn refractivity_moderate_gradient_beneficial() { + assert_eq!(score_refractivity(Some(-160.0), BASELINE_BL, None, None, b10g()), 92); + } + + #[test] + fn refractivity_nil_gradient_returns_50() { + assert_eq!(score_refractivity(None, BASELINE_BL, None, None, b10g()), 50); + } + + #[test] + fn refractivity_nil_bl_no_multiplier() { + assert_eq!(score_refractivity(Some(-160.0), None, None, None, b10g()), 92); + } + + #[test] + fn shallow_bl_boosts_score() { + let shallow = score_refractivity(Some(-100.0), Some(150.0), None, None, b10g()); + let baseline = score_refractivity(Some(-100.0), BASELINE_BL, None, None, b10g()); + assert!(shallow > baseline); + assert!(shallow <= 100); + } + + #[test] + fn deep_bl_penalises_score() { + let deep = score_refractivity(Some(-100.0), Some(2200.0), None, None, b10g()); + let baseline = score_refractivity(Some(-100.0), BASELINE_BL, None, None, b10g()); + assert!(deep < baseline); + assert!(deep >= 0); + } + + #[test] + fn shallow_bl_also_lifts_default_fallback() { + let shallow = score_refractivity(Some(-30.0), Some(150.0), None, None, b10g()); + let baseline = score_refractivity(Some(-30.0), BASELINE_BL, None, None, b10g()); + assert!(shallow > baseline); + } + + #[test] + fn native_duct_boosts_matching_band() { + let boosted = score_refractivity(Some(-80.0), Some(800.0), Some(15.0), None, b10g()); + let plain = score_refractivity(Some(-80.0), Some(800.0), None, None, b10g()); + assert!(boosted > plain); + assert!(boosted <= 100); + } + + #[test] + fn native_duct_below_target_no_boost() { + let unchanged = score_refractivity(Some(-80.0), Some(800.0), Some(3.0), None, b10g()); + let plain = score_refractivity(Some(-80.0), Some(800.0), None, None, b10g()); + assert_eq!(unchanged, plain); + } + + #[test] + fn richardson_gates_the_boost() { + let stable = score_refractivity(Some(-80.0), Some(800.0), Some(15.0), Some(12.0), b10g()); + let plain = score_refractivity(Some(-80.0), Some(800.0), None, None, b10g()); + assert!(stable > plain); + + let turbulent = + score_refractivity(Some(-80.0), Some(800.0), Some(15.0), Some(40.0), b10g()); + assert_eq!(turbulent, plain); + + let nil_r = score_refractivity(Some(-80.0), Some(800.0), Some(15.0), None, b10g()); + let four_arity_equiv = + score_refractivity(Some(-80.0), Some(800.0), Some(15.0), None, b10g()); + assert_eq!(nil_r, four_arity_equiv); + + let boundary = + score_refractivity(Some(-80.0), Some(800.0), Some(15.0), Some(25.0), b10g()); + assert_eq!(boundary, plain); + } + + // ── commercial link boost ──────────────────────────────────── + + #[test] + fn commercial_link_noop_cases() { + assert_eq!(commercial_link_boost(75, 0, 10.0), 75); + assert_eq!(commercial_link_boost(75, 3, 2.0), 75); + } + + #[test] + fn commercial_link_mild_boost() { + let b = commercial_link_boost(75, 3, 5.0); + assert!(b > 75 && b < 90, "{b}"); + } + + #[test] + fn commercial_link_strong_boost_capped() { + let b = commercial_link_boost(75, 3, 12.0); + assert!(b > 85 && b <= 100, "{b}"); + } + + // ── sky ────────────────────────────────────────────────────── + + #[test] + fn sky_table() { + assert_eq!(score_sky(Some(5.0)), 100); + assert_eq!(score_sky(Some(20.0)), 88); + assert_eq!(score_sky(Some(40.0)), 60); + assert_eq!(score_sky(Some(75.0)), 25); + assert_eq!(score_sky(Some(95.0)), 5); + assert_eq!(score_sky(None), 50); + } + + // ── season ─────────────────────────────────────────────────── + + #[test] + fn season_table_basic() { + assert_eq!(score_season(7, None, None, b10g()), 95); + assert_eq!(score_season(3, None, None, b10g()), 22); + assert_eq!(score_season(11, None, None, b24g()), 96); + assert_eq!(score_season(7, None, None, b24g()), 8); + } + + #[test] + fn season_gulf_coast_august_boost() { + let no_region = score_season(8, None, None, b10g()); + let gulf = score_season(8, Some(29.0), Some(-95.0), b10g()); + assert!(gulf > no_region); + } + + #[test] + fn season_corn_belt_august_penalty() { + let no_region = score_season(8, None, None, b10g()); + let iowa = score_season(8, Some(42.0), Some(-93.0), b10g()); + assert!(iowa < no_region); + } + + // ── wind ───────────────────────────────────────────────────── + + #[test] + fn wind_table() { + assert_eq!(score_wind(Some(3.0)), 100); + assert_eq!(score_wind(Some(8.0)), 90); + assert_eq!(score_wind(Some(12.0)), 75); + assert_eq!(score_wind(Some(22.0)), 35); + assert_eq!(score_wind(Some(30.0)), 15); + assert_eq!(score_wind(None), 50); + } + + // ── rain ───────────────────────────────────────────────────── + + #[test] + fn rain_zero_and_nil_full_score() { + assert_eq!(score_rain(Some(0.0), b24g()), 100); + assert_eq!(score_rain(None, b24g()), 100); + } + + #[test] + fn rain_24g_light() { + // gamma = 0.070 * 2.0^1.07 ≈ 0.147 → 75 + assert_eq!(score_rain(Some(2.0), b24g()), 75); + } + + #[test] + fn rain_10g_heavy_still_high() { + // gamma = 0.010 * 5.0^1.28 ≈ 0.076 → 95 + assert_eq!(score_rain(Some(5.0), b10g()), 95); + } + + #[test] + fn rain_75g_heavy_drops() { + // gamma = 0.345 * 10^0.84 ≈ 2.39 → 10 + assert_eq!(score_rain(Some(10.0), b75g()), 10); + } + + // ── pwat ───────────────────────────────────────────────────── + + #[test] + fn pwat_beneficial_table() { + assert_eq!(score_pwat(Some(5.0), b10g()), 55); + assert_eq!(score_pwat(Some(15.0), b10g()), 75); + assert_eq!(score_pwat(Some(25.0), b10g()), 90); + assert_eq!(score_pwat(Some(35.0), b10g()), 70); + assert_eq!(score_pwat(Some(45.0), b10g()), 50); + } + + #[test] + fn pwat_harmful_table() { + assert_eq!(score_pwat(Some(5.0), b24g()), 95); + assert_eq!(score_pwat(Some(15.0), b24g()), 80); + assert_eq!(score_pwat(Some(25.0), b24g()), 60); + assert_eq!(score_pwat(Some(35.0), b24g()), 35); + assert_eq!(score_pwat(Some(45.0), b24g()), 15); + } + + #[test] + fn pwat_nil_returns_60() { + assert_eq!(score_pwat(None, b10g()), 60); + assert_eq!(score_pwat(None, b24g()), 60); + } + + // ── pressure ───────────────────────────────────────────────── + + #[test] + fn pressure_standalone_table() { + assert_eq!(score_pressure(Some(1028.0), None), 30); + assert_eq!(score_pressure(Some(1020.0), None), 30); + assert_eq!(score_pressure(Some(1000.0), None), 55); + } + + #[test] + fn pressure_deltas() { + assert_eq!(score_pressure(Some(1020.0), Some(1015.0)), 80); + assert_eq!(score_pressure(Some(1016.0), Some(1015.0)), 70); + assert_eq!(score_pressure(Some(1015.0), Some(1015.3)), 60); + assert_eq!(score_pressure(Some(1014.0), Some(1015.0)), 65); + assert_eq!(score_pressure(Some(1010.0), Some(1015.0)), 45); + } + + // ── composite ──────────────────────────────────────────────── + + fn canonical_conditions() -> Conditions { + Conditions { + abs_humidity: 10.0, + temp_f: 80.0, + dewpoint_f: 70.0, + wind_speed_kts: Some(5.0), + sky_cover_pct: Some(20.0), + utc_hour: 11, + utc_minute: 15, + month: 6, + longitude: -75.0, + latitude: None, + pressure_mb: Some(1020.0), + prev_pressure_mb: Some(1018.0), + rain_rate_mmhr: Some(0.0), + min_refractivity_gradient: Some(-350.0), + bl_depth_m: Some(400.0), + pwat_mm: Some(25.0), + best_duct_band_ghz: None, + bulk_richardson: None, + } + } + + #[test] + fn composite_returns_0_100_score() { + let r = composite_score(&canonical_conditions(), b10g()); + assert!(r.score >= 0 && r.score <= 100); + } + + #[test] + fn composite_all_factors_populated() { + let r = composite_score(&canonical_conditions(), b10g()); + let f = &r.factors; + for s in [ + f.humidity, f.time_of_day, f.td_depression, f.refractivity, + f.sky, f.season, f.wind, f.rain, f.pwat, f.pressure, + ] { + assert!((0..=100).contains(&s), "factor out of range: {s}"); + } + } + + #[test] + fn precomputed_invariants_match_fresh() { + let c = canonical_conditions(); + let inv = precompute_band_invariants(&c); + let fresh = composite_score(&c, b10g()); + let reused = composite_score_with(&c, b10g(), Some(inv)); + assert_eq!(fresh.factors.time_of_day, reused.factors.time_of_day); + assert_eq!(fresh.factors.sky, reused.factors.sky); + assert_eq!(fresh.factors.wind, reused.factors.wind); + assert_eq!(fresh.factors.pressure, reused.factors.pressure); + assert_eq!(fresh.score, reused.score); + } + + #[test] + fn composite_uses_band_weights() { + let r = composite_score(&canonical_conditions(), b10g()); + let w = DEFAULT_WEIGHTS; + let expected = r.factors.humidity as f64 * w.humidity + + r.factors.time_of_day as f64 * w.time_of_day + + r.factors.td_depression as f64 * w.td_depression + + r.factors.refractivity as f64 * w.refractivity + + r.factors.sky as f64 * w.sky + + r.factors.season as f64 * w.season + + r.factors.wind as f64 * w.wind + + r.factors.rain as f64 * w.rain + + r.factors.pwat as f64 * w.pwat + + r.factors.pressure as f64 * w.pressure; + assert!((r.score as f64 - expected.round()).abs() <= 1.0); + } + + #[test] + fn composite_10g_vs_24g_differ() { + let c = canonical_conditions(); + let r10 = composite_score(&c, b10g()); + let r24 = composite_score(&c, b24g()); + assert_ne!(r10.score, r24.score); + } +} diff --git a/rust/prop_grid_rs/src/scores_file.rs b/rust/prop_grid_rs/src/scores_file.rs new file mode 100644 index 00000000..c5a42663 --- /dev/null +++ b/rust/prop_grid_rs/src/scores_file.rs @@ -0,0 +1,284 @@ +//! Score-grid file format (v1). 1:1 port of +//! `lib/microwaveprop/propagation/scores_file.ex`. +//! +//! Layout (all little-endian): +//! magic : 4 bytes "NTMS" (literal already on disk) +//! version : 1 byte 0x01 +//! band_mhz : 4 bytes uint32 +//! valid_time : 8 bytes int64 unix seconds +//! lat_min : 4 bytes float32 +//! lon_min : 4 bytes float32 +//! step_deg : 4 bytes float32 +//! n_rows : 2 bytes uint16 (lat) +//! n_cols : 2 bytes uint16 (lon) +//! scores : n_rows × n_cols bytes (uint8, 0-100 or 255 = no-data) +//! +//! Writes go through `rename(2)` — we write to `path.tmp.>` +//! then atomically rename to the final path. On the shared NFSv4 mount +//! this is safe for concurrent readers. + +use std::fs::OpenOptions; +use std::io::{BufWriter, Write}; +use std::path::{Path, PathBuf}; +use std::sync::atomic::{AtomicU64, Ordering}; +use std::time::{SystemTime, UNIX_EPOCH}; + +use chrono::{DateTime, Utc}; + +use crate::grid; + +pub const MAGIC: &[u8; 4] = b"NTMS"; +pub const VERSION: u8 = 1; +pub const NO_DATA: u8 = 255; + +static UNIQUE: AtomicU64 = AtomicU64::new(1); + +/// Absolute path `//.ntms`. +pub fn path_for(base_dir: &Path, band_mhz: u32, valid_time: DateTime) -> PathBuf { + let iso = valid_time + .format("%Y-%m-%dT%H:%M:%SZ") + .to_string(); + base_dir + .join(band_mhz.to_string()) + .join(format!("{iso}.ntms")) +} + +#[derive(Debug, thiserror::Error)] +pub enum WriteError { + #[error("io: {0}")] + Io(#[from] std::io::Error), + #[error("grid dimensions exceed uint16 (rows={rows}, cols={cols})")] + GridTooLarge { rows: usize, cols: usize }, +} + +/// Encode the full on-disk representation into a `Vec`. Separated +/// from the file I/O so tests can round-trip without a tempdir. +pub fn encode( + band_mhz: u32, + valid_time: DateTime, + scores: &[ScorePoint], +) -> Result, WriteError> { + let spec = grid::wgrib2_grid_spec(); + let n_rows: u16 = spec + .lat_count + .try_into() + .map_err(|_| WriteError::GridTooLarge { + rows: spec.lat_count, + cols: spec.lon_count, + })?; + let n_cols: u16 = spec + .lon_count + .try_into() + .map_err(|_| WriteError::GridTooLarge { + rows: spec.lat_count, + cols: spec.lon_count, + })?; + + let row_count = n_rows as usize; + let col_count = n_cols as usize; + let mut body = vec![NO_DATA; row_count * col_count]; + + for s in scores { + let row = (((s.lat - spec.lat_start) / spec.lat_step).round()) as isize; + let col = (((s.lon - spec.lon_start) / spec.lon_step).round()) as isize; + if row < 0 || col < 0 || row as usize >= row_count || col as usize >= col_count { + continue; + } + body[(row as usize) * col_count + (col as usize)] = clamp_score(s.score); + } + + let mut out = + Vec::with_capacity(4 + 1 + 4 + 8 + 4 + 4 + 4 + 2 + 2 + body.len()); + out.extend_from_slice(MAGIC); + out.push(VERSION); + out.extend_from_slice(&band_mhz.to_le_bytes()); + out.extend_from_slice(&valid_time.timestamp().to_le_bytes()); + out.extend_from_slice(&(spec.lat_start as f32).to_le_bytes()); + out.extend_from_slice(&(spec.lon_start as f32).to_le_bytes()); + out.extend_from_slice(&(spec.lon_step as f32).to_le_bytes()); + out.extend_from_slice(&n_rows.to_le_bytes()); + out.extend_from_slice(&n_cols.to_le_bytes()); + out.extend_from_slice(&body); + Ok(out) +} + +/// Atomic write to `//.ntms`. Creates the band +/// subdirectory on demand. +pub fn write_atomic( + base_dir: &Path, + band_mhz: u32, + valid_time: DateTime, + scores: &[ScorePoint], +) -> Result<(), WriteError> { + let bytes = encode(band_mhz, valid_time, scores)?; + let final_path = path_for(base_dir, band_mhz, valid_time); + let parent = final_path + .parent() + .ok_or_else(|| std::io::Error::other("path has no parent"))?; + std::fs::create_dir_all(parent)?; + + let nanos = SystemTime::now() + .duration_since(UNIX_EPOCH) + .map(|d| d.as_nanos()) + .unwrap_or(0); + let uniq = UNIQUE.fetch_add(1, Ordering::Relaxed); + let tmp_path = final_path.with_extension(format!("ntms.tmp.{nanos}.{uniq}")); + + { + let file = OpenOptions::new() + .write(true) + .create_new(true) + .open(&tmp_path)?; + let mut w = BufWriter::new(file); + w.write_all(&bytes)?; + w.flush()?; + } + + match std::fs::rename(&tmp_path, &final_path) { + Ok(()) => Ok(()), + Err(e) => { + let _ = std::fs::remove_file(&tmp_path); + Err(WriteError::Io(e)) + } + } +} + +#[derive(Debug, Clone, Copy)] +pub struct ScorePoint { + pub lat: f64, + pub lon: f64, + pub score: i32, +} + +fn clamp_score(s: i32) -> u8 { + if s < 0 { + 0 + } else if s > 100 { + 100 + } else { + s as u8 + } +} + +/// Minimal decoder used only in tests + the shadow comparator. Returns +/// the raw row-major body plus the metadata. Mirrors the Elixir +/// `decode/1` sans Access-behavior niceties. +pub fn decode(bytes: &[u8]) -> Option { + use byteorder::{ByteOrder, LittleEndian}; + if bytes.len() < 4 + 1 + 4 + 8 + 4 + 4 + 4 + 2 + 2 { + return None; + } + if &bytes[0..4] != MAGIC { + return None; + } + if bytes[4] != VERSION { + return None; + } + let band_mhz = LittleEndian::read_u32(&bytes[5..9]); + let valid_time_unix = LittleEndian::read_i64(&bytes[9..17]); + let lat_min = LittleEndian::read_f32(&bytes[17..21]); + let lon_min = LittleEndian::read_f32(&bytes[21..25]); + let step = LittleEndian::read_f32(&bytes[25..29]); + let n_rows = LittleEndian::read_u16(&bytes[29..31]); + let n_cols = LittleEndian::read_u16(&bytes[31..33]); + let body_len = (n_rows as usize) * (n_cols as usize); + if bytes.len() < 33 + body_len { + return None; + } + let valid_time = DateTime::::from_timestamp(valid_time_unix, 0)?; + Some(Decoded { + band_mhz, + valid_time, + lat_min, + lon_min, + step, + n_rows, + n_cols, + body: bytes[33..33 + body_len].to_vec(), + }) +} + +#[derive(Debug)] +pub struct Decoded { + pub band_mhz: u32, + pub valid_time: DateTime, + pub lat_min: f32, + pub lon_min: f32, + pub step: f32, + pub n_rows: u16, + pub n_cols: u16, + pub body: Vec, +} + +#[cfg(test)] +mod tests { + use super::*; + use chrono::TimeZone; + + #[test] + fn round_trip_header() { + let scores = vec![ScorePoint { + lat: 32.0, + lon: -97.0, + score: 73, + }]; + let ts = Utc.with_ymd_and_hms(2026, 4, 19, 15, 0, 0).unwrap(); + let bytes = encode(10_000, ts, &scores).unwrap(); + let decoded = decode(&bytes).unwrap(); + assert_eq!(decoded.band_mhz, 10_000); + assert_eq!(decoded.valid_time, ts); + assert_eq!(decoded.n_rows, 201); + assert_eq!(decoded.n_cols, 473); + // Score at (32, -97): row = (32-25)/0.125 = 56, col = (-97 - -125)/0.125 = 224 + let idx = 56 * 473 + 224; + assert_eq!(decoded.body[idx], 73); + } + + #[test] + fn missing_cells_use_no_data_sentinel() { + let scores: Vec = vec![]; + let ts = Utc.with_ymd_and_hms(2026, 4, 19, 15, 0, 0).unwrap(); + let bytes = encode(10_000, ts, &scores).unwrap(); + let decoded = decode(&bytes).unwrap(); + assert!(decoded.body.iter().all(|&b| b == NO_DATA)); + } + + #[test] + fn clamps_out_of_range_scores() { + let scores = vec![ + ScorePoint { lat: 25.0, lon: -125.0, score: -5 }, + ScorePoint { lat: 25.125, lon: -125.0, score: 150 }, + ]; + let ts = Utc.with_ymd_and_hms(2026, 4, 19, 15, 0, 0).unwrap(); + let bytes = encode(10_000, ts, &scores).unwrap(); + let decoded = decode(&bytes).unwrap(); + assert_eq!(decoded.body[0], 0); + // row = 1, col = 0 → idx = 1 * 473 + assert_eq!(decoded.body[473], 100); + } + + #[test] + fn path_for_layout() { + let ts = Utc.with_ymd_and_hms(2026, 4, 19, 15, 0, 0).unwrap(); + let p = path_for(Path::new("/data/scores"), 10_000, ts); + assert_eq!( + p.to_str().unwrap(), + "/data/scores/10000/2026-04-19T15:00:00Z.ntms" + ); + } + + #[test] + fn atomic_write_and_decode() { + let dir = tempfile::tempdir().unwrap(); + let ts = Utc.with_ymd_and_hms(2026, 4, 19, 15, 0, 0).unwrap(); + let scores = vec![ScorePoint { lat: 32.0, lon: -97.0, score: 73 }]; + + write_atomic(dir.path(), 10_000, ts, &scores).unwrap(); + + let path = path_for(dir.path(), 10_000, ts); + let bytes = std::fs::read(&path).unwrap(); + let decoded = decode(&bytes).unwrap(); + assert_eq!(decoded.band_mhz, 10_000); + assert_eq!(decoded.body[56 * 473 + 224], 73); + } +} diff --git a/rust/prop_grid_rs/src/sounding_params.rs b/rust/prop_grid_rs/src/sounding_params.rs new file mode 100644 index 00000000..cbc73231 --- /dev/null +++ b/rust/prop_grid_rs/src/sounding_params.rs @@ -0,0 +1,125 @@ +//! Minimum-refractivity-gradient derivation from pressure-level +//! profiles. 1:1 port of the subset of +//! `lib/microwaveprop/weather/sounding_params.ex` needed by the +//! propagation grid scorer: `min_refractivity_gradient` per cell. +//! +//! The full Elixir module derives CAPE/LI/K-index/ducts — those are +//! used for contact enrichment and f00's native-duct merge, neither of +//! which is in the Rust f01..f18 scope. + +#[derive(Debug, Clone)] +pub struct Level { + pub pres_mb: f64, + pub hght_m: f64, + pub tmpc: f64, + pub dwpc: Option, +} + +/// Buck equation for saturation vapor pressure (hPa) at temperature in °C. +pub fn sat_vap_pres(t_c: f64) -> f64 { + 6.1121 * ((18.678 - t_c / 234.5) * (t_c / (257.14 + t_c))).exp() +} + +/// Minimum refractivity gradient (dN/dh, N-units per km) from a pressure +/// profile. Returns `None` if fewer than 3 usable levels or adjacent +/// pairs too close in height. Levels are internally sorted descending by +/// pressure (surface first) to match the Elixir orientation. +pub fn min_refractivity_gradient(mut levels: Vec) -> Option { + levels.retain(|l| l.dwpc.is_some()); + if levels.len() < 3 { + return None; + } + levels.sort_by(|a, b| { + b.pres_mb + .partial_cmp(&a.pres_mb) + .unwrap_or(std::cmp::Ordering::Equal) + }); + let sfc_hght = levels[0].hght_m; + + let refract: Vec<(f64, f64)> = levels + .iter() + .map(|l| { + let t_k = l.tmpc + 273.15; + let e = sat_vap_pres(l.dwpc.expect("filtered above")); + let n = 77.6 * l.pres_mb / t_k + 3.73e5 * e / (t_k * t_k); + let h_agl = l.hght_m - sfc_hght; + (h_agl, n) + }) + .collect(); + + let mut min_grad: Option = None; + for pair in refract.windows(2) { + let (h0, n0) = pair[0]; + let (h1, n1) = pair[1]; + let dh_km = (h1 - h0) / 1000.0; + if dh_km.abs() < 0.01 { + continue; + } + let g = (n1 - n0) / dh_km; + min_grad = Some(match min_grad { + Some(prev) if prev <= g => prev, + _ => g, + }); + } + min_grad +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn buck_sat_vap_at_freezing() { + // At 0 °C, Buck gives ~6.112 hPa. + let e = sat_vap_pres(0.0); + assert!((e - 6.112).abs() < 0.02, "{e}"); + } + + #[test] + fn nil_for_too_few_levels() { + let levels = vec![Level { + pres_mb: 1000.0, + hght_m: 0.0, + tmpc: 25.0, + dwpc: Some(20.0), + }]; + assert!(min_refractivity_gradient(levels).is_none()); + } + + #[test] + fn typical_tropospheric_gradient_is_negative() { + // Synthetic profile: refractivity decreases with altitude → negative + // dN/dh. Realistic numbers from a standard atmosphere. + let levels = vec![ + Level { pres_mb: 1000.0, hght_m: 100.0, tmpc: 25.0, dwpc: Some(20.0) }, + Level { pres_mb: 925.0, hght_m: 800.0, tmpc: 20.0, dwpc: Some(14.0) }, + Level { pres_mb: 850.0, hght_m: 1500.0, tmpc: 15.0, dwpc: Some(8.0) }, + Level { pres_mb: 700.0, hght_m: 3000.0, tmpc: 5.0, dwpc: Some(-5.0) }, + ]; + let g = min_refractivity_gradient(levels).unwrap(); + assert!(g < 0.0, "gradient should be negative: {g}"); + assert!(g > -500.0, "no duct expected: {g}"); + } + + #[test] + fn surface_temperature_inversion_yields_stronger_gradient() { + // Surface inversion: cold moist near the ground, warm dry above — + // strong negative refractivity gradient. + let inversion = vec![ + Level { pres_mb: 1000.0, hght_m: 100.0, tmpc: 18.0, dwpc: Some(17.0) }, + Level { pres_mb: 990.0, hght_m: 300.0, tmpc: 22.0, dwpc: Some(5.0) }, + Level { pres_mb: 950.0, hght_m: 700.0, tmpc: 20.0, dwpc: Some(3.0) }, + Level { pres_mb: 850.0, hght_m: 1500.0, tmpc: 15.0, dwpc: Some(0.0) }, + ]; + let strong = min_refractivity_gradient(inversion).unwrap(); + + let neutral = vec![ + Level { pres_mb: 1000.0, hght_m: 100.0, tmpc: 25.0, dwpc: Some(18.0) }, + Level { pres_mb: 925.0, hght_m: 800.0, tmpc: 20.0, dwpc: Some(14.0) }, + Level { pres_mb: 850.0, hght_m: 1500.0, tmpc: 15.0, dwpc: Some(8.0) }, + ]; + let baseline = min_refractivity_gradient(neutral).unwrap(); + + assert!(strong < baseline, "inversion {strong} vs baseline {baseline}"); + } +} diff --git a/rust/prop_grid_rs/tests/scorer_golden.rs b/rust/prop_grid_rs/tests/scorer_golden.rs new file mode 100644 index 00000000..fb300b70 --- /dev/null +++ b/rust/prop_grid_rs/tests/scorer_golden.rs @@ -0,0 +1,84 @@ +//! Golden-fixture regression test for the Rust scorer. +//! +//! Reads `tests/scores.bincode` (produced by `mix rust.golden` from the +//! Elixir scorer) and asserts every (conditions, band) → score tuple +//! matches. Regenerate after any scorer or band_config change: +//! +//! mix rust.golden +//! cp priv/rust_golden/scores.bincode rust/prop_grid_rs/tests/ +//! +//! The format is documented in `lib/mix/tasks/rust.golden.ex`. + +use std::io::{Cursor, Read}; + +use byteorder::{LittleEndian, ReadBytesExt}; +use prop_grid_rs::{band_config, scorer}; + +fn opt(v: f64) -> Option { + if v.is_nan() { + None + } else { + Some(v) + } +} + +fn read_f64(r: &mut Cursor<&[u8]>) -> f64 { + r.read_f64::().expect("f64") +} + +fn read_u8(r: &mut Cursor<&[u8]>) -> u8 { + let mut buf = [0u8; 1]; + r.read_exact(&mut buf).expect("u8"); + buf[0] +} + +#[test] +fn matches_elixir_golden_fixture() { + let bytes = std::fs::read("tests/scores.bincode") + .expect("tests/scores.bincode — run `mix rust.golden` first"); + let mut r = Cursor::new(bytes.as_slice()); + + let n = r.read_u32::().expect("n_samples"); + assert!(n > 0, "empty fixture"); + + let mut mismatches: Vec = Vec::new(); + for i in 0..n { + let c = scorer::Conditions { + abs_humidity: read_f64(&mut r), + temp_f: read_f64(&mut r), + dewpoint_f: read_f64(&mut r), + wind_speed_kts: opt(read_f64(&mut r)), + sky_cover_pct: opt(read_f64(&mut r)), + utc_hour: read_u8(&mut r), + utc_minute: read_u8(&mut r), + month: read_u8(&mut r), + longitude: read_f64(&mut r), + latitude: opt(read_f64(&mut r)), + pressure_mb: opt(read_f64(&mut r)), + prev_pressure_mb: opt(read_f64(&mut r)), + rain_rate_mmhr: opt(read_f64(&mut r)), + min_refractivity_gradient: opt(read_f64(&mut r)), + bl_depth_m: opt(read_f64(&mut r)), + pwat_mm: opt(read_f64(&mut r)), + best_duct_band_ghz: opt(read_f64(&mut r)), + bulk_richardson: opt(read_f64(&mut r)), + }; + let band_mhz = r.read_u32::().expect("band_mhz"); + let expected = r.read_i32::().expect("expected_score"); + + let band = band_config::get(band_mhz).expect("band in config"); + let got = scorer::composite_score(&c, band).score; + if got != expected { + mismatches.push(format!( + "sample {i} band {band_mhz}: elixir={expected} rust={got}" + )); + } + } + + assert!( + mismatches.is_empty(), + "{n} samples checked, {} mismatches:\n{}", + mismatches.len(), + mismatches.join("\n") + ); +} diff --git a/rust/prop_grid_rs/tests/scores.bincode b/rust/prop_grid_rs/tests/scores.bincode new file mode 100644 index 0000000000000000000000000000000000000000..1a8be7f1bb9f6e95246f7bffdbb66d8b7c72cffa GIT binary patch literal 15069 zcmc(lKWh|07>D0WV!RO4L)5}Q@rpx11mht=A&H~$OiUnNE|;7dEG&e8jUa+nZlkSN zT3F}@5YpLhX_*folGbmK()rEII-XzP*j?z_xbzVzF{Xe`aOa$^nd(>Mfk-d z?A6R3vsL&ojw3XV;A6q2@ENgk46j0Q21jUI#u3V=X3OlP{Pf*yjYaq}RkQmn!eOgs z4_JiZ&zh~W2rs9{=;63`WsDvRyK`gKVLM{?+qgLv!Jbdh2>ZV$>^574o0leOgg2dk z5hl{Gm``~%Wjkfd9XD@!eshn?kd?$YhAy1fQ^rQWUoXR9^Ro+S)p}r5;_t>rzef-x z@k=a%Bz~Snki^3XlK2I-3X*uZ3X*uZ3X=Fmwj(6*jzy5f!_Ngt{1#gUNqmn*ki^3X zl6V+F5)UIt;$Z|yJd7ZTrxCowciBlPiGRoN8=^`NqjfGfiJ2z@y&E=zYN#Ueof~~DMuy#PPj-r1Yzj^7)Fr9!w8c2yqZEf zLJ|*GK@#tLXT%Muf+Qa92uXa_SELG(_l2 zcwhEGXcq~Zb;VMYt^SdTwQcB|CDoEntu^@@h%UYx(B=LDQg&;`c9s696 z#OIYYQUytTURfgul6bfyB=PwZ3Zx2>cpAY=`~o{8B=InUBpybP#KQ=Z_>7R=xZ8G| zr*1RN*JjLO&7HJ*7&&fGhAb9h(_WJIDnp~+6+2}(Jb3Ry+G@&CiC>A0evcqX;$Z|y ze21-qBpybP#KQ=Zco;zvzs!yWNqqhY8JU!lco;zv4?h Enum.map(& &1.forecast_hour) |> Enum.sort() == Enum.to_list(1..18) + assert Enum.all?(rows, &(&1.status == "queued")) + + # Raw schemaless selects return NaiveDateTime; compare against the + # equivalent naive form. + expected_naive = DateTime.to_naive(run_time) + + Enum.each(rows, fn r -> + expected = NaiveDateTime.add(expected_naive, r.forecast_hour * 3600, :second) + assert NaiveDateTime.compare(r.valid_time, expected) == :eq + end) + end + + test "is idempotent: reseeding the same run_time inserts nothing" do + run_time = ~U[2026-04-19 16:00:00Z] + assert {:ok, 18} = GridTaskEnqueuer.seed(run_time) + assert {:ok, 0} = GridTaskEnqueuer.seed(run_time) + + count = Repo.one(from(g in "grid_tasks", where: g.run_time == ^run_time, select: count())) + + assert count == 18 + end + + test "never emits a row for fh=0 — Elixir keeps that chain step" do + run_time = ~U[2026-04-19 17:00:00Z] + assert {:ok, 18} = GridTaskEnqueuer.seed(run_time) + + f0_count = Repo.one(from(g in "grid_tasks", where: g.run_time == ^run_time and g.forecast_hour == 0, select: count())) + + assert f0_count == 0 + end +end