diff --git a/config/config.exs b/config/config.exs index 522cea89..f00dbf44 100644 --- a/config/config.exs +++ b/config/config.exs @@ -45,7 +45,18 @@ config :microwaveprop, MicrowavepropWeb.Endpoint, config :microwaveprop, Oban, engine: Oban.Pro.Engines.Smart, repo: Microwaveprop.Repo, - queues: [solar: 1, weather: 20, enqueue: 1, hrrr: 5, terrain: 4, commercial: 2, iemre: 10, propagation: 1, admin: 1], + queues: [ + solar: 1, + weather: 20, + enqueue: 1, + hrrr: 5, + terrain: 4, + commercial: 2, + iemre: 10, + propagation: 1, + admin: 1, + nexrad: 2 + ], plugins: [ {Oban.Plugins.Pruner, max_age: 3600 * 24}, {Oban.Plugins.Lifeline, rescue_after: to_timeout(minute: 30)}, diff --git a/config/runtime.exs b/config/runtime.exs index 1fd043da..2d88779f 100644 --- a/config/runtime.exs +++ b/config/runtime.exs @@ -151,7 +151,8 @@ if config_env() == :prod do ], rtma: 2, backfill_enqueue: 1, - admin: 1 + admin: 1, + nexrad: 2 ], plugins: [ {Oban.Plugins.Pruner, max_age: 3600 * 24}, diff --git a/config/test.exs b/config/test.exs index cb49c0c2..4a409f04 100644 --- a/config/test.exs +++ b/config/test.exs @@ -36,6 +36,7 @@ config :microwaveprop, hrrr_req_options: [plug: {Req.Test, Microwaveprop.Weather config :microwaveprop, iem_req_options: [plug: {Req.Test, Microwaveprop.Weather.IemClient}, retry: false] # Route HTTP requests through Req.Test stubs +config :microwaveprop, nexrad_req_options: [plug: {Req.Test, Microwaveprop.Weather.NexradClient}, retry: false] config :microwaveprop, solar_req_options: [plug: {Req.Test, Microwaveprop.Weather.SolarClient}] config :microwaveprop, srtm_req_options: [plug: {Req.Test, Microwaveprop.Terrain.Srtm}, retry: false] config :microwaveprop, start_freshness_monitor: false diff --git a/lib/microwaveprop/backtest.ex b/lib/microwaveprop/backtest.ex index b5e57fd8..c6f7d94a 100644 --- a/lib/microwaveprop/backtest.ex +++ b/lib/microwaveprop/backtest.ex @@ -385,11 +385,9 @@ defmodule Microwaveprop.Backtest do defp eval_feature_for_contact(_feature, _contact), do: nil # pos1 is sometimes stored with "lon" and sometimes with "lng". - defp pos_to_latlon(%{"lat" => lat, "lon" => lon}) when is_number(lat) and is_number(lon), - do: {lat * 1.0, lon * 1.0} + defp pos_to_latlon(%{"lat" => lat, "lon" => lon}) when is_number(lat) and is_number(lon), do: {lat * 1.0, lon * 1.0} - defp pos_to_latlon(%{"lat" => lat, "lng" => lon}) when is_number(lat) and is_number(lon), - do: {lat * 1.0, lon * 1.0} + defp pos_to_latlon(%{"lat" => lat, "lng" => lon}) when is_number(lat) and is_number(lon), do: {lat * 1.0, lon * 1.0} defp pos_to_latlon(_), do: nil diff --git a/lib/microwaveprop/backtest/features.ex b/lib/microwaveprop/backtest/features.ex index 30748a5b..33a992c5 100644 --- a/lib/microwaveprop/backtest/features.ex +++ b/lib/microwaveprop/backtest/features.ex @@ -24,6 +24,7 @@ defmodule Microwaveprop.Backtest.Features do alias Microwaveprop.Weather alias Microwaveprop.Weather.HrrrClimatology alias Microwaveprop.Weather.HrrrNativeProfile + alias Microwaveprop.Weather.NexradObservation @doc """ Returns a map of all backtestable features: `%{name => fun/3}`. @@ -34,7 +35,8 @@ defmodule Microwaveprop.Backtest.Features do def all_features do excluded = MapSet.new([:duct_usable_for_band, :distance_to_front, :parallel_to_front]) - __MODULE__.__info__(:functions) + :functions + |> __MODULE__.__info__() |> Enum.filter(fn {_name, arity} -> arity == 3 end) |> Enum.reject(fn {name, _} -> MapSet.member?(excluded, name) end) |> Map.new(fn {name, _} -> @@ -51,10 +53,8 @@ defmodule Microwaveprop.Backtest.Features do """ @spec naive_gradient(float, float, DateTime.t()) :: float | nil def naive_gradient(lat, lon, valid_time) do - with %{min_refractivity_gradient: grad} when is_float(grad) <- - Weather.find_nearest_hrrr(lat, lon, valid_time) do - grad - else + case Weather.find_nearest_hrrr(lat, lon, valid_time) do + %{min_refractivity_gradient: grad} when is_float(grad) -> grad _ -> nil end end @@ -67,10 +67,8 @@ defmodule Microwaveprop.Backtest.Features do """ @spec td_depression(float, float, DateTime.t()) :: float | nil def td_depression(lat, lon, valid_time) do - with %{surface_temp_c: t, surface_dewpoint_c: td} when is_float(t) and is_float(td) <- - Weather.find_nearest_hrrr(lat, lon, valid_time) do - t - td - else + case Weather.find_nearest_hrrr(lat, lon, valid_time) do + %{surface_temp_c: t, surface_dewpoint_c: td} when is_float(t) and is_float(td) -> t - td _ -> nil end end @@ -96,10 +94,8 @@ defmodule Microwaveprop.Backtest.Features do """ @spec pressure(float, float, DateTime.t()) :: float | nil def pressure(lat, lon, valid_time) do - with %{surface_pressure_mb: p} when is_float(p) <- - Weather.find_nearest_hrrr(lat, lon, valid_time) do - p - else + case Weather.find_nearest_hrrr(lat, lon, valid_time) do + %{surface_pressure_mb: p} when is_float(p) -> p _ -> nil end end @@ -141,10 +137,8 @@ defmodule Microwaveprop.Backtest.Features do """ @spec bulk_richardson(float, float, DateTime.t()) :: float | nil def bulk_richardson(lat, lon, valid_time) do - with %HrrrNativeProfile{bulk_richardson: ri} when is_float(ri) <- - find_nearest_native(lat, lon, valid_time) do - ri - else + case find_nearest_native(lat, lon, valid_time) do + %HrrrNativeProfile{bulk_richardson: ri} when is_float(ri) -> ri _ -> nil end end @@ -156,10 +150,8 @@ defmodule Microwaveprop.Backtest.Features do """ @spec theta_e_jump(float, float, DateTime.t()) :: float | nil def theta_e_jump(lat, lon, valid_time) do - with %HrrrNativeProfile{theta_e_jump_k: jump} when is_float(jump) <- - find_nearest_native(lat, lon, valid_time) do - jump - else + case find_nearest_native(lat, lon, valid_time) do + %HrrrNativeProfile{theta_e_jump_k: jump} when is_float(jump) -> jump _ -> nil end end @@ -169,10 +161,8 @@ defmodule Microwaveprop.Backtest.Features do """ @spec shear_at_top(float, float, DateTime.t()) :: float | nil def shear_at_top(lat, lon, valid_time) do - with %HrrrNativeProfile{shear_at_top_ms: shear} when is_float(shear) <- - find_nearest_native(lat, lon, valid_time) do - shear - else + case find_nearest_native(lat, lon, valid_time) do + %HrrrNativeProfile{shear_at_top_ms: shear} when is_float(shear) -> shear _ -> nil end end @@ -185,11 +175,12 @@ defmodule Microwaveprop.Backtest.Features do """ @spec duct_thickness(float, float, DateTime.t()) :: float | nil def duct_thickness(lat, lon, valid_time) do - with %HrrrNativeProfile{ducts: ducts} when is_list(ducts) and ducts != [] <- - find_nearest_native(lat, lon, valid_time) do - ducts |> Enum.map(& &1["thickness_m"]) |> Enum.max(fn -> nil end) - else - _ -> nil + case find_nearest_native(lat, lon, valid_time) do + %HrrrNativeProfile{ducts: ducts} when is_list(ducts) and ducts != [] -> + ducts |> Enum.map(& &1["thickness_m"]) |> Enum.max(fn -> nil end) + + _ -> + nil end end @@ -200,10 +191,8 @@ defmodule Microwaveprop.Backtest.Features do """ @spec best_duct_freq(float, float, DateTime.t()) :: float | nil def best_duct_freq(lat, lon, valid_time) do - with %HrrrNativeProfile{best_duct_band_ghz: f} when is_float(f) <- - find_nearest_native(lat, lon, valid_time) do - f - else + case find_nearest_native(lat, lon, valid_time) do + %HrrrNativeProfile{best_duct_band_ghz: f} when is_float(f) -> f _ -> nil end end @@ -216,10 +205,8 @@ defmodule Microwaveprop.Backtest.Features do """ @spec duct_usable_for_band(float, float, DateTime.t(), float) :: float | nil def duct_usable_for_band(lat, lon, valid_time, band_ghz) do - with %HrrrNativeProfile{best_duct_band_ghz: f} when is_float(f) <- - find_nearest_native(lat, lon, valid_time) do - if f <= band_ghz, do: 1.0, else: 0.0 - else + case find_nearest_native(lat, lon, valid_time) do + %HrrrNativeProfile{best_duct_band_ghz: f} when is_float(f) -> if f <= band_ghz, do: 1.0, else: 0.0 _ -> nil end end @@ -292,6 +279,45 @@ defmodule Microwaveprop.Backtest.Features do @doc "Duct usable for 47 GHz." def duct_usable_47ghz(lat, lon, vt), do: duct_usable_for_band(lat, lon, vt, 47.0) + @doc """ + Spatial texture variance from nearest NEXRAD observation. + + Higher variance indicates more boundary-layer turbulence, which is + generally worse for microwave propagation (disrupts stable ducting + layers). Returns nil if no observation found within +/- 15 minutes + and +/- 0.1 degrees. + """ + @spec nexrad_texture(float, float, DateTime.t()) :: float | nil + def nexrad_texture(lat, lon, valid_time) do + dlat = 0.1 + dlon = 0.1 + time_start = DateTime.add(valid_time, -900, :second) + time_end = DateTime.add(valid_time, 900, :second) + + NexradObservation + |> where( + [n], + n.lat >= ^(lat - dlat) and n.lat <= ^(lat + dlat) and + n.lon >= ^(lon - dlon) and n.lon <= ^(lon + dlon) and + n.observed_at >= ^time_start and n.observed_at <= ^time_end + ) + |> order_by([n], + asc: + fragment( + "ABS(? - ?) + ABS(? - ?) + ABS(EXTRACT(EPOCH FROM ? - ?))", + n.lat, + ^lat, + n.lon, + ^lon, + n.observed_at, + ^valid_time + ) + ) + |> limit(1) + |> select([n], n.texture_variance) + |> Repo.one() + end + defp find_nearest_native(lat, lon, valid_time) do dlat = 0.07 dlon = 0.07 diff --git a/lib/microwaveprop/propagation/duct.ex b/lib/microwaveprop/propagation/duct.ex index 922bff18..8a9b984a 100644 --- a/lib/microwaveprop/propagation/duct.ex +++ b/lib/microwaveprop/propagation/duct.ex @@ -114,8 +114,7 @@ defmodule Microwaveprop.Propagation.Duct do where d is duct thickness in meters, ΔM is the M-deficit. Returns a very high frequency (999 GHz) for degenerate ducts. """ - def min_trapped_frequency_ghz(%{thickness_m: d, m_deficit: delta_m}) - when d > 0 and delta_m > 0 do + def min_trapped_frequency_ghz(%{thickness_m: d, m_deficit: delta_m}) when d > 0 and delta_m > 0 do # Maximum trapped wavelength (meters) lambda_max = 2.5 * d * :math.sqrt(delta_m * 1.0e-6) diff --git a/lib/microwaveprop/propagation/inversion.ex b/lib/microwaveprop/propagation/inversion.ex index fd9ff782..20d49262 100644 --- a/lib/microwaveprop/propagation/inversion.ex +++ b/lib/microwaveprop/propagation/inversion.ex @@ -33,13 +33,15 @@ defmodule Microwaveprop.Propagation.Inversion do def find_inversion_top(%{level_count: n}) when n < 2, do: :none def find_inversion_top(%{heights_m: heights, temp_k: temps}) do - levels = Enum.zip(heights, temps) |> Enum.with_index() + levels = heights |> Enum.zip(temps) |> Enum.with_index() # Walk upward looking for the first level where T stops increasing. # "Inversion base" is the level where dT/dz first turns positive. # "Inversion top" is the level where dT/dz turns negative again. case find_inversion_region(levels) do - nil -> :none + nil -> + :none + {base_idx, top_idx} -> {:ok, %{ diff --git a/lib/microwaveprop/propagation/recalibrator.ex b/lib/microwaveprop/propagation/recalibrator.ex new file mode 100644 index 00000000..c493caa4 --- /dev/null +++ b/lib/microwaveprop/propagation/recalibrator.ex @@ -0,0 +1,315 @@ +defmodule Microwaveprop.Propagation.Recalibrator do + @moduledoc """ + Fits new scoring weights via gradient descent on the historical QSO corpus. + + Uses a logistic regression model: `p = sigmoid(w · factors)` where + `factors` is a 10-element vector of individual factor scores (0-100) + from the existing `Scorer` functions, and `w` is the weight vector + to learn. + + QSOs are labeled as positive examples (1.0). Matched random baseline + samples (same locations, ±30 days via `Backtest.random_baseline/2`) + serve as negative examples (0.0). + + The loss is binary cross-entropy. Training uses simple gradient descent + on CPU — fine for 10 parameters. + + Output: a normalized weights map `%{humidity: 0.XX, ...}` summing to 1.0, + along with train and validation losses. + """ + + import Ecto.Query + + alias Microwaveprop.Backtest + alias Microwaveprop.Propagation.BandConfig + alias Microwaveprop.Propagation.Scorer + alias Microwaveprop.Radio.Contact + alias Microwaveprop.Repo + alias Microwaveprop.Weather + + require Logger + + @factor_keys ~w(humidity time_of_day td_depression refractivity sky season wind rain pwat pressure)a + + @doc """ + Fits new weights using gradient descent on the QSO corpus. + + ## Options + + * `:sample_size` - max contacts to load (default: 5000) + * `:learning_rate` - gradient descent step size (default: 0.01) + * `:epochs` - number of training iterations (default: 2000) + + ## Returns + + A map with `:weights`, `:train_loss`, `:val_loss`, and `:initial_loss`. + """ + @spec fit(keyword()) :: %{weights: map(), train_loss: float(), val_loss: float(), initial_loss: float()} + def fit(opts \\ []) do + sample_size = Keyword.get(opts, :sample_size, 5000) + learning_rate = Keyword.get(opts, :learning_rate, 0.01) + epochs = Keyword.get(opts, :epochs, 2000) + + Logger.info("Recalibrator: loading contacts and computing factor vectors...") + + # Load contacts with HRRR profiles + {positives, contacts} = load_positive_samples(sample_size) + Logger.info("Recalibrator: #{length(positives)} positive samples from #{length(contacts)} contacts") + + # Generate matched negatives + negatives = generate_negative_samples(contacts, length(positives)) + Logger.info("Recalibrator: #{length(negatives)} negative samples") + + if positives == [] or negatives == [] do + Logger.warning("Recalibrator: insufficient data, returning current weights") + weights = BandConfig.weights() + + %{ + weights: weights, + train_loss: 0.0, + val_loss: 0.0, + initial_loss: 0.0 + } + else + run_training(positives, negatives, learning_rate, epochs) + end + end + + @doc """ + Computes the 10-element factor score vector from an HRRR profile and timestamp. + + Uses `Scorer` functions with a default 10 GHz band config (beneficial humidity). + Returns a list of 10 floats, each in [0, 100]. + """ + @spec compute_factors(map(), DateTime.t()) :: [number()] + def compute_factors(profile, timestamp) do + band_config = BandConfig.get(10_000) + + temp_f = Scorer.c_to_f(profile.surface_temp_c) + dewpoint_f = Scorer.c_to_f(profile.surface_dewpoint_c) + + abs_humidity = + if profile.surface_temp_c && profile.surface_dewpoint_c do + Scorer.absolute_humidity(profile.surface_temp_c, profile.surface_dewpoint_c) + else + 10.0 + end + + {tod_score, _label} = + Scorer.score_time_of_day( + timestamp.hour, + timestamp.minute, + timestamp.month, + profile.lon || -97.0 + ) + + [ + Scorer.score_humidity(abs_humidity, band_config), + tod_score, + Scorer.score_td_depression(temp_f || 70, dewpoint_f || 60, band_config), + Scorer.score_refractivity(profile.min_refractivity_gradient, profile.hpbl_m, band_config), + Scorer.score_sky(nil), + Scorer.score_season(timestamp.month, profile.lat, profile.lon, band_config), + Scorer.score_wind(nil), + Scorer.score_rain(0.0, band_config), + Scorer.score_pwat(profile.pwat_mm, band_config), + Scorer.score_pressure(profile.surface_pressure_mb, nil) + ] + end + + @doc """ + Runs gradient descent on pre-built factor vectors. + + `positives` and `negatives` are lists of 10-element factor score lists + (each score in [0, 100]). Returns the same result shape as `fit/1`. + + This is the core training function, separated from data loading so it + can be tested independently. + """ + @spec train([list()], [list()], keyword()) :: %{ + weights: map(), + train_loss: float(), + val_loss: float(), + initial_loss: float() + } + def train(positives, negatives, opts \\ []) do + learning_rate = Keyword.get(opts, :learning_rate, 0.01) + epochs = Keyword.get(opts, :epochs, 2000) + run_training(positives, negatives, learning_rate, epochs) + end + + # ── Private ────────────────────────────────────────────────────── + + defp load_positive_samples(sample_size) do + contacts = + Contact + |> where([c], not is_nil(c.pos1)) + |> where([c], c.distance_km < 3000) + |> order_by([c], desc: c.qso_timestamp) + |> limit(^sample_size) + |> Repo.all() + + factor_vectors = + contacts + |> Enum.map(&contact_to_factors/1) + |> Enum.reject(&is_nil/1) + + {factor_vectors, contacts} + end + + defp contact_to_factors(%Contact{pos1: pos, qso_timestamp: ts}) do + lat = pos["lat"] + lon = pos["lon"] || pos["lng"] + + if lat && lon do + case Weather.find_nearest_hrrr(lat, lon, ts) do + nil -> nil + profile -> compute_factors(profile, ts) + end + end + end + + defp generate_negative_samples(contacts, n) do + baselines = Backtest.random_baseline(n, sample_size: length(contacts)) + + baselines + |> Enum.map(fn {lat, lon, time} -> + case Weather.find_nearest_hrrr(lat, lon, time) do + nil -> nil + profile -> compute_factors(profile, time) + end + end) + |> Enum.reject(&is_nil/1) + end + + defp run_training(positives, negatives, learning_rate, epochs) do + # Build feature matrix and label vector + all_features = positives ++ negatives + labels = List.duplicate(1.0, length(positives)) ++ List.duplicate(0.0, length(negatives)) + + n = length(all_features) + + # Normalize factor scores to [0, 1] + x = all_features |> Nx.tensor(type: :f32) |> Nx.divide(100.0) + y = labels |> Nx.tensor(type: :f32) |> Nx.reshape({n, 1}) + + # Cross-validation: split by month index + # Use indices to split - train on ~80%, validate on ~20% + {train_x, train_y, val_x, val_y} = split_train_val(x, y, all_features, positives, negatives) + + # Initialize weights uniformly (10 factors, each 0.1) + w = 0.1 |> List.duplicate(10) |> Nx.tensor(type: :f32) |> Nx.reshape({10, 1}) + bias = Nx.tensor([0.0], type: :f32) + + # Compute initial loss + initial_loss = compute_loss(train_x, train_y, w, bias) + + # Training loop + {final_w, final_bias, train_loss} = + Enum.reduce(1..epochs, {w, bias, initial_loss}, fn epoch, {w_acc, b_acc, _prev_loss} -> + # Forward pass + logits = Nx.add(Nx.dot(train_x, w_acc), b_acc) + preds = Nx.sigmoid(logits) + + # Gradient of BCE loss + error = Nx.subtract(preds, train_y) + n_train = Nx.axis_size(train_x, 0) + + grad_w = train_x |> Nx.transpose() |> Nx.dot(error) |> Nx.divide(n_train) + grad_b = Nx.mean(error) + + # Update + new_w = Nx.subtract(w_acc, Nx.multiply(learning_rate, grad_w)) + new_b = Nx.subtract(b_acc, Nx.multiply(learning_rate, grad_b)) + + loss = compute_loss(train_x, train_y, new_w, new_b) + + if rem(epoch, 500) == 0 do + Logger.info("Recalibrator: epoch #{epoch}, train_loss=#{loss |> Nx.to_number() |> Float.round(6)}") + end + + {new_w, new_b, loss} + end) + + # Validation loss + val_loss = compute_loss(val_x, val_y, final_w, final_bias) + + # Extract and normalize weights + raw_weights = final_w |> Nx.abs() |> Nx.squeeze() + weight_sum = Nx.sum(raw_weights) + normalized = Nx.divide(raw_weights, weight_sum) + + weight_values = Nx.to_flat_list(normalized) + + weights_map = + @factor_keys + |> Enum.zip(weight_values) + |> Map.new(fn {k, v} -> {k, Float.round(v, 4)} end) + + %{ + weights: weights_map, + train_loss: train_loss |> Nx.to_number() |> to_float(), + val_loss: val_loss |> Nx.to_number() |> to_float(), + initial_loss: initial_loss |> Nx.to_number() |> to_float() + } + end + + defp compute_loss(x, y, w, bias) do + logits = Nx.add(Nx.dot(x, w), bias) + preds = Nx.sigmoid(logits) + + # Binary cross-entropy: -mean(y * log(p) + (1-y) * log(1-p)) + eps = 1.0e-7 + preds_clipped = Nx.clip(preds, eps, 1.0 - eps) + + Nx.negate( + Nx.mean( + Nx.add( + Nx.multiply(y, Nx.log(preds_clipped)), + Nx.multiply(Nx.subtract(1.0, y), Nx.log(Nx.subtract(1.0, preds_clipped))) + ) + ) + ) + end + + defp split_train_val(x, y, _all_features, positives, negatives) do + # Simple 80/20 split + n_pos = length(positives) + n_neg = length(negatives) + + train_pos = max(1, round(n_pos * 0.8)) + train_neg = max(1, round(n_neg * 0.8)) + + # Positives come first in the combined tensor, then negatives. + # Take 80% of positives and 80% of negatives for training. + train_pos_x = Nx.slice(x, [0, 0], [train_pos, 10]) + train_neg_x = Nx.slice(x, [n_pos, 0], [train_neg, 10]) + train_x = Nx.concatenate([train_pos_x, train_neg_x], axis: 0) + + train_pos_y = Nx.slice(y, [0, 0], [train_pos, 1]) + train_neg_y = Nx.slice(y, [n_pos, 0], [train_neg, 1]) + train_y = Nx.concatenate([train_pos_y, train_neg_y], axis: 0) + + val_pos_count = n_pos - train_pos + val_neg_count = n_neg - train_neg + + if val_pos_count <= 0 or val_neg_count <= 0 do + # Not enough data for validation — use training data for both + {train_x, train_y, train_x, train_y} + else + val_pos_x = Nx.slice(x, [train_pos, 0], [val_pos_count, 10]) + val_neg_x = Nx.slice(x, [n_pos + train_neg, 0], [val_neg_count, 10]) + val_x = Nx.concatenate([val_pos_x, val_neg_x], axis: 0) + + val_pos_y = Nx.slice(y, [train_pos, 0], [val_pos_count, 1]) + val_neg_y = Nx.slice(y, [n_pos + train_neg, 0], [val_neg_count, 1]) + val_y = Nx.concatenate([val_pos_y, val_neg_y], axis: 0) + + {train_x, train_y, val_x, val_y} + end + end + + defp to_float(n) when is_float(n), do: n + defp to_float(n) when is_integer(n), do: n * 1.0 + defp to_float(n), do: n * 1.0 +end diff --git a/lib/microwaveprop/propagation/scorer.ex b/lib/microwaveprop/propagation/scorer.ex index 501f908f..81e7d008 100644 --- a/lib/microwaveprop/propagation/scorer.ex +++ b/lib/microwaveprop/propagation/scorer.ex @@ -10,6 +10,7 @@ defmodule Microwaveprop.Propagation.Scorer do alias Microwaveprop.Propagation.BandConfig # ── Temperature conversion helpers ──────────────────────────────── + alias Microwaveprop.Propagation.Region @doc "Converts Fahrenheit to Celsius. Returns nil for nil input." @spec f_to_c(number() | nil) :: float() | nil @@ -220,8 +221,8 @@ defmodule Microwaveprop.Propagation.Scorer do region_mult = if is_number(lat) and is_number(lon) do - region = Microwaveprop.Propagation.Region.for_point(lat, lon) - Microwaveprop.Propagation.Region.seasonal_adjustment(region, month) + region = Region.for_point(lat, lon) + Region.seasonal_adjustment(region, month) else 1.0 end diff --git a/lib/microwaveprop/propagation/scorer_diff.ex b/lib/microwaveprop/propagation/scorer_diff.ex new file mode 100644 index 00000000..66af7be7 --- /dev/null +++ b/lib/microwaveprop/propagation/scorer_diff.ex @@ -0,0 +1,206 @@ +defmodule Microwaveprop.Propagation.ScorerDiff do + @moduledoc """ + Compare two weight sets on the same propagation scoring data. + + Loads the most recent HRRR frame's worth of grid points from + `propagation_scores` (which stores the per-factor scores in JSONB), + re-scores each with both weight sets, and returns diff statistics. + """ + + import Ecto.Query + + alias Microwaveprop.Propagation.GridScore + alias Microwaveprop.Repo + + @regression_threshold 5 + @worst_count 10 + + @doc """ + Compare old and new weights on recent propagation scores. + + Options: + * `:limit` — max number of scores to process (default: all for the latest valid_time) + + Returns a map with: + * `:summary` — aggregate stats (mean_diff, max_diff, regressions, improvements, total) + * `:band_diffs` — per-band breakdown keyed by band_mhz + * `:worst_regressions` — up to 10 rows with largest score decrease + """ + @spec compare(map(), map(), keyword()) :: %{ + summary: map(), + band_diffs: map(), + worst_regressions: [map()] + } + def compare(old_weights, new_weights, opts \\ []) do + scores = load_scores(opts) + diffs = compute_diffs(scores, old_weights, new_weights) + + %{ + summary: summarize(diffs), + band_diffs: band_breakdown(diffs), + worst_regressions: worst_regressions(diffs) + } + end + + @doc """ + Compute a weighted sum from a factors map and a weights map. + + Both maps use string keys. Returns a rounded integer. + """ + @spec weighted_sum(map(), map()) :: integer() + def weighted_sum(factors, weights) do + weights + |> Enum.reduce(0.0, fn {factor, weight}, acc -> + score = Map.get(factors, factor, 0) + acc + score * weight + end) + |> round() + end + + # ── Private ──────────────────────────────────────────────────────── + + defp load_scores(opts) do + limit = Keyword.get(opts, :limit) + + latest_time = + GridScore + |> select([g], max(g.valid_time)) + |> Repo.one() + + if is_nil(latest_time) do + [] + else + query = + GridScore + |> where([g], g.valid_time == ^latest_time) + |> order_by([g], [g.band_mhz, g.lat, g.lon]) + + query = if limit, do: limit(query, ^limit), else: query + + Repo.all(query) + end + end + + defp compute_diffs(scores, old_weights, new_weights) do + Enum.map(scores, fn score -> + factors = score.factors + old_score = weighted_sum(factors, old_weights) + new_score = weighted_sum(factors, new_weights) + diff = new_score - old_score + + %{ + lat: score.lat, + lon: score.lon, + band_mhz: score.band_mhz, + old_score: old_score, + new_score: new_score, + diff: diff + } + end) + end + + defp summarize([]), do: %{mean_diff: 0.0, max_diff: 0, regressions: 0, improvements: 0, total: 0} + + defp summarize(diffs) do + total = length(diffs) + abs_diffs = Enum.map(diffs, fn d -> abs(d.diff) end) + mean_diff = Float.round(Enum.sum(abs_diffs) / total, 2) + max_diff = Enum.max(abs_diffs) + regressions = Enum.count(diffs, fn d -> d.diff < -@regression_threshold end) + improvements = Enum.count(diffs, fn d -> d.diff > @regression_threshold end) + + %{ + mean_diff: mean_diff, + max_diff: max_diff, + regressions: regressions, + improvements: improvements, + total: total + } + end + + defp band_breakdown(diffs) do + diffs + |> Enum.group_by(& &1.band_mhz) + |> Map.new(fn {band_mhz, band_diffs} -> + abs_diffs = Enum.map(band_diffs, fn d -> abs(d.diff) end) + mean = Float.round(Enum.sum(abs_diffs) / length(abs_diffs), 2) + + regressions = Enum.filter(band_diffs, fn d -> d.diff < 0 end) + + max_regression = + case regressions do + [] -> 0 + list -> list |> Enum.map(fn d -> abs(d.diff) end) |> Enum.max() + end + + {band_mhz, %{mean_diff: mean, max_regression: max_regression}} + end) + end + + defp worst_regressions(diffs) do + diffs + |> Enum.filter(fn d -> d.diff < 0 end) + |> Enum.sort_by(fn d -> d.diff end) + |> Enum.take(@worst_count) + |> Enum.map(fn d -> + %{ + lat: d.lat, + lon: d.lon, + band_mhz: d.band_mhz, + old_score: d.old_score, + new_score: d.new_score + } + end) + end + + @doc """ + Format comparison results as a markdown report string. + """ + @spec to_markdown(%{summary: map(), band_diffs: map(), worst_regressions: [map()]}) :: String.t() + def to_markdown(result) do + s = result.summary + + summary = """ + ## Scorer Weight Comparison + + | Metric | Value | + |--------|-------| + | Total grid points | #{s.total} | + | Mean absolute diff | #{s.mean_diff} | + | Max absolute diff | #{s.max_diff} | + | Regressions (>#{@regression_threshold} pts) | #{s.regressions} | + | Improvements (>#{@regression_threshold} pts) | #{s.improvements} | + """ + + band_header = """ + + ## Per-Band Breakdown + + | Band (MHz) | Mean Diff | Max Regression | + |-----------|-----------|----------------| + """ + + band_rows = + result.band_diffs + |> Enum.sort_by(fn {band, _} -> band end) + |> Enum.map_join("\n", fn {band, stats} -> + "| #{band} | #{stats.mean_diff} | #{stats.max_regression} |" + end) + + worst_header = """ + + ## Worst Regressions + + | Lat | Lon | Band | Old | New | Diff | + |-----|-----|------|-----|-----|------| + """ + + worst_rows = + Enum.map_join(result.worst_regressions, "\n", fn r -> + diff = r.new_score - r.old_score + "| #{r.lat} | #{r.lon} | #{r.band_mhz} | #{r.old_score} | #{r.new_score} | #{diff} |" + end) + + summary <> band_header <> band_rows <> worst_header <> worst_rows + end +end diff --git a/lib/microwaveprop/release.ex b/lib/microwaveprop/release.ex index b5f71889..2aef471c 100644 --- a/lib/microwaveprop/release.ex +++ b/lib/microwaveprop/release.ex @@ -13,11 +13,13 @@ defmodule Microwaveprop.Release do bin/microwaveprop eval "Microwaveprop.Release.climatology" bin/microwaveprop eval "Microwaveprop.Release.native_backfill(500)" bin/microwaveprop eval "Microwaveprop.Release.native_derive" + bin/microwaveprop eval "Microwaveprop.Release.recalibrate" + bin/microwaveprop eval "Microwaveprop.Release.scorer_diff(\"{\\\"humidity\\\":0.18,...}\")" """ - @app :microwaveprop - alias Microwaveprop.Workers.AdminTaskWorker + @app :microwaveprop + def migrate do load_app() @@ -76,9 +78,7 @@ defmodule Microwaveprop.Release do for [year, month, day, hour, cnt] <- rows do args = %{year: year, month: month, day: day, hour: hour} - case Oban.insert( - Microwaveprop.Workers.HrrrNativeGridWorker.new(args, unique: [period: :infinity]) - ) do + case Oban.insert(Microwaveprop.Workers.HrrrNativeGridWorker.new(args, unique: [period: :infinity])) do {:ok, _} -> IO.puts(" #{year}-#{month}-#{day} #{hour}Z (#{cnt} contacts)") {:error, _} -> IO.puts(" #{year}-#{month}-#{day} #{hour}Z (skipped)") end @@ -87,6 +87,18 @@ defmodule Microwaveprop.Release do IO.puts("Done.") end + def recalibrate do + start_app() + {:ok, job} = Oban.insert(AdminTaskWorker.new(%{task: "recalibrate"})) + IO.puts("Enqueued recalibration (job #{job.id})") + end + + def scorer_diff(new_weights_json) when is_binary(new_weights_json) do + start_app() + {:ok, job} = Oban.insert(AdminTaskWorker.new(%{task: "scorer_diff", new_weights: new_weights_json})) + IO.puts("Enqueued scorer diff (job #{job.id})") + end + def native_derive(limit \\ 10_000) do start_app() {:ok, job} = Oban.insert(AdminTaskWorker.new(%{task: "native_derive", limit: limit})) diff --git a/lib/microwaveprop/weather/frontal_analysis.ex b/lib/microwaveprop/weather/frontal_analysis.ex index a9bcdde7..ef9307de 100644 --- a/lib/microwaveprop/weather/frontal_analysis.ex +++ b/lib/microwaveprop/weather/frontal_analysis.ex @@ -150,7 +150,8 @@ defmodule Microwaveprop.Weather.FrontalAnalysis do mask_binary = Nx.to_flat_list(mask) bearing_list = Nx.to_flat_list(bearing_deg) - Enum.zip(mask_binary, bearing_list) + mask_binary + |> Enum.zip(bearing_list) |> Enum.with_index() |> Enum.flat_map(fn {{is_front, bearing}, idx} -> if is_front == 1 do diff --git a/lib/microwaveprop/weather/hrrr_native_client.ex b/lib/microwaveprop/weather/hrrr_native_client.ex index faacb0c0..5d2c7fda 100644 --- a/lib/microwaveprop/weather/hrrr_native_client.ex +++ b/lib/microwaveprop/weather/hrrr_native_client.ex @@ -77,7 +77,8 @@ defmodule Microwaveprop.Weather.HrrrNativeClient do """ def build_native_profile(parsed) when is_map(parsed) do levels = - Enum.map(@native_levels, fn level -> + @native_levels + |> Enum.map(fn level -> level_str = "#{level} hybrid level" %{ @@ -105,9 +106,9 @@ defmodule Microwaveprop.Weather.HrrrNativeClient do u_wind_ms: Enum.map(levels, & &1.ugrd), v_wind_ms: Enum.map(levels, & &1.vgrd), tke_m2s2: Enum.map(levels, & &1.tke), - surface_temp_k: parsed["TMP:surface"] || List.first(levels) |> safe_get(:tmp), - surface_spfh: parsed["SPFH:2 m above ground"] || List.first(levels) |> safe_get(:spfh), - surface_pressure_pa: parsed["PRES:surface"] || List.first(levels) |> safe_get(:pres) + surface_temp_k: parsed["TMP:surface"] || levels |> List.first() |> safe_get(:tmp), + surface_spfh: parsed["SPFH:2 m above ground"] || levels |> List.first() |> safe_get(:spfh), + surface_pressure_pa: parsed["PRES:surface"] || levels |> List.first() |> safe_get(:pres) } end @@ -203,9 +204,12 @@ defmodule Microwaveprop.Weather.HrrrNativeClient do defp nearest_grid_cell(grid_data, lat, lon) do grid_data - |> Enum.min_by(fn {{glat, glon}, _} -> - :math.pow(glat - lat, 2) + :math.pow(glon - lon, 2) - end, fn -> nil end) + |> Enum.min_by( + fn {{glat, glon}, _} -> + :math.pow(glat - lat, 2) + :math.pow(glon - lon, 2) + end, + fn -> nil end + ) |> case do nil -> nil {_point, parsed} -> parsed diff --git a/lib/microwaveprop/weather/hrrr_native_profile.ex b/lib/microwaveprop/weather/hrrr_native_profile.ex index 08409396..1bd0b20d 100644 --- a/lib/microwaveprop/weather/hrrr_native_profile.ex +++ b/lib/microwaveprop/weather/hrrr_native_profile.ex @@ -80,7 +80,7 @@ defmodule Microwaveprop.Weather.HrrrNativeProfile do |> Enum.map(&{&1, get_field(changeset, &1)}) |> Enum.reject(fn {_k, v} -> is_nil(v) end) - expected = level_count || (arrays |> Enum.map(fn {_, v} -> length(v) end) |> Enum.max(fn -> 0 end)) + expected = level_count || arrays |> Enum.map(fn {_, v} -> length(v) end) |> Enum.max(fn -> 0 end) Enum.reduce(arrays, changeset, fn {key, values}, acc -> if length(values) == expected do diff --git a/lib/microwaveprop/weather/ncei_metar_client.ex b/lib/microwaveprop/weather/ncei_metar_client.ex index 150e8fda..e26371a6 100644 --- a/lib/microwaveprop/weather/ncei_metar_client.ex +++ b/lib/microwaveprop/weather/ncei_metar_client.ex @@ -56,8 +56,8 @@ defmodule Microwaveprop.Weather.NceiMetarClient do # Fixed-width header: WBAN(5) ICAO(4) FAA(4) DATE(8) HHMM(4) ... # Then METAR text after the "5-MIN" marker with <<_wban::binary-5, icao::binary-4, _faa::binary-4, date::binary-8, hhmm::binary-4, _rest::binary>> <- line, - {:ok, observed_at} <- parse_datetime(date, hhmm), - metar_text <- extract_metar(line) do + {:ok, observed_at} <- parse_datetime(date, hhmm) do + metar_text = extract_metar(line) obs = parse_metar_fields(metar_text) [Map.merge(obs, %{icao: String.trim(icao), observed_at: observed_at})] else @@ -106,30 +106,36 @@ defmodule Microwaveprop.Weather.NceiMetarClient do # Extract precise temp/dew from RMK T-group: T02110094 → 21.1°C / 9.4°C defp extract_temp_f(parts) do - with t_group when is_binary(t_group) <- Enum.find(parts, &String.starts_with?(&1, "T0")) do - parse_t_group_temp(t_group) - else + case Enum.find(parts, &String.starts_with?(&1, "T0")) do + t_group when is_binary(t_group) -> + parse_t_group_temp(t_group) + _ -> # Fall back to temp/dew field: "21/09" - with td when is_binary(td) <- Enum.find(parts, &Regex.match?(~r/^M?\d+\/M?\d+$/, &1)) do - [t, _d] = String.split(td, "/") - parse_metar_temp(t) |> c_to_f() - else - _ -> nil + case Enum.find(parts, &Regex.match?(~r/^M?\d+\/M?\d+$/, &1)) do + td when is_binary(td) -> + [t, _d] = String.split(td, "/") + t |> parse_metar_temp() |> c_to_f() + + _ -> + nil end end end defp extract_dewpoint_f(parts) do - with t_group when is_binary(t_group) <- Enum.find(parts, &String.starts_with?(&1, "T0")) do - parse_t_group_dew(t_group) - else + case Enum.find(parts, &String.starts_with?(&1, "T0")) do + t_group when is_binary(t_group) -> + parse_t_group_dew(t_group) + _ -> - with td when is_binary(td) <- Enum.find(parts, &Regex.match?(~r/^M?\d+\/M?\d+$/, &1)) do - [_t, d] = String.split(td, "/") - parse_metar_temp(d) |> c_to_f() - else - _ -> nil + case Enum.find(parts, &Regex.match?(~r/^M?\d+\/M?\d+$/, &1)) do + td when is_binary(td) -> + [_t, d] = String.split(td, "/") + d |> parse_metar_temp() |> c_to_f() + + _ -> + nil end end end @@ -137,14 +143,16 @@ defmodule Microwaveprop.Weather.NceiMetarClient do # T02110094 → temp = +21.1°C, dew = +9.4°C (leading 0=+, 1=-) defp parse_t_group_temp(<<"T", sign::binary-1, d1::binary-2, d2::binary-1, _rest::binary>>) do val = String.to_integer(d1) + String.to_integer(d2) / 10.0 - (if sign == "1", do: -val, else: val) |> c_to_f() + if_result = if sign == "1", do: -val, else: val + c_to_f(if_result) end defp parse_t_group_temp(_), do: nil defp parse_t_group_dew(<<"T", _::binary-4, sign::binary-1, d1::binary-2, d2::binary-1, _rest::binary>>) do val = String.to_integer(d1) + String.to_integer(d2) / 10.0 - (if sign == "1", do: -val, else: val) |> c_to_f() + if_result = if sign == "1", do: -val, else: val + c_to_f(if_result) end defp parse_t_group_dew(_), do: nil @@ -174,18 +182,20 @@ defmodule Microwaveprop.Weather.NceiMetarClient do end defp extract_altimeter(parts) do - with alt when is_binary(alt) <- Enum.find(parts, &String.starts_with?(&1, "A2")) do - case Float.parse(String.slice(alt, 1..-1//1)) do - {val, _} -> val / 100.0 - :error -> nil - end - else - _ -> nil + case Enum.find(parts, &String.starts_with?(&1, "A2")) do + alt when is_binary(alt) -> + case Float.parse(String.slice(alt, 1..-1//1)) do + {val, _} -> val / 100.0 + :error -> nil + end + + _ -> + nil end end defp extract_sky(parts) do sky_tokens = Enum.filter(parts, &Regex.match?(~r/^(CLR|FEW|SCT|BKN|OVC|VV)\d*/, &1)) - if sky_tokens != [], do: Enum.join(sky_tokens, " "), else: nil + if sky_tokens == [], do: nil, else: Enum.join(sky_tokens, " ") end end diff --git a/lib/microwaveprop/weather/nexrad_client.ex b/lib/microwaveprop/weather/nexrad_client.ex new file mode 100644 index 00000000..199df5aa --- /dev/null +++ b/lib/microwaveprop/weather/nexrad_client.ex @@ -0,0 +1,331 @@ +defmodule Microwaveprop.Weather.NexradClient do + @moduledoc """ + Fetches and processes IEM CONUS n0q composite reflectivity PNGs. + + The n0q product is a palettized 8-bit PNG, 12200 x 5400 pixels, + covering CONUS at 0.005 degrees/pixel (-126W to -65W, 23N to 50N). + Available every 5 minutes from the IEM archive. + + Pixel value 0 = no echo; values 1-255 map linearly to roughly + -30 to +95 dBZ. For each point of interest we extract a ~25 km + box (~50x50 pixels) and compute summary statistics. + """ + + require Logger + + # Image geometry + @lon_min -126.0 + @lat_max 50.0 + @deg_per_pixel 0.005 + + # Default box half-size in pixels (~25 km at mid-latitudes) + @default_box_half 25 + + @base_url "https://mesonet.agron.iastate.edu/archive/data" + + @doc "Fetch and process one n0q frame. Returns list of observation maps for points of interest." + @spec fetch_frame(DateTime.t(), [{float(), float()}]) :: {:ok, [map()]} | {:error, term()} + def fetch_frame(timestamp, points_of_interest) do + rounded = round_to_5min(timestamp) + url = frame_url(rounded) + + req_opts = Application.get_env(:microwaveprop, :nexrad_req_options, []) + + case Req.get(url, [receive_timeout: 60_000, retry: false] ++ req_opts) do + {:ok, %{status: 200, body: body}} when is_binary(body) -> + process_frame(body, rounded, points_of_interest) + + {:ok, %{status: status}} -> + {:error, "NEXRAD n0q HTTP #{status}"} + + {:error, reason} -> + {:error, reason} + end + end + + @doc "Round a DateTime down to the nearest 5-minute boundary." + @spec round_to_5min(DateTime.t()) :: DateTime.t() + def round_to_5min(%DateTime{} = dt) do + rounded_minute = div(dt.minute, 5) * 5 + %{dt | minute: rounded_minute, second: 0, microsecond: {0, 0}} + end + + @doc "Build the IEM archive URL for an n0q frame at the given (already-rounded) timestamp." + @spec frame_url(DateTime.t()) :: String.t() + def frame_url(%DateTime{} = dt) do + date_path = Calendar.strftime(dt, "%Y/%m/%d") + file_ts = Calendar.strftime(dt, "%Y%m%d%H%M") + "#{@base_url}/#{date_path}/GIS/uscomp/n0q_#{file_ts}.png" + end + + @doc """ + Convert (lat, lon) to pixel coordinates in the n0q image. + + Returns `{x, y}` where x is the column (0 = leftmost = -126W) + and y is the row (0 = topmost = 50N). + """ + @spec latlon_to_pixel(float(), float()) :: {integer(), integer()} + def latlon_to_pixel(lat, lon) do + x = round((lon - @lon_min) / @deg_per_pixel) + y = round((@lat_max - lat) / @deg_per_pixel) + {x, y} + end + + @doc """ + Convert a pixel value (0-255) to dBZ. + + Value 0 = no echo (returns 0.0). Values 1-255 map linearly to + approximately -30 to +95 dBZ. + """ + @spec pixel_to_dbz(integer()) :: float() + def pixel_to_dbz(0), do: 0.0 + + def pixel_to_dbz(value) when value >= 1 and value <= 255 do + # Linear mapping: 1 -> -30 dBZ, 255 -> 95 dBZ + -30.0 + (value - 1) / 254.0 * 125.0 + end + + @doc """ + Extract summary statistics from a box of pixels centered at (cx, cy). + + `pixels` is a flat binary where each byte is one pixel value. + `width` is the image width (number of columns per row). + `half` is the half-size of the extraction box in pixels. + + Returns a map with :mean_reflectivity_dbz, :max_reflectivity_dbz, + :texture_variance, and :pixel_count. + """ + @spec extract_box(binary(), integer(), integer(), integer(), integer()) :: map() + def extract_box(pixels, width, cx, cy, half \\ @default_box_half) do + height = div(byte_size(pixels), width) + + x_min = max(cx - half, 0) + x_max = min(cx + half, width - 1) + y_min = max(cy - half, 0) + y_max = min(cy + half, height - 1) + + # Collect non-zero dBZ values + dbz_values = + for y <- y_min..y_max, + x <- x_min..x_max, + offset = y * width + x, + offset >= 0, + offset < byte_size(pixels), + <<_::binary-size(offset), pixel_val::8, _::binary>> = pixels, + pixel_val > 0 do + pixel_to_dbz(pixel_val) + end + + count = length(dbz_values) + + compute_box_stats(dbz_values, count) + end + + # -- Private -- + + defp compute_box_stats(_values, 0) do + %{mean_reflectivity_dbz: 0.0, max_reflectivity_dbz: 0.0, texture_variance: 0.0, pixel_count: 0} + end + + defp compute_box_stats(values, count) do + mean = Enum.sum(values) / count + max_val = Enum.max(values) + variance = sample_variance(values, mean, count) + + %{ + mean_reflectivity_dbz: Float.round(mean, 2), + max_reflectivity_dbz: Float.round(max_val, 2), + texture_variance: Float.round(variance, 2), + pixel_count: count + } + end + + defp sample_variance(_values, _mean, count) when count < 2, do: 0.0 + + defp sample_variance(values, mean, count) do + values + |> Enum.map(fn v -> (v - mean) * (v - mean) end) + |> Enum.sum() + |> Kernel./(count - 1) + end + + defp process_frame(png_binary, observed_at, points) do + case decode_png_to_pixels(png_binary) do + {:ok, pixels, width} -> + observations = + Enum.map(points, fn {lat, lon} -> + {cx, cy} = latlon_to_pixel(lat, lon) + stats = extract_box(pixels, width, cx, cy) + + Map.merge(stats, %{ + observed_at: observed_at, + lat: lat, + lon: lon + }) + end) + + {:ok, observations} + + {:error, reason} -> + {:error, reason} + end + end + + # Minimal PNG decoder for indexed-color (palette) 8-bit images. + # + # Extracts raw palette index bytes (not RGB), which is what we need + # since the n0q product encodes dBZ as palette index values. + # + # PNG structure: + # 8-byte signature + # Chunks: [length(4) | type(4) | data(length) | crc(4)] + # IHDR chunk has width, height, bit_depth, color_type + # IDAT chunks contain zlib-compressed scanline data + # Each scanline: filter_byte | pixel_data[width] + defp decode_png_to_pixels(png_binary) do + # Validate PNG signature + <<137, 80, 78, 71, 13, 10, 26, 10, rest::binary>> = png_binary + + # Parse chunks + chunks = parse_chunks(rest, []) + + # Get dimensions from IHDR + [{_type, ihdr_data} | _] = chunks + <> = ihdr_data + + # Concatenate all IDAT data and decompress + idat_data = + chunks + |> Enum.filter(fn {type, _} -> type == "IDAT" end) + |> Enum.map(fn {_, data} -> data end) + |> IO.iodata_to_binary() + + decompressed = :zlib.uncompress(idat_data) + + # Undo per-row PNG filters, extract raw pixel index bytes + pixels = unfilter_rows(decompressed, width, height) + + {:ok, pixels, width} + rescue + e -> + {:error, "PNG decode failed: #{inspect(e)}"} + end + + defp parse_chunks(<<>>, acc), do: Enum.reverse(acc) + + defp parse_chunks(<>, acc) do + <> = rest + parse_chunks(remaining, [{type, data} | acc]) + end + + # Undo PNG row filters for 8-bit indexed images (1 byte per pixel). + # Returns a flat binary of pixel index values. + defp unfilter_rows(data, width, height) do + # bytes_per_pixel = 1 for indexed color + # Each row is: 1 filter byte + width pixel bytes + stride = width + + {rows, _} = + Enum.reduce(1..height, {[], data}, fn _row_idx, {acc, remaining} -> + <> = remaining + prev_row = List.first(acc, :binary.copy(<<0>>, stride)) + filtered = apply_filter(filter, row_data, prev_row, 1) + {[filtered | acc], rest} + end) + + rows + |> Enum.reverse() + |> IO.iodata_to_binary() + end + + # PNG filter types for 8-bit indexed (bpp=1) + defp apply_filter(0, row, _prev, _bpp), do: row + + defp apply_filter(1, row, _prev, bpp) do + # Sub: each byte is diff from byte bpp positions to the left + unfilter_sub(row, bpp) + end + + defp apply_filter(2, row, prev, _bpp) do + # Up: each byte is diff from byte directly above + unfilter_up(row, prev) + end + + defp apply_filter(3, row, prev, bpp) do + # Average: each byte is diff from average of left and above + unfilter_average(row, prev, bpp) + end + + defp apply_filter(4, row, prev, bpp) do + # Paeth: each byte uses Paeth predictor + unfilter_paeth(row, prev, bpp) + end + + defp unfilter_sub(row, bpp) do + row_bytes = :binary.bin_to_list(row) + + {result, _} = + Enum.reduce(row_bytes, {[], 0}, fn byte, {acc, idx} -> + left = if idx >= bpp, do: Enum.at(acc, idx - bpp), else: 0 + val = rem(byte + left, 256) + {acc ++ [val], idx + 1} + end) + + :binary.list_to_bin(result) + end + + defp unfilter_up(row, prev) do + row_bytes = :binary.bin_to_list(row) + prev_bytes = :binary.bin_to_list(prev) + + row_bytes + |> Enum.zip(prev_bytes) + |> Enum.map(fn {a, b} -> rem(a + b, 256) end) + |> :binary.list_to_bin() + end + + defp unfilter_average(row, prev, bpp) do + row_bytes = :binary.bin_to_list(row) + prev_bytes = :binary.bin_to_list(prev) + + {result, _} = + Enum.reduce(Enum.with_index(row_bytes), {[], prev_bytes}, fn {byte, idx}, {acc, prev_list} -> + left = if idx >= bpp, do: Enum.at(acc, idx - bpp), else: 0 + above = Enum.at(prev_list, idx, 0) + val = rem(byte + div(left + above, 2), 256) + {acc ++ [val], prev_list} + end) + + :binary.list_to_bin(result) + end + + defp unfilter_paeth(row, prev, bpp) do + row_bytes = :binary.bin_to_list(row) + prev_bytes = :binary.bin_to_list(prev) + + {result, _} = + Enum.reduce(Enum.with_index(row_bytes), {[], prev_bytes}, fn {byte, idx}, {acc, prev_list} -> + left = if idx >= bpp, do: Enum.at(acc, idx - bpp), else: 0 + above = Enum.at(prev_list, idx, 0) + upper_left = if idx >= bpp, do: Enum.at(prev_list, idx - bpp, 0), else: 0 + pr = paeth_predictor(left, above, upper_left) + val = rem(byte + pr, 256) + {acc ++ [val], prev_list} + end) + + :binary.list_to_bin(result) + end + + defp paeth_predictor(a, b, c) do + p = a + b - c + pa = abs(p - a) + pb = abs(p - b) + pc = abs(p - c) + + cond do + pa <= pb and pa <= pc -> a + pb <= pc -> b + true -> c + end + end +end diff --git a/lib/microwaveprop/weather/nexrad_observation.ex b/lib/microwaveprop/weather/nexrad_observation.ex new file mode 100644 index 00000000..4457c7e6 --- /dev/null +++ b/lib/microwaveprop/weather/nexrad_observation.ex @@ -0,0 +1,31 @@ +defmodule Microwaveprop.Weather.NexradObservation do + @moduledoc false + use Ecto.Schema + + import Ecto.Changeset + + @primary_key {:id, :binary_id, autogenerate: true} + @foreign_key_type :binary_id + + schema "nexrad_observations" do + field :observed_at, :utc_datetime + field :lat, :float + field :lon, :float + field :mean_reflectivity_dbz, :float + field :max_reflectivity_dbz, :float + field :texture_variance, :float + field :pixel_count, :integer + + timestamps(type: :utc_datetime) + end + + @required_fields ~w(observed_at lat lon)a + @optional_fields ~w(mean_reflectivity_dbz max_reflectivity_dbz texture_variance pixel_count)a + + def changeset(observation, attrs) do + observation + |> cast(attrs, @required_fields ++ @optional_fields) + |> validate_required(@required_fields) + |> unique_constraint([:lat, :lon, :observed_at]) + end +end diff --git a/lib/microwaveprop/workers/admin_task_worker.ex b/lib/microwaveprop/workers/admin_task_worker.ex index 852e0399..5ad07221 100644 --- a/lib/microwaveprop/workers/admin_task_worker.ex +++ b/lib/microwaveprop/workers/admin_task_worker.ex @@ -14,7 +14,7 @@ defmodule Microwaveprop.Workers.AdminTaskWorker do """ use Oban.Worker, queue: :admin, max_attempts: 1, unique: [period: 60] - require Logger + import Ecto.Query alias Microwaveprop.Backtest alias Microwaveprop.Backtest.Features @@ -24,7 +24,7 @@ defmodule Microwaveprop.Workers.AdminTaskWorker do alias Microwaveprop.Weather.HrrrNativeProfile alias Microwaveprop.Weather.ThetaE - import Ecto.Query + require Logger @impl Oban.Worker def perform(%Oban.Job{args: %{"task" => "backtest_all"} = args}) do @@ -48,7 +48,7 @@ defmodule Microwaveprop.Workers.AdminTaskWorker do sample_size = Map.get(args, "sample_size", 5000) fun_atom = String.to_atom(feature_name) - unless function_exported?(Features, fun_atom, 3) do + if !function_exported?(Features, fun_atom, 3) do Logger.error("AdminTask: unknown feature #{feature_name}") {:error, "unknown feature: #{feature_name}"} end @@ -183,6 +183,57 @@ defmodule Microwaveprop.Workers.AdminTaskWorker do :ok end + def perform(%Oban.Job{args: %{"task" => "recalibrate"} = args}) do + sample_size = Map.get(args, "sample_size", 5000) + epochs = Map.get(args, "epochs", 2000) + learning_rate = Map.get(args, "learning_rate", 0.01) + + Logger.info("AdminTask: running weight recalibration (sample=#{sample_size}, epochs=#{epochs}, lr=#{learning_rate})") + + result = + Microwaveprop.Propagation.Recalibrator.fit( + sample_size: sample_size, + epochs: epochs, + learning_rate: learning_rate + ) + + Logger.info("AdminTask: recalibration complete") + + Logger.info(" train_loss=#{result.train_loss}, val_loss=#{result.val_loss}, initial_loss=#{result.initial_loss}") + + for {factor, weight} <- Enum.sort(result.weights) do + Logger.info(" #{factor}: #{weight}") + end + + :ok + end + + def perform(%Oban.Job{args: %{"task" => "scorer_diff", "new_weights" => new_weights_json}}) do + alias Microwaveprop.Propagation.BandConfig + alias Microwaveprop.Propagation.ScorerDiff + + new_weights = Jason.decode!(new_weights_json) + old_weights = Map.new(BandConfig.weights(), fn {k, v} -> {to_string(k), v} end) + + Logger.info("AdminTask: running scorer diff comparison") + + result = ScorerDiff.compare(old_weights, new_weights) + markdown = ScorerDiff.to_markdown(result) + + path = "priv/scorer_diff_reports/latest.md" + File.mkdir_p!(Path.dirname(path)) + File.write!(path, markdown) + + Logger.info( + "AdminTask: scorer diff complete — #{result.summary.total} points, " <> + "#{result.summary.regressions} regressions, #{result.summary.improvements} improvements" + ) + + Logger.info("AdminTask: wrote report to #{path}") + + :ok + end + def perform(%Oban.Job{args: %{"task" => task}}) do Logger.error("AdminTask: unknown task #{task}") {:error, "unknown task: #{task}"} diff --git a/lib/microwaveprop/workers/hrrr_native_grid_worker.ex b/lib/microwaveprop/workers/hrrr_native_grid_worker.ex index eac6f4aa..e010ddb0 100644 --- a/lib/microwaveprop/workers/hrrr_native_grid_worker.ex +++ b/lib/microwaveprop/workers/hrrr_native_grid_worker.ex @@ -87,7 +87,8 @@ defmodule Microwaveprop.Workers.HrrrNativeGridWorker do |> where([p], p.valid_time == ^valid_time) |> where( [p], - fragment("ROW(?, ?) IN (SELECT * FROM UNNEST(?::float[], ?::float[]))", + fragment( + "ROW(?, ?) IN (SELECT * FROM UNNEST(?::float[], ?::float[]))", p.lat, p.lon, ^Enum.map(points, &elem(&1, 0)), @@ -112,14 +113,12 @@ defmodule Microwaveprop.Workers.HrrrNativeGridWorker do _ = Logger.info("HRRR native downloaded #{byte_size(grib_binary)} bytes, extracting #{length(points)} points..."), {:ok, profiles} <- HrrrNativeClient.extract_native_profiles(grib_binary, points) do rows = - profiles - |> Enum.flat_map(fn {{lat, lon}, profile} -> + Enum.flat_map(profiles, fn {{lat, lon}, profile} -> if profile[:level_count] && profile.level_count > 0 do now = DateTime.truncate(DateTime.utc_now(), :second) [ - profile - |> Map.merge(%{ + Map.merge(profile, %{ id: Ecto.UUID.bingenerate(), valid_time: valid_time, run_time: valid_time, diff --git a/lib/microwaveprop/workers/nexrad_worker.ex b/lib/microwaveprop/workers/nexrad_worker.ex new file mode 100644 index 00000000..3425e398 --- /dev/null +++ b/lib/microwaveprop/workers/nexrad_worker.ex @@ -0,0 +1,96 @@ +defmodule Microwaveprop.Workers.NexradWorker do + @moduledoc """ + Fetches one n0q NEXRAD composite frame and stores per-point + observations in `nexrad_observations`. + + Jobs are unique on `{year, month, day, hour, minute}` so backfill + sweeps that enqueue duplicate timestamps collapse automatically. + """ + + use Oban.Worker, + queue: :nexrad, + max_attempts: 3, + unique: [ + period: :infinity, + states: [:available, :scheduled, :executing, :retryable], + keys: [:year, :month, :day, :hour, :minute] + ] + + import Ecto.Query + + alias Microwaveprop.Radio.Contact + alias Microwaveprop.Repo + alias Microwaveprop.Weather.NexradClient + alias Microwaveprop.Weather.NexradObservation + + require Logger + + @impl Oban.Worker + def perform(%Oban.Job{args: args}) do + %{"year" => year, "month" => month, "day" => day, "hour" => hour, "minute" => minute} = args + {:ok, date} = Date.new(year, month, day) + {:ok, time} = Time.new(hour, minute, 0) + {:ok, timestamp} = DateTime.new(date, time, "Etc/UTC") + + points = points_of_interest_for_time(timestamp) + + if points == [] do + Logger.info("NexradWorker: no points for #{timestamp}, skipping") + :ok + else + fetch_and_upsert(timestamp, points) + end + end + + @doc false + def points_of_interest_for_time(timestamp) do + # Match contacts within +/- 30 minutes of this frame + time_start = DateTime.add(timestamp, -1800, :second) + time_end = DateTime.add(timestamp, 1800, :second) + + Contact + |> where([c], not is_nil(c.pos1)) + |> where([c], c.qso_timestamp >= ^time_start and c.qso_timestamp <= ^time_end) + |> select([c], c.pos1) + |> Repo.all() + |> Enum.flat_map(fn pos -> + case {pos["lat"], pos["lon"] || pos["lng"]} do + {lat, lon} when is_number(lat) and is_number(lon) -> [{snap(lat), snap(lon)}] + _ -> [] + end + end) + |> Enum.uniq() + end + + defp snap(x), do: Float.round(x * 1.0, 3) + + defp fetch_and_upsert(timestamp, points) do + case NexradClient.fetch_frame(timestamp, points) do + {:ok, observations} -> + now = DateTime.truncate(DateTime.utc_now(), :second) + + rows = + Enum.map(observations, fn obs -> + obs + |> Map.put(:id, Ecto.UUID.bingenerate()) + |> Map.put(:inserted_at, now) + |> Map.put(:updated_at, now) + end) + + {inserted, _} = + Repo.insert_all( + NexradObservation, + rows, + on_conflict: {:replace_all_except, [:id, :inserted_at]}, + conflict_target: [:lat, :lon, :observed_at] + ) + + Logger.info("NexradWorker: upserted #{inserted} observations for #{timestamp}") + :ok + + {:error, reason} -> + Logger.warning("NexradWorker failed for #{timestamp}: #{inspect(reason)}") + {:error, reason} + end + end +end diff --git a/lib/microwaveprop_web/router.ex b/lib/microwaveprop_web/router.ex index f4bb5992..fe526f6d 100644 --- a/lib/microwaveprop_web/router.ex +++ b/lib/microwaveprop_web/router.ex @@ -104,7 +104,7 @@ defmodule MicrowavepropWeb.Router do ecto_psql_extras_options: [long_running_queries: [threshold: "200 milliseconds"]], on_mount: [{MicrowavepropWeb.UserAuth, :require_admin}] - oban_dashboard "/admin/oban", on_mount: [{MicrowavepropWeb.UserAuth, :require_admin}] + oban_dashboard("/admin/oban", on_mount: [{MicrowavepropWeb.UserAuth, :require_admin}]) end # Enable Swoosh mailbox preview in development diff --git a/lib/mix/tasks/backtest.ex b/lib/mix/tasks/backtest.ex index 2ce05374..08a873e3 100644 --- a/lib/mix/tasks/backtest.ex +++ b/lib/mix/tasks/backtest.ex @@ -29,6 +29,7 @@ defmodule Mix.Tasks.Backtest do use Mix.Task alias Microwaveprop.Backtest + alias Microwaveprop.Backtest.Features @impl Mix.Task def run(argv) do @@ -52,7 +53,7 @@ defmodule Mix.Tasks.Backtest do baseline_size = Keyword.get(opts, :baseline, sample_size) out_path = Keyword.get(opts, :out) - features = Microwaveprop.Backtest.Features.all_features() + features = Features.all_features() Mix.shell().info("Running consolidated backtest for #{map_size(features)} features...") results = @@ -104,8 +105,8 @@ defmodule Mix.Tasks.Backtest do defp resolve_feature(spec) do case String.split(spec, ".") do [name] -> - fun = resolve_function_atom!(Microwaveprop.Backtest.Features, name) - feature_fun = &apply(Microwaveprop.Backtest.Features, fun, [&1, &2, &3]) + fun = resolve_function_atom!(Features, name) + feature_fun = &apply(Features, fun, [&1, &2, &3]) {feature_fun, "Microwaveprop.Backtest.Features.#{fun}"} parts -> @@ -127,13 +128,12 @@ defmodule Mix.Tasks.Backtest do fun else exported = - module.__info__(:functions) + :functions + |> module.__info__() |> Enum.filter(fn {_f, arity} -> arity == 3 end) |> Enum.map_join(", ", fn {f, _} -> to_string(f) end) - Mix.raise( - "Feature #{inspect(module)}.#{normalized}/3 is not defined.\nAvailable 3-arity functions: #{exported}" - ) + Mix.raise("Feature #{inspect(module)}.#{normalized}/3 is not defined.\nAvailable 3-arity functions: #{exported}") end end end diff --git a/lib/mix/tasks/hrrr_climatology.ex b/lib/mix/tasks/hrrr_climatology.ex index 553218df..40edcdcd 100644 --- a/lib/mix/tasks/hrrr_climatology.ex +++ b/lib/mix/tasks/hrrr_climatology.ex @@ -37,9 +37,7 @@ defmodule Mix.Tasks.HrrrClimatology do timeout: 120_000 ) - Mix.shell().info( - "Building climatology (min_samples=#{min_samples}, #{length(combos)} batches)..." - ) + Mix.shell().info("Building climatology (min_samples=#{min_samples}, #{length(combos)} batches)...") total = combos diff --git a/lib/mix/tasks/hrrr_native_backfill.ex b/lib/mix/tasks/hrrr_native_backfill.ex index b03084a8..8f5e9a70 100644 --- a/lib/mix/tasks/hrrr_native_backfill.ex +++ b/lib/mix/tasks/hrrr_native_backfill.ex @@ -47,24 +47,25 @@ defmodule Mix.Tasks.HrrrNativeBackfill do defp pad(n), do: n |> Integer.to_string() |> String.pad_leading(2, "0") defp top_hours_by_contact_count(limit) do - from(c in Contact, - where: not is_nil(c.pos1), - select: %{ - year: fragment("EXTRACT(YEAR FROM ?)::int", c.qso_timestamp), - month: fragment("EXTRACT(MONTH FROM ?)::int", c.qso_timestamp), - day: fragment("EXTRACT(DAY FROM ?)::int", c.qso_timestamp), - hour: fragment("EXTRACT(HOUR FROM ?)::int", c.qso_timestamp), - contacts: count(c.id) - }, - group_by: [ - fragment("EXTRACT(YEAR FROM ?)::int", c.qso_timestamp), - fragment("EXTRACT(MONTH FROM ?)::int", c.qso_timestamp), - fragment("EXTRACT(DAY FROM ?)::int", c.qso_timestamp), - fragment("EXTRACT(HOUR FROM ?)::int", c.qso_timestamp) - ], - order_by: [desc: count(c.id)], - limit: ^limit + Repo.all( + from(c in Contact, + where: not is_nil(c.pos1), + select: %{ + year: fragment("EXTRACT(YEAR FROM ?)::int", c.qso_timestamp), + month: fragment("EXTRACT(MONTH FROM ?)::int", c.qso_timestamp), + day: fragment("EXTRACT(DAY FROM ?)::int", c.qso_timestamp), + hour: fragment("EXTRACT(HOUR FROM ?)::int", c.qso_timestamp), + contacts: count(c.id) + }, + group_by: [ + fragment("EXTRACT(YEAR FROM ?)::int", c.qso_timestamp), + fragment("EXTRACT(MONTH FROM ?)::int", c.qso_timestamp), + fragment("EXTRACT(DAY FROM ?)::int", c.qso_timestamp), + fragment("EXTRACT(HOUR FROM ?)::int", c.qso_timestamp) + ], + order_by: [desc: count(c.id)], + limit: ^limit + ) ) - |> Repo.all() end end diff --git a/lib/mix/tasks/import_contest_logs.ex b/lib/mix/tasks/import_contest_logs.ex index 5c67c08c..e9ab1f75 100644 --- a/lib/mix/tasks/import_contest_logs.ex +++ b/lib/mix/tasks/import_contest_logs.ex @@ -41,7 +41,7 @@ defmodule Mix.Tasks.ImportContestLogs do |> List.flatten() |> Enum.join("\n") - total_lines = combined |> String.split("\n") |> length() |> then(&(&1 - 1)) + total_lines = combined |> String.split("\n") |> length() |> Kernel.-(1) Mix.shell().info("Combined #{length(argv)} files: #{total_lines} data rows") Mix.shell().info("Running preview (validate + deduplicate)...") diff --git a/lib/mix/tasks/nexrad_backfill.ex b/lib/mix/tasks/nexrad_backfill.ex new file mode 100644 index 00000000..291b53ef --- /dev/null +++ b/lib/mix/tasks/nexrad_backfill.ex @@ -0,0 +1,72 @@ +defmodule Mix.Tasks.NexradBackfill do + @shortdoc "Enqueue NexradWorker jobs for the top-N hours by contact count" + @moduledoc """ + Backfills the `nexrad_observations` table by enqueueing one + `NexradWorker` job per distinct `(year, month, day, hour, minute=0)` + where we have contacts, prioritized by contact count so the most + data-dense hours land first. + + Each n0q frame is 2-4 MB, so a 200-hour backfill is ~400-800 MB. + + mix nexrad_backfill --limit 200 + + Jobs are deduplicated by Oban's unique constraint on + `{year, month, day, hour, minute}` so running the task twice is safe. + """ + use Mix.Task + + import Ecto.Query + + alias Microwaveprop.Radio.Contact + alias Microwaveprop.Repo + alias Microwaveprop.Workers.NexradWorker + + @impl Mix.Task + def run(argv) do + Mix.Task.run("app.start") + Oban.pause_all_queues(Oban) + + {opts, _, _} = OptionParser.parse(argv, switches: [limit: :integer]) + limit = Keyword.get(opts, :limit, 200) + + hours = top_hours_by_contact_count(limit) + Mix.shell().info("Enqueueing #{length(hours)} NexradWorker jobs") + + Enum.each(hours, fn %{year: y, month: m, day: d, hour: h, contacts: n} -> + args = %{"year" => y, "month" => m, "day" => d, "hour" => h, "minute" => 0} + + case Oban.insert(NexradWorker.new(args)) do + {:ok, _job} -> + Mix.shell().info(" #{y}-#{pad(m)}-#{pad(d)} #{pad(h)}Z (#{n} contacts)") + + {:error, reason} -> + Mix.shell().error(" #{y}-#{pad(m)}-#{pad(d)} #{pad(h)}Z failed: #{inspect(reason)}") + end + end) + end + + defp pad(n), do: n |> Integer.to_string() |> String.pad_leading(2, "0") + + defp top_hours_by_contact_count(limit) do + Repo.all( + from(c in Contact, + where: not is_nil(c.pos1), + select: %{ + year: fragment("EXTRACT(YEAR FROM ?)::int", c.qso_timestamp), + month: fragment("EXTRACT(MONTH FROM ?)::int", c.qso_timestamp), + day: fragment("EXTRACT(DAY FROM ?)::int", c.qso_timestamp), + hour: fragment("EXTRACT(HOUR FROM ?)::int", c.qso_timestamp), + contacts: count(c.id) + }, + group_by: [ + fragment("EXTRACT(YEAR FROM ?)::int", c.qso_timestamp), + fragment("EXTRACT(MONTH FROM ?)::int", c.qso_timestamp), + fragment("EXTRACT(DAY FROM ?)::int", c.qso_timestamp), + fragment("EXTRACT(HOUR FROM ?)::int", c.qso_timestamp) + ], + order_by: [desc: count(c.id)], + limit: ^limit + ) + ) + end +end diff --git a/lib/mix/tasks/recalibrate_scorer.ex b/lib/mix/tasks/recalibrate_scorer.ex new file mode 100644 index 00000000..1efe08bd --- /dev/null +++ b/lib/mix/tasks/recalibrate_scorer.ex @@ -0,0 +1,92 @@ +defmodule Mix.Tasks.RecalibrateScorer do + @shortdoc "Recalibrate propagation scorer weights via gradient descent" + @moduledoc """ + Fits new scoring weights using logistic regression on the historical + QSO corpus. Compares QSO conditions (positive examples) against + random-time baselines (negative examples) to learn which factors + best discriminate actual propagation events. + + ## Usage + + mix recalibrate_scorer + mix recalibrate_scorer --sample 5000 --epochs 2000 --lr 0.01 + + ## Options + + * `--sample` - max contacts to load (default: 5000) + * `--epochs` - training iterations (default: 2000) + * `--lr` - learning rate (default: 0.01) + """ + use Mix.Task + + alias Microwaveprop.Propagation.BandConfig + alias Microwaveprop.Propagation.Recalibrator + + @impl Mix.Task + def run(argv) do + Mix.Task.run("app.start") + Oban.pause_all_queues(Oban) + + {opts, _, _} = + OptionParser.parse(argv, + switches: [sample: :integer, epochs: :integer, lr: :float] + ) + + sample_size = Keyword.get(opts, :sample, 5000) + epochs = Keyword.get(opts, :epochs, 2000) + learning_rate = Keyword.get(opts, :lr, 0.01) + + Mix.shell().info("Recalibrating scorer weights...") + Mix.shell().info(" sample_size: #{sample_size}") + Mix.shell().info(" epochs: #{epochs}") + Mix.shell().info(" learning_rate: #{learning_rate}") + Mix.shell().info("") + + result = + Recalibrator.fit( + sample_size: sample_size, + epochs: epochs, + learning_rate: learning_rate + ) + + current_weights = BandConfig.weights() + + Mix.shell().info("=== Results ===") + Mix.shell().info("") + Mix.shell().info("Train loss: #{format_float(result.train_loss)}") + Mix.shell().info("Val loss: #{format_float(result.val_loss)}") + Mix.shell().info("Initial loss: #{format_float(result.initial_loss)}") + Mix.shell().info("") + + Mix.shell().info("=== Weight Comparison ===") + Mix.shell().info("") + + header = + String.pad_trailing("Factor", 18) <> String.pad_trailing("Current", 10) <> String.pad_trailing("New", 10) <> "Delta" + + Mix.shell().info(header) + Mix.shell().info(String.duplicate("-", 48)) + + for key <- result.weights |> Map.keys() |> Enum.sort() do + current = Map.get(current_weights, key, 0.0) + new_val = result.weights[key] + delta = new_val - current + + sign = if delta >= 0, do: "+", else: "" + + line = + String.pad_trailing(to_string(key), 18) <> + String.pad_trailing(format_float(current), 10) <> + String.pad_trailing(format_float(new_val), 10) <> + "#{sign}#{format_float(delta)}" + + Mix.shell().info(line) + end + + Mix.shell().info("") + Mix.shell().info("To apply these weights, update @weights in lib/microwaveprop/propagation/band_config.ex") + end + + defp format_float(f) when is_float(f), do: :erlang.float_to_binary(f, decimals: 4) + defp format_float(f), do: to_string(f) +end diff --git a/lib/mix/tasks/scorer_diff.ex b/lib/mix/tasks/scorer_diff.ex new file mode 100644 index 00000000..6deb3467 --- /dev/null +++ b/lib/mix/tasks/scorer_diff.ex @@ -0,0 +1,94 @@ +defmodule Mix.Tasks.ScorerDiff do + @shortdoc "Compare two weight sets on recent propagation scores" + @moduledoc """ + Loads the most recent HRRR frame's propagation scores and compares + the current weights against a new set, printing a markdown diff report. + + ## Usage + + mix scorer_diff --new-weights '{"humidity":0.18,"time_of_day":0.20,"td_depression":0.10,"refractivity":0.08,"sky":0.08,"season":0.08,"wind":0.05,"rain":0.08,"pressure":0.05,"pwat":0.10}' + + ## Options + + * `--new-weights` — JSON string of new weights (required) + * `--out` — optional file path to write the report + """ + use Mix.Task + + alias Microwaveprop.Propagation.BandConfig + alias Microwaveprop.Propagation.ScorerDiff + + @impl Mix.Task + def run(argv) do + Mix.Task.run("app.start") + Oban.pause_all_queues(Oban) + + {opts, _, _} = + OptionParser.parse(argv, + switches: [new_weights: :string, out: :string] + ) + + new_weights_json = Keyword.fetch!(opts, :new_weights) + out_path = Keyword.get(opts, :out) + + new_weights = Jason.decode!(new_weights_json) + validate_weights!(new_weights) + + old_weights = stringify_weights(BandConfig.weights()) + + Mix.shell().info("Comparing current weights vs new weights...") + Mix.shell().info("Old: #{inspect(old_weights)}") + Mix.shell().info("New: #{inspect(new_weights)}") + Mix.shell().info("") + + result = ScorerDiff.compare(old_weights, new_weights) + markdown = ScorerDiff.to_markdown(result) + + IO.puts(markdown) + + if out_path do + File.mkdir_p!(Path.dirname(out_path)) + File.write!(out_path, markdown) + Mix.shell().info("Wrote report to #{out_path}") + end + end + + defp stringify_weights(weights) do + Map.new(weights, fn {k, v} -> {to_string(k), v} end) + end + + defp validate_weights!(weights) do + validate_weight_keys!(weights) + validate_weight_sum!(weights) + end + + defp validate_weight_keys!(weights) do + expected_keys = + MapSet.new(~w(humidity time_of_day td_depression refractivity sky season wind rain pressure pwat)) + + provided_keys = MapSet.new(Map.keys(weights)) + + missing = MapSet.difference(expected_keys, provided_keys) + extra = MapSet.difference(provided_keys, expected_keys) + + errors = + [] + |> prepend_if(MapSet.size(missing) > 0, "Missing: #{inspect(MapSet.to_list(missing))}") + |> prepend_if(MapSet.size(extra) > 0, "Extra: #{inspect(MapSet.to_list(extra))}") + + if errors != [] do + Mix.raise("Invalid weight keys. " <> Enum.join(errors, " ")) + end + end + + defp validate_weight_sum!(weights) do + sum = weights |> Map.values() |> Enum.sum() |> Float.round(4) + + if !(sum >= 0.99 and sum <= 1.01) do + Mix.raise("Weights must sum to ~1.0, got #{sum}") + end + end + + defp prepend_if(list, true, item), do: list ++ [item] + defp prepend_if(list, false, _item), do: list +end diff --git a/priv/repo/migrations/20260410200000_create_nexrad_observations.exs b/priv/repo/migrations/20260410200000_create_nexrad_observations.exs new file mode 100644 index 00000000..23e383ff --- /dev/null +++ b/priv/repo/migrations/20260410200000_create_nexrad_observations.exs @@ -0,0 +1,21 @@ +defmodule Microwaveprop.Repo.Migrations.CreateNexradObservations do + use Ecto.Migration + + def change do + create table(:nexrad_observations, primary_key: false) do + add :id, :binary_id, primary_key: true, null: false + add :observed_at, :utc_datetime, null: false + add :lat, :float, null: false + add :lon, :float, null: false + add :mean_reflectivity_dbz, :float + add :max_reflectivity_dbz, :float + add :texture_variance, :float + add :pixel_count, :integer + + timestamps(type: :utc_datetime) + end + + create unique_index(:nexrad_observations, [:lat, :lon, :observed_at]) + create index(:nexrad_observations, [:observed_at]) + end +end diff --git a/test/microwaveprop/propagation/recalibrator_test.exs b/test/microwaveprop/propagation/recalibrator_test.exs new file mode 100644 index 00000000..bbf4b2e6 --- /dev/null +++ b/test/microwaveprop/propagation/recalibrator_test.exs @@ -0,0 +1,145 @@ +defmodule Microwaveprop.Propagation.RecalibratorTest do + use Microwaveprop.DataCase, async: true + + alias Microwaveprop.Propagation.Recalibrator + alias Microwaveprop.Radio.Contact + alias Microwaveprop.Weather.HrrrProfile + + @factor_keys ~w(humidity time_of_day td_depression refractivity sky season wind rain pwat pressure)a + + defp create_contact(attrs) do + ts = Map.get(attrs, :qso_timestamp, ~U[2024-07-15 06:00:00Z]) + lat = Map.get(attrs, :lat, 32.9) + lon = Map.get(attrs, :lon, -97.0) + + contact_attrs = %{ + station1: Map.fetch!(attrs, :station1), + station2: "K5TR", + qso_timestamp: ts, + mode: "CW", + band: Decimal.new("10000"), + grid1: "EM12", + grid2: "EM00", + pos1: %{"lat" => lat, "lon" => lon}, + pos2: %{"lat" => 30.3, "lon" => -97.7}, + distance_km: Decimal.new("295") + } + + {:ok, contact} = + %Contact{} + |> Contact.changeset(contact_attrs) + |> Repo.insert() + + contact + end + + defp create_hrrr_profile(attrs) do + {:ok, profile} = + %HrrrProfile{} + |> HrrrProfile.changeset(%{ + valid_time: Map.fetch!(attrs, :valid_time), + lat: Map.fetch!(attrs, :lat), + lon: Map.fetch!(attrs, :lon), + surface_temp_c: Map.get(attrs, :surface_temp_c, 28.0), + surface_dewpoint_c: Map.get(attrs, :surface_dewpoint_c, 22.0), + surface_pressure_mb: Map.get(attrs, :surface_pressure_mb, 1005.0), + min_refractivity_gradient: Map.get(attrs, :min_refractivity_gradient, -120.0), + hpbl_m: Map.get(attrs, :hpbl_m, 800.0), + pwat_mm: Map.get(attrs, :pwat_mm, 35.0) + }) + |> Repo.insert() + + profile + end + + # Synthetic factor vectors for training tests. + # Positives: high scores (good propagation conditions). + # Negatives: low scores (poor propagation conditions). + defp synthetic_positives do + for _ <- 1..20 do + [85, 90, 80, 75, 88, 70, 95, 100, 75, 82] + |> Enum.map(fn base -> base + :rand.uniform(10) - 5 end) + |> Enum.map(&min(100, max(0, &1))) + end + end + + defp synthetic_negatives do + for _ <- 1..20 do + [40, 35, 45, 50, 30, 55, 50, 60, 50, 45] + |> Enum.map(fn base -> base + :rand.uniform(10) - 5 end) + |> Enum.map(&min(100, max(0, &1))) + end + end + + describe "train/3" do + test "returns a weights map with all 10 factor keys summing to ~1.0" do + result = Recalibrator.train(synthetic_positives(), synthetic_negatives(), epochs: 50, learning_rate: 0.01) + + assert is_map(result.weights) + + for key <- @factor_keys do + assert Map.has_key?(result.weights, key), "missing weight for #{key}" + weight = result.weights[key] + assert is_float(weight), "weight for #{key} should be float, got #{inspect(weight)}" + assert weight >= 0.0, "weight for #{key} should be non-negative" + end + + sum = result.weights |> Map.values() |> Enum.sum() + assert_in_delta sum, 1.0, 0.001, "weights should sum to 1.0, got #{sum}" + end + + test "train_loss decreases from initial loss" do + result = Recalibrator.train(synthetic_positives(), synthetic_negatives(), epochs: 200, learning_rate: 0.01) + + assert is_float(result.train_loss) + assert is_float(result.val_loss) + assert is_float(result.initial_loss) + + assert result.train_loss < result.initial_loss, + "train_loss (#{result.train_loss}) should be less than initial_loss (#{result.initial_loss})" + end + end + + describe "fit/1" do + test "returns valid result with insufficient HRRR data (falls back to current weights)" do + # Contacts without matching HRRR profiles for baselines + create_contact(%{station1: "W5AA"}) + create_contact(%{station1: "W5BB", lat: 33.0, lon: -96.0, qso_timestamp: ~U[2024-07-16 06:00:00Z]}) + + result = Recalibrator.fit(sample_size: 5, epochs: 20, learning_rate: 0.01) + + assert is_map(result.weights) + assert map_size(result.weights) == 10 + + sum = result.weights |> Map.values() |> Enum.sum() + assert_in_delta sum, 1.0, 0.01 + end + end + + describe "compute_factors/2" do + test "builds a factor vector from an HRRR profile and timestamp" do + profile = + create_hrrr_profile(%{ + valid_time: ~U[2024-07-15 06:00:00Z], + lat: 32.9, + lon: -97.0, + surface_temp_c: 28.0, + surface_dewpoint_c: 22.0, + surface_pressure_mb: 1005.0, + min_refractivity_gradient: -120.0, + hpbl_m: 800.0, + pwat_mm: 35.0 + }) + + factors = Recalibrator.compute_factors(profile, ~U[2024-07-15 06:00:00Z]) + + assert is_list(factors) + assert length(factors) == 10 + + Enum.each(factors, fn f -> + assert is_number(f) + assert f >= 0 and f <= 100 + end) + end + end +end diff --git a/test/microwaveprop/propagation/scorer_diff_test.exs b/test/microwaveprop/propagation/scorer_diff_test.exs new file mode 100644 index 00000000..62c53683 --- /dev/null +++ b/test/microwaveprop/propagation/scorer_diff_test.exs @@ -0,0 +1,207 @@ +defmodule Microwaveprop.Propagation.ScorerDiffTest do + use Microwaveprop.DataCase, async: true + + alias Microwaveprop.Propagation.GridScore + alias Microwaveprop.Propagation.ScorerDiff + + @valid_time ~U[2026-04-10 12:00:00Z] + + @old_weights %{ + "humidity" => 0.18, + "time_of_day" => 0.10, + "td_depression" => 0.10, + "refractivity" => 0.08, + "sky" => 0.08, + "season" => 0.08, + "wind" => 0.05, + "rain" => 0.08, + "pressure" => 0.15, + "pwat" => 0.10 + } + + # Shift weight from pressure to time_of_day + @new_weights %{ + "humidity" => 0.18, + "time_of_day" => 0.20, + "td_depression" => 0.10, + "refractivity" => 0.08, + "sky" => 0.08, + "season" => 0.08, + "wind" => 0.05, + "rain" => 0.08, + "pressure" => 0.05, + "pwat" => 0.10 + } + + defp insert_score!(attrs) do + %GridScore{} + |> GridScore.changeset(attrs) + |> Repo.insert!() + end + + defp make_factors(overrides \\ %{}) do + base = %{ + "humidity" => 80, + "time_of_day" => 60, + "td_depression" => 70, + "refractivity" => 55, + "sky" => 88, + "season" => 65, + "wind" => 90, + "rain" => 95, + "pressure" => 40, + "pwat" => 75 + } + + Map.merge(base, overrides) + end + + setup do + # Insert scores at two different bands for the same valid_time + insert_score!(%{ + lat: 32.875, + lon: -97.0, + valid_time: @valid_time, + band_mhz: 10_000, + score: 70, + factors: make_factors() + }) + + insert_score!(%{ + lat: 33.0, + lon: -97.125, + valid_time: @valid_time, + band_mhz: 10_000, + score: 50, + factors: make_factors(%{"time_of_day" => 20, "pressure" => 90}) + }) + + insert_score!(%{ + lat: 32.875, + lon: -97.0, + valid_time: @valid_time, + band_mhz: 24_000, + score: 65, + factors: make_factors(%{"humidity" => 40, "season" => 90}) + }) + + # Insert an older score that should not be included + insert_score!(%{ + lat: 32.875, + lon: -97.0, + valid_time: ~U[2026-04-10 11:00:00Z], + band_mhz: 10_000, + score: 55, + factors: make_factors(%{"humidity" => 30}) + }) + + :ok + end + + describe "compare/3" do + test "returns summary with mean_diff, max_diff, regressions, improvements" do + result = ScorerDiff.compare(@old_weights, @new_weights) + + assert is_map(result.summary) + assert is_float(result.summary.mean_diff) + assert is_number(result.summary.max_diff) + assert is_integer(result.summary.regressions) + assert is_integer(result.summary.improvements) + end + + test "returns band_diffs keyed by band_mhz" do + result = ScorerDiff.compare(@old_weights, @new_weights) + + assert Map.has_key?(result.band_diffs, 10_000) + assert Map.has_key?(result.band_diffs, 24_000) + assert is_float(result.band_diffs[10_000].mean_diff) + end + + test "returns worst_regressions as a list" do + result = ScorerDiff.compare(@old_weights, @new_weights) + + assert is_list(result.worst_regressions) + + for entry <- result.worst_regressions do + assert Map.has_key?(entry, :lat) + assert Map.has_key?(entry, :lon) + assert Map.has_key?(entry, :band_mhz) + assert Map.has_key?(entry, :old_score) + assert Map.has_key?(entry, :new_score) + end + end + + test "computes correct scores from factors and weights" do + result = ScorerDiff.compare(@old_weights, @new_weights) + + # For the first 10 GHz point (lat 32.875, lon -97.0): + # factors: humidity=80, time_of_day=60, td_depression=70, refractivity=55, + # sky=88, season=65, wind=90, rain=95, pressure=40, pwat=75 + # + # old: 80*0.18 + 60*0.10 + 70*0.10 + 55*0.08 + 88*0.08 + 65*0.08 + 90*0.05 + 95*0.08 + 40*0.15 + 75*0.10 + # = 14.4 + 6.0 + 7.0 + 4.4 + 7.04 + 5.2 + 4.5 + 7.6 + 6.0 + 7.5 = 69.64 -> 70 + # + # new: 80*0.18 + 60*0.20 + 70*0.10 + 55*0.08 + 88*0.08 + 65*0.08 + 90*0.05 + 95*0.08 + 40*0.05 + 75*0.10 + # = 14.4 + 12.0 + 7.0 + 4.4 + 7.04 + 5.2 + 4.5 + 7.6 + 2.0 + 7.5 = 71.64 -> 72 + # + # diff = 72 - 70 = +2 (improvement) + + # For the second 10 GHz point (lat 33.0, lon -97.125): + # factors: time_of_day=20, pressure=90 (others same) + # + # old: 80*0.18 + 20*0.10 + 70*0.10 + 55*0.08 + 88*0.08 + 65*0.08 + 90*0.05 + 95*0.08 + 90*0.15 + 75*0.10 + # = 14.4 + 2.0 + 7.0 + 4.4 + 7.04 + 5.2 + 4.5 + 7.6 + 13.5 + 7.5 = 73.14 -> 73 + # + # new: 80*0.18 + 20*0.20 + 70*0.10 + 55*0.08 + 88*0.08 + 65*0.08 + 90*0.05 + 95*0.08 + 90*0.05 + 75*0.10 + # = 14.4 + 4.0 + 7.0 + 4.4 + 7.04 + 5.2 + 4.5 + 7.6 + 4.5 + 7.5 = 66.14 -> 66 + # + # diff = 66 - 73 = -7 (regression!) + + # We should see at least 1 regression (diff > 5 in magnitude) + assert result.summary.regressions >= 1 + + # Total count should be 3 (only the most recent valid_time) + assert result.summary.total == 3 + end + + test "identical weights produce zero diff" do + result = ScorerDiff.compare(@old_weights, @old_weights) + + assert result.summary.mean_diff == 0.0 + assert result.summary.max_diff == 0 + assert result.summary.regressions == 0 + assert result.summary.improvements == 0 + end + + test "only loads scores from the most recent valid_time" do + result = ScorerDiff.compare(@old_weights, @new_weights) + + # 3 scores at @valid_time (2 for 10GHz, 1 for 24GHz), 1 at earlier time + assert result.summary.total == 3 + end + + test "accepts limit option" do + result = ScorerDiff.compare(@old_weights, @new_weights, limit: 1) + + # Should only process 1 score + assert result.summary.total == 1 + end + end + + describe "weighted_sum/2" do + test "computes weighted sum from factors and weights" do + factors = %{"humidity" => 100, "time_of_day" => 50} + weights = %{"humidity" => 0.5, "time_of_day" => 0.5} + + assert ScorerDiff.weighted_sum(factors, weights) == 75 + end + + test "rounds to nearest integer" do + factors = %{"humidity" => 80, "time_of_day" => 60} + weights = %{"humidity" => 0.6, "time_of_day" => 0.4} + + # 80*0.6 + 60*0.4 = 48 + 24 = 72 + assert ScorerDiff.weighted_sum(factors, weights) == 72 + end + end +end diff --git a/test/microwaveprop/weather/hrrr_native_client_test.exs b/test/microwaveprop/weather/hrrr_native_client_test.exs index 7acf07f8..ad851d70 100644 --- a/test/microwaveprop/weather/hrrr_native_client_test.exs +++ b/test/microwaveprop/weather/hrrr_native_client_test.exs @@ -34,7 +34,7 @@ defmodule Microwaveprop.Weather.HrrrNativeClientTest do # Pick level 5 — should have all 7 vars level5 = Enum.filter(messages, fn %{level: l} -> l == "5 hybrid level" end) - vars = Enum.map(level5, & &1.var) |> Enum.sort() + vars = level5 |> Enum.map(& &1.var) |> Enum.sort() assert vars == ~w(HGT PRES SPFH TKE TMP UGRD VGRD) end diff --git a/test/microwaveprop/weather/hrrr_native_profile_test.exs b/test/microwaveprop/weather/hrrr_native_profile_test.exs index 3cec335b..3b863354 100644 --- a/test/microwaveprop/weather/hrrr_native_profile_test.exs +++ b/test/microwaveprop/weather/hrrr_native_profile_test.exs @@ -43,7 +43,7 @@ defmodule Microwaveprop.Weather.HrrrNativeProfileTest do attrs = Map.put(@valid_attrs, :temp_k, [295.0, 292.0]) changeset = HrrrNativeProfile.changeset(%HrrrNativeProfile{}, attrs) refute changeset.valid? - assert errors_on(changeset)[:temp_k] != nil + assert errors_on(changeset)[:temp_k] end test "enforces (lat, lon, valid_time) uniqueness" do diff --git a/test/microwaveprop/weather/ncei_metar_client_test.exs b/test/microwaveprop/weather/ncei_metar_client_test.exs index 970d6765..645f4239 100644 --- a/test/microwaveprop/weather/ncei_metar_client_test.exs +++ b/test/microwaveprop/weather/ncei_metar_client_test.exs @@ -38,7 +38,9 @@ defmodule Microwaveprop.Weather.NceiMetarClientTest do end test "handles negative temperatures via M prefix" do - line = "03927KDFW DFW20260115120010103/15/26 12:00:31 5-MIN KDFW 151200Z 36005KT 10SM CLR M02/M08 A3032 560 47 1400 360/05 RMK AO2" + line = + "03927KDFW DFW20260115120010103/15/26 12:00:31 5-MIN KDFW 151200Z 36005KT 10SM CLR M02/M08 A3032 560 47 1400 360/05 RMK AO2" + [obs] = NceiMetarClient.parse(line) # M02 = -2°C = 28.4°F, M08 = -8°C = 17.6°F diff --git a/test/microwaveprop/weather/nexrad_client_test.exs b/test/microwaveprop/weather/nexrad_client_test.exs new file mode 100644 index 00000000..5e086413 --- /dev/null +++ b/test/microwaveprop/weather/nexrad_client_test.exs @@ -0,0 +1,158 @@ +defmodule Microwaveprop.Weather.NexradClientTest do + use ExUnit.Case, async: true + + alias Microwaveprop.Weather.NexradClient + + describe "round_to_5min/1" do + test "rounds down to nearest 5-minute boundary" do + dt = ~U[2022-08-20 14:07:32Z] + assert NexradClient.round_to_5min(dt) == ~U[2022-08-20 14:05:00Z] + end + + test "already on a 5-minute boundary stays unchanged" do + dt = ~U[2022-08-20 14:10:00Z] + assert NexradClient.round_to_5min(dt) == ~U[2022-08-20 14:10:00Z] + end + + test "59 minutes rounds down to 55" do + dt = ~U[2022-08-20 14:59:59Z] + assert NexradClient.round_to_5min(dt) == ~U[2022-08-20 14:55:00Z] + end + end + + describe "frame_url/1" do + test "builds correct IEM n0q URL from datetime" do + dt = ~U[2022-08-20 14:05:00Z] + url = NexradClient.frame_url(dt) + + assert url == + "https://mesonet.agron.iastate.edu/archive/data/2022/08/20/GIS/uscomp/n0q_202208201405.png" + end + + test "pads single-digit month and day" do + dt = ~U[2022-01-05 03:00:00Z] + url = NexradClient.frame_url(dt) + + assert url == + "https://mesonet.agron.iastate.edu/archive/data/2022/01/05/GIS/uscomp/n0q_202201050300.png" + end + end + + describe "latlon_to_pixel/2" do + test "converts coordinates to pixel positions" do + # Top-left corner: lat=50, lon=-126 + {x, y} = NexradClient.latlon_to_pixel(50.0, -126.0) + assert x == 0 + assert y == 0 + end + + test "bottom-right corner" do + # lon = -65 -> x = (-65 - -126) / 0.005 = 61 / 0.005 = 12200 + # lat = 23 -> y = (50 - 23) / 0.005 = 27 / 0.005 = 5400 + {x, y} = NexradClient.latlon_to_pixel(23.0, -65.0) + assert x == 12_200 + assert y == 5400 + end + + test "DFW area" do + # DFW: lat ~32.9, lon ~-97.0 + {x, y} = NexradClient.latlon_to_pixel(32.9, -97.0) + + # x = (-97 - -126) / 0.005 = 29 / 0.005 = 5800 + assert x == 5800 + # y = (50 - 32.9) / 0.005 = 17.1 / 0.005 = 3420 + assert y == 3420 + end + end + + describe "extract_box/4" do + test "extracts stats from a box of pixel values" do + # Create a simple 100x100 image (all zeros except a small region) + width = 100 + height = 100 + pixels = :binary.copy(<<0>>, width * height) + + # Put some non-zero values around pixel (50, 50) - a 5x5 block + pixels = + Enum.reduce(48..52, pixels, fn y, acc -> + Enum.reduce(48..52, acc, fn x, inner_acc -> + offset = y * width + x + <> = inner_acc + <> + end) + end) + + result = NexradClient.extract_box(pixels, width, 50, 50, 5) + + assert result.pixel_count == 25 + assert result.max_reflectivity_dbz > 0 + assert result.mean_reflectivity_dbz > 0 + assert result.texture_variance == 0.0 + end + + test "returns zeros for an empty box" do + width = 100 + height = 100 + pixels = :binary.copy(<<0>>, width * height) + + result = NexradClient.extract_box(pixels, width, 50, 50, 5) + + assert result.pixel_count == 0 + assert result.mean_reflectivity_dbz == 0.0 + assert result.max_reflectivity_dbz == 0.0 + assert result.texture_variance == 0.0 + end + + test "computes texture variance for mixed values" do + width = 100 + height = 100 + pixels = :binary.copy(<<0>>, width * height) + + # Place alternating values: half 100, half 200 + pixels = + Enum.reduce(48..52, pixels, fn y, acc -> + Enum.reduce(48..52, acc, fn x, inner_acc -> + offset = y * width + x + value = if rem(x + y, 2) == 0, do: 100, else: 200 + <> = inner_acc + <> + end) + end) + + result = NexradClient.extract_box(pixels, width, 50, 50, 5) + + assert result.pixel_count == 25 + # Variance should be non-zero since we have mixed values + assert result.texture_variance > 0.0 + end + + test "handles box at image boundary by clamping" do + width = 100 + height = 100 + pixels = :binary.copy(<<50>>, width * height) + + # Corner case: near edge + result = NexradClient.extract_box(pixels, width, 2, 2, 25) + + # Should not crash; pixel_count depends on how many pixels are in bounds + assert result.pixel_count > 0 + end + end + + describe "pixel_to_dbz/1" do + test "value 0 maps to 0 (no echo)" do + assert NexradClient.pixel_to_dbz(0) == 0.0 + end + + test "value 255 maps to approximately 95 dBZ" do + dbz = NexradClient.pixel_to_dbz(255) + assert_in_delta dbz, 95.0, 1.0 + end + + test "mid-range value produces reasonable dBZ" do + dbz = NexradClient.pixel_to_dbz(128) + assert dbz > 0.0 + assert dbz < 95.0 + end + end +end diff --git a/test/microwaveprop/weather/nexrad_observation_test.exs b/test/microwaveprop/weather/nexrad_observation_test.exs new file mode 100644 index 00000000..41417bda --- /dev/null +++ b/test/microwaveprop/weather/nexrad_observation_test.exs @@ -0,0 +1,63 @@ +defmodule Microwaveprop.Weather.NexradObservationTest do + use Microwaveprop.DataCase, async: true + + alias Microwaveprop.Weather.NexradObservation + + @valid_attrs %{ + observed_at: ~U[2022-08-20 20:00:00Z], + lat: 32.875, + lon: -97.0, + mean_reflectivity_dbz: 12.5, + max_reflectivity_dbz: 35.0, + texture_variance: 48.7, + pixel_count: 180 + } + + describe "changeset/2" do + test "valid attributes produce a valid changeset" do + changeset = NexradObservation.changeset(%NexradObservation{}, @valid_attrs) + assert changeset.valid? + end + + test "requires observed_at, lat, and lon" do + changeset = NexradObservation.changeset(%NexradObservation{}, %{}) + + assert %{ + observed_at: ["can't be blank"], + lat: ["can't be blank"], + lon: ["can't be blank"] + } = errors_on(changeset) + end + + test "numeric fields are optional" do + attrs = Map.take(@valid_attrs, [:observed_at, :lat, :lon]) + changeset = NexradObservation.changeset(%NexradObservation{}, attrs) + assert changeset.valid? + end + + test "persists to the database" do + changeset = NexradObservation.changeset(%NexradObservation{}, @valid_attrs) + assert {:ok, obs} = Repo.insert(changeset) + assert obs.lat == 32.875 + assert obs.lon == -97.0 + assert obs.mean_reflectivity_dbz == 12.5 + assert obs.max_reflectivity_dbz == 35.0 + assert obs.texture_variance == 48.7 + assert obs.pixel_count == 180 + end + + test "enforces unique constraint on lat, lon, observed_at" do + assert {:ok, _} = + %NexradObservation{} + |> NexradObservation.changeset(@valid_attrs) + |> Repo.insert() + + assert {:error, changeset} = + %NexradObservation{} + |> NexradObservation.changeset(@valid_attrs) + |> Repo.insert() + + assert %{lat: ["has already been taken"]} = errors_on(changeset) + end + end +end