1846 lines
83 KiB
Rust
1846 lines
83 KiB
Rust
//! Deterministic equilibrium hydrology — lake fill, overflow re-routing, and
|
||
//! gorge carving (T-1177, body-map-viewer workshop measurement ①).
|
||
//!
|
||
//! **What this is.** A pure function of `(seed, heightmap, moisture/climate)`
|
||
//! that computes a *settled end-state* for standing water: every depression
|
||
//! (basin) on the elevation field either (a) fills to a spill level and
|
||
//! overflows onward toward the sea — possibly chaining through further
|
||
//! basins — or (b) is classified **endorheic** (holds its own closed water
|
||
//! cycle, no outflow) when the body's moisture/climate makes evaporative
|
||
//! balance plausible. Where the cheapest overflow path must cross a ridge
|
||
//! higher than the lake's own spill level, the solver carves a deterministic
|
||
//! channel (gorge) — geological backstory computed at generation time, not a
|
||
//! stateful erosion simulation.
|
||
//!
|
||
//! **This is NOT a simulation.** There is no tick loop, no iterative erosion,
|
||
//! no rainfall-accumulation-over-time model. One elevation field in, one
|
||
//! settled hydrology result out, same every time for the same inputs
|
||
//! (D-010). See the "what this prototype does NOT do" list in
|
||
//! `docs/workshops/body-map-viewer/measurements/t1177-hydrology.md` for the
|
||
//! full list of deliberately deferred behaviour.
|
||
//!
|
||
//! **Algorithm (priority-flood family, Barnes/Planchon-Darboux class):**
|
||
//! 1. Priority-flood fill from the boundary/sea inward: a min-heap seeded at
|
||
//! every below-`sea_level` cell floods the grid, raising each interior
|
||
//! cell to the highest of (its own elevation, the water level that reached
|
||
//! it) — the standard "flood from the outside in" formulation. This
|
||
//! produces, for every cell, its *filled* elevation (drainage.rs's
|
||
//! existing `depression_fill` is an iterative relaxation of the same
|
||
//! problem but does not expose spill points or basin membership — this
|
||
//! module needs both, so it re-derives the fill with a heap that *does*
|
||
//! track them, rather than layering basin-detection on top of the
|
||
//! existing black-box result).
|
||
//! 2. Every interior cell where `filled > original` is a *lake cell*; the
|
||
//! heap visit order also yields each basin's spill point (the last
|
||
//! boundary cell that let water into it) essentially for free.
|
||
//! 3. Basins are grouped by contiguous lake-cell connectivity (flood-fill).
|
||
//! Each basin's spill level = the filled elevation at its spill point.
|
||
//! 4. Endorheic classification (per-basin, moisture-governed — see
|
||
//! `is_endorheic`): large + dry basins are declared closed-cycle sinks
|
||
//! (TERMINAL). Everything else overflows: a cheapest-path search
|
||
//! (Dijkstra/A*-family over elevation-crossing cost, same `BinaryHeap`
|
||
//! pattern as `road_graph.rs`'s `astar`) from the spill point to the
|
||
//! nearest lower-or-equal terrain (another basin's interior or the sea)
|
||
//! finds the overflow channel. Cells on that channel above the spill
|
||
//! level are carved down to it — the gorge.
|
||
//! 5. Chaining: an overflowing basin's outlet may land inside *another*
|
||
//! basin's footprint; that basin is processed in the same pass (basins
|
||
//! are handled in ascending spill-level order, so a lower basin is always
|
||
//! resolved before anything can overflow into it a second time).
|
||
//!
|
||
//! Determinism (D-010): elevation and ordering comparisons are on the same
|
||
//! `i64`-scaled integers `drainage.rs` uses; heap tie-breaks are `(cost, cell
|
||
//! index)` so equal-cost frontier cells always resolve in the same order;
|
||
//! every collection that participates in output ordering is a `Vec` keyed by
|
||
//! cell index or basin id (row-major/discovery order), never a
|
||
//! `HashMap`/`HashSet` iteration. **One float comparison exists** —
|
||
//! `is_endorheic`'s `area_frac >= ENDORHEIC_AREA_FLOOR` gate (`f64`, derived
|
||
//! from `basin_cells[b].len() as f64 / n as f64`) — but it is
|
||
//! deterministic-by-derivation: both operands are pure functions of the same
|
||
//! inputs (basin cell count and grid size), computed identically every run,
|
||
//! so the comparison itself always yields the same boolean for the same
|
||
//! input grid. It is not IEEE-754-hazardous in the way accumulated
|
||
//! floating-point summation across runs/platforms can be; it is a single,
|
||
//! reproducible division and comparison.
|
||
|
||
use std::cmp::Reverse;
|
||
use std::collections::{BinaryHeap, VecDeque};
|
||
|
||
use crate::atlas::drainage::d8_offset;
|
||
|
||
/// Scale factor for converting f32 elevation to integer for deterministic
|
||
/// comparison — matches `drainage::ELEV_SCALE` (kept as a private duplicate
|
||
/// rather than exported from `drainage.rs`, since the two modules' integer
|
||
/// domains are independent: this module's `i64` values are always genuine
|
||
/// elevations in the same units, never `drainage.rs`'s post-relaxation
|
||
/// `depression_fill` values).
|
||
const ELEV_SCALE: f64 = 1_000_000.0;
|
||
|
||
/// A cell is a river cell when its flow accumulation exceeds this threshold —
|
||
/// mirrors `drainage::RIVER_THRESHOLD`; only overflow channels that also
|
||
/// qualify as "a real waterway" (i.e. would show up on the river network)
|
||
/// are carved as gorges. A trickling overflow with negligible flow doesn't
|
||
/// get a canyon.
|
||
const RIVER_THRESHOLD: i32 = 200;
|
||
|
||
/// D8 neighbor offsets, same fixed priority order as `drainage::D8` — needed
|
||
/// locally so `d8_offset` calls stay index-compatible with the shared table
|
||
/// without re-exporting the private array.
|
||
const D8_LEN: u8 = 8;
|
||
|
||
/// One basin's settled hydrology outcome.
|
||
#[derive(Debug, Clone, PartialEq, Eq)]
|
||
pub enum BasinOutcome {
|
||
/// The basin overflows: `spill_cell` is where water leaves the basin,
|
||
/// `outlet_path` is the carved channel (spill cell → downstream target,
|
||
/// inclusive of both ends), and `downstream_target` is either another
|
||
/// basin (chained) or the sea.
|
||
Overflow {
|
||
outlet_path: Vec<usize>,
|
||
downstream_target: DownstreamTarget,
|
||
},
|
||
/// The basin holds a closed water cycle — no outflow (TERMINAL).
|
||
Endorheic { reason: EndorheicReason },
|
||
}
|
||
|
||
/// What an overflowing basin's outlet channel terminates at.
|
||
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
|
||
pub enum DownstreamTarget {
|
||
/// Reaches a sub-`sea_level` cell.
|
||
Sea,
|
||
/// Flows into another basin's footprint (chained overflow), by basin id.
|
||
Basin(u32),
|
||
/// **Success case.** Reaches open low ground — at or below the
|
||
/// originating basin's own spill level, but not sea and not another
|
||
/// basin's footprint (the water simply spreads here without needing a
|
||
/// lake label of its own). `outlet_path` on the `BasinOutcome::Overflow`
|
||
/// that carries this variant is a real, complete carved path, same as
|
||
/// the `Sea`/`Basin` cases — this is a genuine overflow terminus, not a
|
||
/// failure. Named for the geomorphological term (an open, uncontained
|
||
/// spillway channel), not to be confused with `EdgeUnreachable` below.
|
||
OpenSpillway,
|
||
/// **Failure case.** The search exhausted its budget (`w * h` node pops)
|
||
/// without reaching ANY valid terminus (`Sea`, `Basin`, or
|
||
/// `OpenSpillway`) — a pathological body that genuinely has no
|
||
/// reachable lower ground within the grid. `outlet_path` on the
|
||
/// `BasinOutcome::Overflow` that carries this variant is the
|
||
/// best-effort partial result (see `cheapest_overflow_path`'s fallback),
|
||
/// not a complete channel. The lake-fill analogue of
|
||
/// `RIVER_DOWNSTREAM_EDGE_DRAIN` in spirit (a grid-topology dead end),
|
||
/// but unlike that sentinel this one is a genuine non-terminus, not a
|
||
/// polar-row artifact.
|
||
EdgeUnreachable,
|
||
}
|
||
|
||
/// Why a basin was classified endorheic — carried for diagnostics/tuning,
|
||
/// not currently read by any downstream consumer.
|
||
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
|
||
pub enum EndorheicReason {
|
||
/// Basin area + dryness cleared the moisture-governed threshold (see
|
||
/// `is_endorheic`).
|
||
MoistureGoverned,
|
||
}
|
||
|
||
/// One filled basin (lake) in the settled hydrology result.
|
||
#[derive(Debug, Clone)]
|
||
pub struct Basin {
|
||
pub basin_id: u32,
|
||
/// Row-major cell indices belonging to this basin's lake footprint
|
||
/// (cells where `filled_elevation > original_elevation`), sorted
|
||
/// ascending for determinism.
|
||
pub cells: Vec<usize>,
|
||
/// The basin's spill level, in the same i64-scaled units as
|
||
/// `drainage::analyze`'s internal representation (divide by
|
||
/// `ELEV_SCALE` to recover the [0,1] float).
|
||
pub spill_level_scaled: i64,
|
||
/// The cell that is the basin's lowest boundary point — the pour point
|
||
/// water first breaches when the basin overtops.
|
||
pub spill_cell: usize,
|
||
pub outcome: BasinOutcome,
|
||
}
|
||
|
||
/// Full result of the equilibrium hydrology solve for one body.
|
||
#[derive(Debug, Clone)]
|
||
pub struct HydrologyResult {
|
||
pub basins: Vec<Basin>,
|
||
/// Per-cell filled elevation (i64-scaled) — the settled water surface
|
||
/// where a lake exists, or the original terrain elevation otherwise.
|
||
/// Row-major, `w × h`.
|
||
pub filled_scaled: Vec<i64>,
|
||
/// Per-cell channel depth in the same i64-scaled elevation units: 0
|
||
/// everywhere except carved gorge cells, where it is
|
||
/// `original_scaled[cell] - carved_floor_scaled[cell]` — see the module
|
||
/// docs' cliff-representation section. Row-major, `w × h`.
|
||
pub channel_depth_scaled: Vec<i64>,
|
||
/// True for every cell that is part of a carved gorge channel (a cliff
|
||
/// edge — the rim/floor discontinuity the workshop's red flag 4 asks
|
||
/// about). Row-major, `w × h`.
|
||
pub cliff_edge: Vec<bool>,
|
||
/// T-1188 lake-margin normalization: for every lake cell, the MAXIMUM
|
||
/// settled depth (`filled_scaled[c] - original_scaled[c]`, i64-scaled)
|
||
/// anywhere in that cell's basin — broadcast to every cell in the basin
|
||
/// (not just the deepest one), `0` for non-lake cells. Row-major, `w ×
|
||
/// h`, computed once here (basin membership — `Basin.cells` — and both
|
||
/// elevation fields are already in scope at the end of `solve()`) rather
|
||
/// than re-derived by every caller.
|
||
///
|
||
/// **Why per-basin, not a single body-wide constant:** the PR #206
|
||
/// eyeball found `district_profile`'s original absolute-depth
|
||
/// `LAKE_MARGIN_DEPTH_CEILING` (a single fixed fraction calibrated
|
||
/// against one body's p90) compressed the bulk of any more-skewed
|
||
/// depth distribution into single-digit `lake_margin_q` values — and,
|
||
/// separately, that heightmap resolution (~40–78 km/px) means
|
||
/// within-basin absolute-depth variation is sometimes sub-texel-tiny
|
||
/// (a real basin can be "shallow" in absolute normalized-elevation
|
||
/// terms while still having a perfectly good internal shore→deep
|
||
/// gradient relative to ITS OWN range). Normalizing against each
|
||
/// basin's own max depth makes every non-degenerate lake use the full
|
||
/// 0–100 tone range regardless of the body's absolute elevation scale —
|
||
/// the perceptual point of the ticket (a lake's edge should read
|
||
/// lighter than its centre, on every lake, not just deep ones).
|
||
pub basin_max_depth_scaled: Vec<i64>,
|
||
}
|
||
|
||
/// Convert an `i64`-scaled elevation value (as carried on
|
||
/// [`HydrologyResult::filled_scaled`] and friends) back to the `[0.0, 1.0]`
|
||
/// normalized fraction the raw heightmap and `sea_level` are expressed in —
|
||
/// the exact inverse of the `(e as f64 * ELEV_SCALE) as i64` conversion
|
||
/// `solve()` applies at its own entry point. `pub(crate)` so callers outside
|
||
/// this module (T-1184: [`crate::atlas::features::TerrainAnalysis::with_hydrology`])
|
||
/// never need to know or duplicate the scale constant — the module that owns
|
||
/// the scaling owns the inverse too.
|
||
pub(crate) fn scaled_to_fraction(scaled: i64) -> f32 {
|
||
(scaled as f64 / ELEV_SCALE) as f32
|
||
}
|
||
|
||
/// Moisture/climate inputs governing the endorheic-vs-overflow decision.
|
||
/// Deliberately minimal and explicitly tunable — see the module docs and the
|
||
/// results doc's "endorheic criterion" section for the rationale and the
|
||
/// knobs future tuning is expected to touch.
|
||
#[derive(Debug, Clone, Copy)]
|
||
pub struct ClimateInputs {
|
||
/// Body-level moisture ceiling, 0–100 (mirrors
|
||
/// `district_profile::derive_moisture_q`'s hydrosphere-derived `ceiling`
|
||
/// before the per-district spatial gradient — a single body-wide value
|
||
/// is enough for a first-pass basin classifier; a per-basin-centroid
|
||
/// value is a documented future refinement, not built here).
|
||
pub moisture_q: i32,
|
||
}
|
||
|
||
/// Solve equilibrium hydrology for one elevation grid. `elevation` is
|
||
/// row-major, shape `height × width`, values in `[0.0, 1.0]`. `sea_level` is
|
||
/// the fraction below which terrain is ocean. `climate` governs the
|
||
/// endorheic-vs-overflow decision (see `is_endorheic`).
|
||
///
|
||
/// Determinism (D-010): pure function of its arguments, no RNG, no
|
||
/// wall-clock, no `HashMap` iteration — same inputs always produce a
|
||
/// byte-identical `HydrologyResult` (verified by the `determinism` test).
|
||
pub fn solve(
|
||
elevation: &[f32],
|
||
width: u32,
|
||
height: u32,
|
||
sea_level: f32,
|
||
climate: ClimateInputs,
|
||
) -> HydrologyResult {
|
||
let w = width as usize;
|
||
let h = height as usize;
|
||
let n = w * h;
|
||
|
||
let original: Vec<i64> = elevation
|
||
.iter()
|
||
.map(|&e| (e as f64 * ELEV_SCALE) as i64)
|
||
.collect();
|
||
let sea_scaled = (sea_level as f64 * ELEV_SCALE) as i64;
|
||
|
||
// ---- Step 1: priority-flood fill from the sea/boundary inward -------
|
||
let filled = priority_flood_fill(&original, w, h, sea_scaled);
|
||
|
||
// ---- Step 2: lake cells + basin grouping -----------------------------
|
||
let is_lake: Vec<bool> = (0..n).map(|i| filled[i] > original[i]).collect();
|
||
let (basin_of, basin_count) = label_lake_basins(&is_lake, w, h);
|
||
|
||
// ---- Step 3: per-basin spill point + spill level ---------------------
|
||
// The spill point is the lake cell adjacent to the LOWEST-original-
|
||
// elevation rim cell outside the basin — i.e. the true topographic
|
||
// saddle water would first overtop, not merely any lake cell touching
|
||
// non-basin terrain. A uniformly-filled lake (the common case: every
|
||
// interior cell floods to the same spill level) makes ALL boundary
|
||
// cells share the same `filled` value, so ranking candidates by
|
||
// `filled` (the former approach) degenerates to an arbitrary
|
||
// lowest-index tie-break — the actual rim geometry never enters the
|
||
// comparison. Ranking by the OUTSIDE neighbor's `original` (unfilled)
|
||
// elevation instead finds the genuine lowest point on the rim, which is
|
||
// what `cheapest_overflow_path`'s Dijkstra search needs as a useful
|
||
// starting point (starting from an arbitrary flat-lake cell instead of
|
||
// the real saddle can send the search off in an unrelated direction —
|
||
// caught by the `wide_catchment_forces_a_carved_gorge` regression test).
|
||
// The basin's spill LEVEL remains `filled[spill_cell]` (the lake's own
|
||
// uniform surface, i.e. the water level at overtopping) — only which
|
||
// cell is chosen as the search's origin changes.
|
||
let mut spill_cell: Vec<Option<usize>> = vec![None; basin_count];
|
||
let mut spill_rim_elev: Vec<i64> = vec![i64::MAX; basin_count];
|
||
let mut basin_cells: Vec<Vec<usize>> = vec![Vec::new(); basin_count];
|
||
for (idx, entry) in basin_of.iter().enumerate().take(n) {
|
||
if let Some(b) = entry {
|
||
basin_cells[*b as usize].push(idx);
|
||
}
|
||
}
|
||
for idx in 0..n {
|
||
let Some(b) = basin_of[idx] else { continue };
|
||
let r = (idx / w) as i32;
|
||
let c = (idx % w) as i32;
|
||
for k in 0..D8_LEN {
|
||
let (dr, dc) = d8_offset(k);
|
||
let nr = r + dr;
|
||
if nr < 0 || nr >= h as i32 {
|
||
continue;
|
||
}
|
||
let nc = (c + dc).rem_euclid(w as i32) as usize;
|
||
let ni = nr as usize * w + nc;
|
||
let neighbor_in_basin = basin_of[ni] == Some(b);
|
||
if neighbor_in_basin {
|
||
continue;
|
||
}
|
||
// Rank by the OUTSIDE neighbor's true (unfilled) elevation —
|
||
// the actual rim height at this crossing point.
|
||
let candidate_rim = original[ni];
|
||
if candidate_rim < spill_rim_elev[b as usize]
|
||
|| (candidate_rim == spill_rim_elev[b as usize]
|
||
&& Some(idx) < spill_cell[b as usize])
|
||
{
|
||
spill_rim_elev[b as usize] = candidate_rim;
|
||
spill_cell[b as usize] = Some(idx);
|
||
}
|
||
}
|
||
}
|
||
// Spill LEVEL is the lake's own filled surface at the chosen cell (the
|
||
// water level at overtopping) — always well-defined once `spill_cell`
|
||
// is set, since every basin cell shares the same `filled` value in the
|
||
// uniform-lake case, and even in non-uniform edge cases `filled[cell]`
|
||
// is exactly the water level AT that cell, which is what carving
|
||
// compares against.
|
||
let spill_level: Vec<i64> = (0..basin_count)
|
||
.map(|b| spill_cell[b].map(|c| filled[c]).unwrap_or(i64::MAX))
|
||
.collect();
|
||
|
||
// ---- Step 4: flow accumulation (reused for the gorge river-threshold
|
||
// gate and for river-network compatibility) ---------------------------
|
||
let fdir = flow_direction_on(&filled, w, h);
|
||
let accum = flow_accumulation_on(&fdir, w, h);
|
||
|
||
// ---- Step 5: per-basin outcome, ascending spill-level order ----------
|
||
// Ascending order guarantees a lower basin is always resolved before a
|
||
// higher one's overflow could chain into it a second time within this
|
||
// pass (a basin's own outcome never depends on a *higher* basin).
|
||
let mut order: Vec<usize> = (0..basin_count).collect();
|
||
order.sort_by_key(|&b| (spill_level[b], b));
|
||
|
||
let mut channel_depth_scaled = vec![0i64; n];
|
||
let mut cliff_edge = vec![false; n];
|
||
let mut outcomes: Vec<Option<BasinOutcome>> = vec![None; basin_count];
|
||
|
||
for b in order {
|
||
let Some(spill) = spill_cell[b] else {
|
||
// No boundary found (whole-grid lake edge case) — treat as
|
||
// unreachable overflow.
|
||
outcomes[b] = Some(BasinOutcome::Overflow {
|
||
outlet_path: vec![],
|
||
downstream_target: DownstreamTarget::EdgeUnreachable,
|
||
});
|
||
continue;
|
||
};
|
||
let area_frac = basin_cells[b].len() as f64 / n as f64;
|
||
if is_endorheic(area_frac, climate.moisture_q) {
|
||
outcomes[b] = Some(BasinOutcome::Endorheic {
|
||
reason: EndorheicReason::MoistureGoverned,
|
||
});
|
||
continue;
|
||
}
|
||
|
||
// Cheapest overflow path: Dijkstra from the spill cell to the
|
||
// nearest cell whose ORIGINAL elevation is <= the basin's spill
|
||
// level (i.e. water reaching it needs no further carving) OR that
|
||
// belongs to a different, already-resolved (lower) basin, OR the
|
||
// sea. Cost = elevation the path must climb above the spill level
|
||
// (carving cost), so the search naturally prefers the topographic
|
||
// saddle.
|
||
let (path, target) = cheapest_overflow_path(
|
||
spill,
|
||
spill_level[b],
|
||
&original,
|
||
&basin_of,
|
||
b as u32,
|
||
sea_scaled,
|
||
w,
|
||
h,
|
||
);
|
||
|
||
// Carve: every cell on the path whose original elevation exceeds
|
||
// the spill level is cut down to it. Only carve where the channel
|
||
// also clears the river threshold (a real waterway) — a spill path
|
||
// through nearly flat ground with negligible flow doesn't need a
|
||
// canyon, it just floods.
|
||
let accum_gate = path
|
||
.iter()
|
||
.map(|&c| accum.get(c).copied().unwrap_or(0))
|
||
.max()
|
||
.unwrap_or(0);
|
||
if accum_gate > RIVER_THRESHOLD {
|
||
for &cell in &path {
|
||
if original[cell] > spill_level[b] {
|
||
channel_depth_scaled[cell] = original[cell] - spill_level[b];
|
||
cliff_edge[cell] = true;
|
||
}
|
||
}
|
||
}
|
||
|
||
outcomes[b] = Some(BasinOutcome::Overflow {
|
||
outlet_path: path,
|
||
downstream_target: target,
|
||
});
|
||
}
|
||
|
||
let basins: Vec<Basin> = (0..basin_count)
|
||
.map(|b| Basin {
|
||
basin_id: b as u32,
|
||
cells: basin_cells[b].clone(),
|
||
spill_level_scaled: spill_level[b],
|
||
spill_cell: spill_cell[b].unwrap_or(0),
|
||
outcome: outcomes[b].clone().unwrap_or(BasinOutcome::Endorheic {
|
||
reason: EndorheicReason::MoistureGoverned,
|
||
}),
|
||
})
|
||
.collect();
|
||
|
||
// T-1188: per-basin max depth, broadcast to every cell in that basin.
|
||
// `basin_cells[b]` is exactly the membership list `Basin.cells` above is
|
||
// built from — walked directly here rather than through `basins` to
|
||
// avoid a second clone of every basin's cell list.
|
||
let mut basin_max_depth_scaled = vec![0i64; n];
|
||
for cells in basin_cells.iter().take(basin_count) {
|
||
let max_depth = cells
|
||
.iter()
|
||
.map(|&c| (filled[c] - original[c]).max(0))
|
||
.max()
|
||
.unwrap_or(0);
|
||
for &c in cells {
|
||
basin_max_depth_scaled[c] = max_depth;
|
||
}
|
||
}
|
||
|
||
HydrologyResult {
|
||
basins,
|
||
filled_scaled: filled,
|
||
channel_depth_scaled,
|
||
cliff_edge,
|
||
basin_max_depth_scaled,
|
||
}
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Basin-outlet → D8 river-network wiring (T-1185, D-227 amendment (4))
|
||
// ---------------------------------------------------------------------------
|
||
|
||
/// Extend a D8-derived [`RiverNetwork`] with **outlet continuations** for
|
||
/// every [`BasinOutcome::Overflow`] basin in `result` — the wiring the D-227
|
||
/// amendment (4) endorheic cue rests on ("outflow-course PRESENCE, no wire
|
||
/// bit, no 18th zone"). Additive per T-1170 Ruling 7b: `Endorheic` basins
|
||
/// contribute nothing (their absence of a continuation IS the cue, already
|
||
/// read correctly by every existing consumer — [`crate::atlas::river_course::
|
||
/// build_edges`] already no-ops on cells it never sees), and every appended
|
||
/// cell is a plain new [`RiverNetwork::river_cells`] entry with a real D8
|
||
/// [`RiverNetwork::river_downstream`] pointer, indistinguishable downstream
|
||
/// from a naturally-extracted river cell — no new wire field, no new sentinel.
|
||
///
|
||
/// **Edge identity (D-010 determinism).** Each new river cell's identity is
|
||
/// its own `(row, col)` position (the same convention every existing river
|
||
/// cell already uses — [`crate::atlas::river_course::pack_cell_id`] derives
|
||
/// `edge_id` from the upstream cell's packed position) — never an insertion
|
||
/// index or iteration order. Basins are walked in ascending `basin_id` order
|
||
/// (itself a deterministic row-major discovery order, [`label_lake_basins`]),
|
||
/// and `outlet_path` is walked in its own fixed emission order (a Dijkstra
|
||
/// `came`-chain reconstruction, deterministic per basin) — so two independent
|
||
/// solves of the same input produce byte-identical appended cells in
|
||
/// byte-identical order. No hashing, no `HashMap`, nothing keyed on wall-clock
|
||
/// or thread scheduling.
|
||
///
|
||
/// **Cell selection — the spill cell ALWAYS gets a real outflow pointer
|
||
/// (PR #202 review, Hoshe finding 1).** `outlet_path` is "inclusive of both
|
||
/// ends": index 0 is the basin's `spill_cell` (a `Lake` gridunit at the
|
||
/// basin's own boundary — morphology sourcing and river-cell membership are
|
||
/// independent, T-1184, so a cell can legitimately be both), index
|
||
/// `len - 1` is the downstream terminus. The spill cell is the ONE cell that
|
||
/// visually anchors the lake's exit: every course drawn on the map starts at
|
||
/// a `RiverEdge::upstream` position ([`crate::atlas::river_course::
|
||
/// invent_course`]'s `anchor_a`), so unless the spill cell itself carries an
|
||
/// outflow-pointing entry, NOTHING in the served network visually touches
|
||
/// this lake — a real `Overflow` basin would read exactly like `Endorheic`
|
||
/// (no course leaving the lake), which is the CUE this ticket exists to
|
||
/// deliver, not an incidental gap. This is why the spill cell is never
|
||
/// skipped, unconditionally: it always gets an entry pointing along
|
||
/// `outlet_path` (D8 direction toward `outlet_path[1]`, or straight to the
|
||
/// terminus sentinel on a length-1 residual path — see below), whether or
|
||
/// not it was already a `river_cells` member. If it WAS already present (a
|
||
/// real, if rare, case — verified on GJ1c: 2/51 real `Overflow` basins have
|
||
/// a pre-existing-river-cell spill point), its EXISTING entry is
|
||
/// overwritten in place (same array index, same position, new downstream
|
||
/// direction/class/seaward) rather than appended — appending would create a
|
||
/// second `river_cells` entry at the same `(row, col)`, which
|
||
/// [`crate::atlas::river_course::build_edges`] would turn into two
|
||
/// `RiverEdge`s sharing the same `edge_id` (`pack_cell_id` is a pure
|
||
/// function of position, Ruling 2d), corrupting the
|
||
/// one-edge-per-upstream-cell invariant `edge_ids_are_unique` guards.
|
||
/// Overwriting is correct, not merely safe: the pre-existing pointer was
|
||
/// computed by `extract_river_network`'s D8 walk on the ORIGINAL
|
||
/// (unfilled) surface, which has no knowledge of the basin's spill
|
||
/// direction — the hydrology solve is the more authoritative answer for
|
||
/// what water actually does at this specific cell once the basin is full,
|
||
/// so it wins.
|
||
///
|
||
/// Every remaining interior path cell (index 1 through `len - 2`) becomes a
|
||
/// new river cell whose downstream pointer is the real D8 direction toward
|
||
/// the NEXT path cell (mirroring `extract_river_network`'s own `fdir[i]`
|
||
/// convention exactly — `outlet_path` is itself a chain of true D8
|
||
/// neighbors, guaranteed by construction: [`cheapest_overflow_path`]'s
|
||
/// Dijkstra only ever relaxes D8-adjacent cells). The final path cell's
|
||
/// downstream pointer is set from `downstream_target`: [`DownstreamTarget::
|
||
/// Sea`] maps to [`RIVER_DOWNSTREAM_MOUTH`] (a real river reaching the sea
|
||
/// IS a mouth, the same semantics `extract_river_network` already assigns),
|
||
/// and [`DownstreamTarget::Basin`]/[`DownstreamTarget::OpenSpillway`] both
|
||
/// map to [`RIVER_DOWNSTREAM_EDGE_DRAIN`] (the course simply ends at the
|
||
/// last station — there is no meaningful "direction" once the outlet has
|
||
/// reached another lake's footprint or open ground, the same semantics
|
||
/// `build_edges` already gives a grid-artifact edge-drain: no further chord
|
||
/// to invent).
|
||
///
|
||
/// **Interior cells (index ≥ 1) already present in `river_cells` STOP the
|
||
/// walk** (an outlet path can legitimately re-enter the D8-extracted
|
||
/// network, e.g. `DownstreamTarget::Basin` chaining through a stretch of
|
||
/// terrain the original extraction already classified as a river) — this is
|
||
/// safe (unlike the spill-cell case above) precisely BECAUSE the spill cell
|
||
/// always got its own real entry first: the chain from the lake is never
|
||
/// silently dropped, only its LATER re-entry into pre-existing machinery is
|
||
/// deduplicated. The cell immediately before the collision (which may BE the
|
||
/// spill cell itself, on a length-2 path) correctly points its D8 direction
|
||
/// at the existing river cell, which already has its own onward pointer —
|
||
/// the course continues through machinery that was already there.
|
||
///
|
||
/// **`EdgeUnreachable` paths and a length-1 residual (spill-only) path still
|
||
/// get the spill cell's real entry, just with no interior D8 hop.** A path
|
||
/// of length ≤ 1 (`outlet_path == [spill]` — the `EdgeUnreachable`
|
||
/// best-effort case, or any basin whose search terminates immediately) has
|
||
/// no `outlet_path[1]` to point toward, so the spill cell's downstream
|
||
/// sentinel is taken directly from `downstream_target` (the same terminus
|
||
/// mapping the multi-cell path's LAST cell uses) instead of a D8 direction —
|
||
/// still a real, readable outflow entry, never silently dropped.
|
||
///
|
||
/// **`river_class` scope note.** Every appended cell is classified `0`
|
||
/// (stream — [`crate::atlas::river_course`]'s narrowest, most conservative
|
||
/// meander-amplitude multiplier), a fixed default rather than a value
|
||
/// derived from the outlet channel's own flow accumulation. This is a
|
||
/// deliberate scope boundary, not an oversight: this ticket wires PRESENCE
|
||
/// (does a course leave this lake, yes/no — the D-227 amendment (4) cue),
|
||
/// not a believability pass on outlet-channel *classification*. A basin's
|
||
/// outlet plausibly deserves a wider/trunk-scaled meander (it concentrates
|
||
/// the whole basin's catchment through one channel) — that is a real,
|
||
/// separate refinement, deferred rather than guessed at here.
|
||
pub fn extend_river_network_with_basin_outlets(
|
||
result: &HydrologyResult,
|
||
width: u32,
|
||
height: u32,
|
||
mut network: crate::atlas::body_world_state::RiverNetwork,
|
||
) -> crate::atlas::body_world_state::RiverNetwork {
|
||
use crate::atlas::body_world_state::RIVER_DOWNSTREAM_MOUTH;
|
||
|
||
let w = width as usize;
|
||
let n = (width as usize) * (height as usize);
|
||
|
||
// Position -> array-index for the network's pre-existing river cells — a
|
||
// dense `Vec` (same shape/determinism story as `is_lake`/`basin_of`
|
||
// above: O(1) point lookup, never iterated, so there is no D-010
|
||
// HashMap-iteration-order concern to even raise). Built once, outside
|
||
// the basin loop (basins can number in the thousands on a real body, per
|
||
// the T-1177 population survey: 22,270 basins across 267 bodies).
|
||
// Doubles as the membership test the interior-cell dedup walk needs; the
|
||
// INDEX half is what the spill-cell overwrite case (Hoshe finding 1)
|
||
// needs to mutate the correct existing entry in place rather than
|
||
// appending a duplicate.
|
||
let mut river_cell_index: Vec<Option<usize>> = vec![None; n];
|
||
for (idx, &(r, c)) in network.river_cells.iter().enumerate() {
|
||
river_cell_index[r as usize * w + c as usize] = Some(idx);
|
||
}
|
||
|
||
// Basins are already stored in ascending `basin_id` order (the `(0..
|
||
// basin_count).map(...)` construction above) — iterate as-is rather than
|
||
// re-sorting, preserving the deterministic row-major discovery order
|
||
// `label_lake_basins` assigned.
|
||
for basin in &result.basins {
|
||
let BasinOutcome::Overflow {
|
||
outlet_path,
|
||
downstream_target,
|
||
} = &basin.outcome
|
||
else {
|
||
continue; // Endorheic — no continuation; absence IS the cue.
|
||
};
|
||
if outlet_path.is_empty() {
|
||
continue; // Degenerate — no spill cell at all to anchor from.
|
||
}
|
||
|
||
// The spill cell's downstream sentinel (fn doc's "Cell selection"
|
||
// section): real D8 direction toward `outlet_path[1]` when an
|
||
// interior cell exists, otherwise the terminus mapping directly (a
|
||
// length-1 residual/EdgeUnreachable path) — either way, ALWAYS a
|
||
// real entry, never skipped.
|
||
let spill = outlet_path[0];
|
||
let (sr, sc) = (spill / w, spill % w);
|
||
let spill_sentinel = if outlet_path.len() >= 2 {
|
||
let next = outlet_path[1];
|
||
let (nr, nc) = (next / w, next % w);
|
||
d8_direction_between((sr, sc), (nr, nc), w)
|
||
} else {
|
||
terminus_sentinel(downstream_target)
|
||
};
|
||
let spill_seaward = seaward_for(spill_sentinel, sr, sc, RIVER_DOWNSTREAM_MOUTH);
|
||
|
||
match river_cell_index[spill] {
|
||
Some(existing_idx) => {
|
||
// Overwrite in place (fn doc: the hydrology solve is the
|
||
// more authoritative answer for this cell's true downstream
|
||
// direction than the original flat D8 extraction) — never
|
||
// append, which would duplicate `edge_id` at this position.
|
||
network.river_class[existing_idx] = 0;
|
||
network.river_downstream[existing_idx] = spill_sentinel;
|
||
network.river_seaward[existing_idx] = spill_seaward;
|
||
}
|
||
None => {
|
||
let new_idx = network.river_cells.len();
|
||
network.river_cells.push((sr as u16, sc as u16));
|
||
network.river_class.push(0);
|
||
network.river_downstream.push(spill_sentinel);
|
||
network.river_seaward.push(spill_seaward);
|
||
river_cell_index[spill] = Some(new_idx);
|
||
}
|
||
}
|
||
|
||
if outlet_path.len() < 2 {
|
||
continue; // No interior cell beyond the spill entry above.
|
||
}
|
||
|
||
for i in 1..outlet_path.len() {
|
||
let cell = outlet_path[i];
|
||
let (r, c) = (cell / w, cell % w);
|
||
if river_cell_index[cell].is_some() {
|
||
// Re-entered the pre-existing D8 network (e.g. chained into
|
||
// another basin's already-extracted river reach) — the prior
|
||
// cell's downstream pointer (the spill entry above, or a
|
||
// PREVIOUS loop iteration) already points here, so the chain
|
||
// is visually continuous without adding a duplicate entry.
|
||
// Safe here (unlike the spill cell) because the spill entry
|
||
// above already guarantees the lake's own outflow is never
|
||
// silently dropped — this only dedups a LATER re-entry.
|
||
break;
|
||
}
|
||
|
||
// Downstream pointer: real D8 direction toward the NEXT path
|
||
// cell, or the terminus sentinel at the path's own end.
|
||
let sentinel = if i + 1 < outlet_path.len() {
|
||
let next = outlet_path[i + 1];
|
||
let (nr, nc) = (next / w, next % w);
|
||
d8_direction_between((r, c), (nr, nc), w)
|
||
} else {
|
||
terminus_sentinel(downstream_target)
|
||
};
|
||
|
||
let new_idx = network.river_cells.len();
|
||
network.river_cells.push((r as u16, c as u16));
|
||
network.river_class.push(0); // stream — fixed default, see fn doc's scope note
|
||
network.river_downstream.push(sentinel);
|
||
network
|
||
.river_seaward
|
||
.push(seaward_for(sentinel, r, c, RIVER_DOWNSTREAM_MOUTH));
|
||
river_cell_index[cell] = Some(new_idx);
|
||
}
|
||
}
|
||
|
||
network
|
||
}
|
||
|
||
/// Terminus sentinel for a [`BasinOutcome::Overflow`]'s downstream end, from
|
||
/// [`DownstreamTarget`] — shared by both the spill-cell entry (a length-1
|
||
/// residual path with no interior cell to point at) and the last interior
|
||
/// path cell's entry (the common multi-cell case).
|
||
/// [`DownstreamTarget::Sea`] maps to [`RIVER_DOWNSTREAM_MOUTH`] (a real
|
||
/// river reaching the sea IS a mouth, the same semantics
|
||
/// `extract_river_network` already assigns); every other variant
|
||
/// ([`DownstreamTarget::Basin`], [`DownstreamTarget::OpenSpillway`],
|
||
/// [`DownstreamTarget::EdgeUnreachable`]) maps to
|
||
/// [`RIVER_DOWNSTREAM_EDGE_DRAIN`] — no further chord to invent once the
|
||
/// outlet has reached another lake's footprint, open ground, or exhausted
|
||
/// its search budget.
|
||
fn terminus_sentinel(target: &DownstreamTarget) -> u8 {
|
||
use crate::atlas::body_world_state::{RIVER_DOWNSTREAM_EDGE_DRAIN, RIVER_DOWNSTREAM_MOUTH};
|
||
match target {
|
||
DownstreamTarget::Sea => RIVER_DOWNSTREAM_MOUTH,
|
||
DownstreamTarget::Basin(_)
|
||
| DownstreamTarget::OpenSpillway
|
||
| DownstreamTarget::EdgeUnreachable => RIVER_DOWNSTREAM_EDGE_DRAIN,
|
||
}
|
||
}
|
||
|
||
/// [`RiverNetwork::river_seaward`] value for a cell whose downstream
|
||
/// sentinel is `sentinel` — a real seaward neighbor when `sentinel ==
|
||
/// mouth_sentinel` (see [`extend_river_network_with_basin_outlets`]'s doc on
|
||
/// why `(r, c)` itself, not a cell beyond it, is the correct seaward value
|
||
/// here — mirrors `extract_river_network`'s own `(nr, nc)` capture), the
|
||
/// unreadable `(0, 0)` placeholder otherwise (matches every non-MOUTH
|
||
/// `river_seaward` entry across the codebase).
|
||
fn seaward_for(sentinel: u8, r: usize, c: usize, mouth_sentinel: u8) -> (u16, u16) {
|
||
if sentinel == mouth_sentinel {
|
||
(r as u16, c as u16)
|
||
} else {
|
||
(0, 0)
|
||
}
|
||
}
|
||
|
||
/// The D8 direction index `k` (matching [`d8_offset`]'s table) such that
|
||
/// stepping from `(r, c)` by `d8_offset(k)` (with horizontal wraparound,
|
||
/// [`crate::atlas::drainage`]'s convention) reaches `(nr, nc)`. Panics if the
|
||
/// two cells are not true D8 neighbors — a programmer error (every caller
|
||
/// derives `(nr, nc)` from a Dijkstra path that only ever relaxes true D8
|
||
/// neighbors, so this is a self-consistency assertion, not a runtime
|
||
/// condition that can legitimately fail on real hydrology output).
|
||
fn d8_direction_between(from: (usize, usize), to: (usize, usize), w: usize) -> u8 {
|
||
let dr = to.0 as i32 - from.0 as i32;
|
||
// Shortest signed wrap-around column delta in {-1, 0, 1} — the grid wraps
|
||
// horizontally (equirectangular), same as every D8 walk in this module.
|
||
let raw_dc = to.1 as i32 - from.1 as i32;
|
||
let dc = if raw_dc == 0 {
|
||
0
|
||
} else if raw_dc == 1 || raw_dc == -(w as i32 - 1) {
|
||
1
|
||
} else if raw_dc == -1 || raw_dc == (w as i32 - 1) {
|
||
-1
|
||
} else {
|
||
panic!(
|
||
"d8_direction_between: ({from:?}) -> ({to:?}) is not a D8 neighbor step (dc={raw_dc})"
|
||
);
|
||
};
|
||
for k in 0..D8_LEN {
|
||
if d8_offset(k) == (dr, dc) {
|
||
return k;
|
||
}
|
||
}
|
||
panic!(
|
||
"d8_direction_between: ({from:?}) -> ({to:?}) is not a D8 neighbor step (dr={dr}, dc={dc})"
|
||
);
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Priority-flood fill (Barnes/Planchon-Darboux class)
|
||
// ---------------------------------------------------------------------------
|
||
|
||
/// Classic priority-flood depression fill: seed a min-heap with every
|
||
/// boundary/sea cell, then repeatedly pop the lowest-water-level frontier
|
||
/// cell and raise each unvisited neighbor to `max(neighbor_original,
|
||
/// popped_level)`. Produces the same "filled to spill level" surface as
|
||
/// `drainage::depression_fill`'s iterative relaxation, but in one
|
||
/// `O(n log n)` pass instead of up to 10 full-grid passes, and exposes basin
|
||
/// membership as a side effect of visit order (used by `label_lake_basins`
|
||
/// downstream via the `filled > original` comparison — no separate flag
|
||
/// needed).
|
||
///
|
||
/// Determinism: heap ordering is `(level, cell_index)` — `Reverse` for a
|
||
/// min-heap, cell index as the tie-break, so equal-level frontier cells
|
||
/// always pop in the same order regardless of insertion order.
|
||
fn priority_flood_fill(original: &[i64], w: usize, h: usize, sea_scaled: i64) -> Vec<i64> {
|
||
let n = w * h;
|
||
let mut filled = original.to_vec();
|
||
let mut visited = vec![false; n];
|
||
let mut heap: BinaryHeap<Reverse<(i64, usize)>> = BinaryHeap::new();
|
||
|
||
// Seed: every below-sea-level cell AND every grid-edge row cell (the
|
||
// equirectangular grid has no south/north neighbor past the poles, so
|
||
// the top/bottom rows are boundary same as a coastline — this mirrors
|
||
// `drainage.rs`'s row-bounded D8 walk, which never routes off the
|
||
// top/bottom edge either).
|
||
for idx in 0..n {
|
||
let r = idx / w;
|
||
if original[idx] <= sea_scaled || r == 0 || r == h - 1 {
|
||
visited[idx] = true;
|
||
filled[idx] = original[idx];
|
||
heap.push(Reverse((filled[idx], idx)));
|
||
}
|
||
}
|
||
|
||
while let Some(Reverse((level, idx))) = heap.pop() {
|
||
let r = (idx / w) as i32;
|
||
let c = (idx % w) as i32;
|
||
for k in 0..D8_LEN {
|
||
let (dr, dc) = d8_offset(k);
|
||
let nr = r + dr;
|
||
if nr < 0 || nr >= h as i32 {
|
||
continue;
|
||
}
|
||
let nc = (c + dc).rem_euclid(w as i32) as usize;
|
||
let ni = nr as usize * w + nc;
|
||
if visited[ni] {
|
||
continue;
|
||
}
|
||
visited[ni] = true;
|
||
let new_level = original[ni].max(level);
|
||
filled[ni] = new_level;
|
||
heap.push(Reverse((new_level, ni)));
|
||
}
|
||
}
|
||
|
||
filled
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Lake basin labeling
|
||
// ---------------------------------------------------------------------------
|
||
|
||
/// Group connected `is_lake` cells into basins via row-major-seeded BFS
|
||
/// flood fill. Returns `(basin_of, basin_count)` — `basin_of[idx]` is `None`
|
||
/// for non-lake cells. Basin ids are assigned in the order their seed cell
|
||
/// is first reached in row-major scan order, so relabeling is stable across
|
||
/// runs (D-010).
|
||
fn label_lake_basins(is_lake: &[bool], w: usize, h: usize) -> (Vec<Option<u32>>, usize) {
|
||
let n = w * h;
|
||
let mut basin_of: Vec<Option<u32>> = vec![None; n];
|
||
let mut next_id: u32 = 0;
|
||
|
||
for start in 0..n {
|
||
if !is_lake[start] || basin_of[start].is_some() {
|
||
continue;
|
||
}
|
||
let id = next_id;
|
||
next_id += 1;
|
||
let mut queue = VecDeque::new();
|
||
queue.push_back(start);
|
||
basin_of[start] = Some(id);
|
||
while let Some(idx) = queue.pop_front() {
|
||
let r = (idx / w) as i32;
|
||
let c = (idx % w) as i32;
|
||
for k in 0..D8_LEN {
|
||
let (dr, dc) = d8_offset(k);
|
||
let nr = r + dr;
|
||
if nr < 0 || nr >= h as i32 {
|
||
continue;
|
||
}
|
||
let nc = (c + dc).rem_euclid(w as i32) as usize;
|
||
let ni = nr as usize * w + nc;
|
||
if is_lake[ni] && basin_of[ni].is_none() {
|
||
basin_of[ni] = Some(id);
|
||
queue.push_back(ni);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
(basin_of, next_id as usize)
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Endorheic classification
|
||
// ---------------------------------------------------------------------------
|
||
|
||
/// Whether a basin of the given size (as a fraction of the body's total
|
||
/// surface) is declared **endorheic** (closed water cycle, no outflow)
|
||
/// rather than overflowing.
|
||
///
|
||
/// **Tunable — this is a first-pass heuristic, not a calibrated physical
|
||
/// model.** The criterion: a basin is endorheic when it is large enough that
|
||
/// its own evaporative surface plausibly balances its catchment's inflow
|
||
/// (bigger basin = more evaporating surface per unit of catchment,
|
||
/// Great-Salt-Lake/Caspian-style), UNLESS the body is wet enough that
|
||
/// evaporation is unlikely to keep pace (`moisture_q` high) — a wet body's
|
||
/// basins overflow far more often than a dry one's. Concretely:
|
||
///
|
||
/// - `area_frac >= ENDORHEIC_AREA_FLOOR` (the basin is large enough to matter
|
||
/// at all — tiny depressions always overflow, they're deterministically
|
||
/// more likely to be swamped by their catchment).
|
||
/// - AND `moisture_q <= ENDORHEIC_MOISTURE_CEILING` (dry-to-moderate bodies
|
||
/// only; a saturated/oceanic body's basins are assumed to always find an
|
||
/// outflow).
|
||
///
|
||
/// Both constants and the boolean-AND shape are placeholders pending
|
||
/// real-world calibration (see the results doc's "endorheic criterion"
|
||
/// section) — a follow-up could scale the area floor continuously against
|
||
/// moisture instead of two independent gates.
|
||
const ENDORHEIC_AREA_FLOOR: f64 = 0.004; // ~0.4% of body surface
|
||
const ENDORHEIC_MOISTURE_CEILING: i32 = 60;
|
||
|
||
fn is_endorheic(area_frac: f64, moisture_q: i32) -> bool {
|
||
area_frac >= ENDORHEIC_AREA_FLOOR && moisture_q <= ENDORHEIC_MOISTURE_CEILING
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Cheapest overflow path (Dijkstra, carving-cost weighted)
|
||
// ---------------------------------------------------------------------------
|
||
|
||
/// Dijkstra from `spill` outward, cost = cumulative elevation carved above
|
||
/// `spill_level` (0 for any step that stays at or below spill level).
|
||
/// Terminates at the first cell that is: (a) below `sea_scaled`
|
||
/// (`DownstreamTarget::Sea`), (b) inside a different basin
|
||
/// (`DownstreamTarget::Basin(id)`), or (c) has original elevation `<=
|
||
/// spill_level` and is not part of basin `own_basin` — open low ground, a
|
||
/// genuine SUCCESS terminus reported as `DownstreamTarget::OpenSpillway`
|
||
/// with a real, complete `outlet_path` (the water simply spreads here
|
||
/// without needing a lake label of its own; not sea, not another basin, but
|
||
/// still a valid place for the overflow to end).
|
||
///
|
||
/// Search budget: capped at `w * h` node pops (a full-grid worst case), so a
|
||
/// pathological body can't spin forever — returns
|
||
/// `DownstreamTarget::EdgeUnreachable` (the FAILURE case: exhaustion, no
|
||
/// terminus of any kind found) with only the best-effort partial path if
|
||
/// the budget runs out.
|
||
#[allow(clippy::too_many_arguments)]
|
||
fn cheapest_overflow_path(
|
||
spill: usize,
|
||
spill_level: i64,
|
||
original: &[i64],
|
||
basin_of: &[Option<u32>],
|
||
own_basin: u32,
|
||
sea_scaled: i64,
|
||
w: usize,
|
||
h: usize,
|
||
) -> (Vec<usize>, DownstreamTarget) {
|
||
let n = w * h;
|
||
let mut best_cost = vec![i64::MAX; n];
|
||
let mut came: Vec<usize> = vec![usize::MAX; n];
|
||
best_cost[spill] = 0;
|
||
let mut heap: BinaryHeap<Reverse<(i64, usize)>> = BinaryHeap::new();
|
||
heap.push(Reverse((0, spill)));
|
||
|
||
let mut budget = n;
|
||
while let Some(Reverse((cost, idx))) = heap.pop() {
|
||
if cost > best_cost[idx] {
|
||
continue;
|
||
}
|
||
budget = budget.saturating_sub(1);
|
||
if budget == 0 {
|
||
break;
|
||
}
|
||
|
||
// Termination check (skip at the spill cell itself — it's inside
|
||
// its own basin by definition).
|
||
if idx != spill {
|
||
if original[idx] <= sea_scaled {
|
||
return (reconstruct_path(&came, idx), DownstreamTarget::Sea);
|
||
}
|
||
if let Some(other) = basin_of[idx] {
|
||
if other != own_basin {
|
||
return (reconstruct_path(&came, idx), DownstreamTarget::Basin(other));
|
||
}
|
||
} else if original[idx] <= spill_level {
|
||
// Open low ground, not another basin, not sea: a genuine
|
||
// SUCCESS terminus — the water simply spreads here without
|
||
// needing a lake label of its own. Reported as
|
||
// `OpenSpillway`, distinct from `EdgeUnreachable` (that
|
||
// variant is reserved for genuine search exhaustion, below)
|
||
// — this arm always carries a real, complete path.
|
||
return (reconstruct_path(&came, idx), DownstreamTarget::OpenSpillway);
|
||
}
|
||
}
|
||
|
||
let r = (idx / w) as i32;
|
||
let c = (idx % w) as i32;
|
||
for k in 0..D8_LEN {
|
||
let (dr, dc) = d8_offset(k);
|
||
let nr = r + dr;
|
||
if nr < 0 || nr >= h as i32 {
|
||
continue;
|
||
}
|
||
let nc = (c + dc).rem_euclid(w as i32) as usize;
|
||
let ni = nr as usize * w + nc;
|
||
let step_cost = (original[ni] - spill_level).max(0);
|
||
let tentative = cost + step_cost;
|
||
if tentative < best_cost[ni] {
|
||
best_cost[ni] = tentative;
|
||
came[ni] = idx;
|
||
heap.push(Reverse((tentative, ni)));
|
||
}
|
||
}
|
||
}
|
||
|
||
// Exhausted without reaching a terminus — best-effort: no path.
|
||
(vec![spill], DownstreamTarget::EdgeUnreachable)
|
||
}
|
||
|
||
fn reconstruct_path(came: &[usize], goal: usize) -> Vec<usize> {
|
||
let mut path = Vec::new();
|
||
let mut cur = goal;
|
||
loop {
|
||
path.push(cur);
|
||
let prev = came[cur];
|
||
if prev == usize::MAX {
|
||
break;
|
||
}
|
||
cur = prev;
|
||
}
|
||
path.reverse();
|
||
path
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Local D8 flow direction / accumulation on the FILLED surface
|
||
// ---------------------------------------------------------------------------
|
||
//
|
||
// Duplicated (not imported) from `drainage.rs` deliberately: this module's
|
||
// flow analysis runs on the *filled* (post-lake) surface to gate gorge
|
||
// carving by "is this overflow channel a real waterway", which is a
|
||
// different input than `drainage::analyze`'s own internal
|
||
// `depression_fill` — sharing the function would require exposing
|
||
// `drainage.rs` internals that are `fn`-private by design (not `pub(crate)`
|
||
// like `d8_offset`). Small (30 lines) and self-contained; not worth a
|
||
// visibility change to the existing, already-reviewed module for a
|
||
// same-module measurement.
|
||
|
||
fn flow_direction_on(filled: &[i64], w: usize, h: usize) -> Vec<i8> {
|
||
let mut fdir = vec![-1i8; w * h];
|
||
for r in 0..h {
|
||
for c in 0..w {
|
||
let elev = filled[r * w + c];
|
||
let mut best_drop = 0i64;
|
||
let mut best_k: i8 = -1;
|
||
for k in 0..D8_LEN {
|
||
let (dr, dc) = d8_offset(k);
|
||
let nr = r as i32 + dr;
|
||
let nc = (c as i32 + dc).rem_euclid(w as i32) as usize;
|
||
if nr < 0 || nr >= h as i32 {
|
||
continue;
|
||
}
|
||
let drop = elev - filled[nr as usize * w + nc];
|
||
if drop > best_drop {
|
||
best_drop = drop;
|
||
best_k = k as i8;
|
||
}
|
||
}
|
||
fdir[r * w + c] = best_k;
|
||
}
|
||
}
|
||
fdir
|
||
}
|
||
|
||
fn flow_accumulation_on(fdir: &[i8], w: usize, h: usize) -> Vec<i32> {
|
||
let n = w * h;
|
||
let mut in_degree = vec![0i32; n];
|
||
for r in 0..h {
|
||
for c in 0..w {
|
||
let k = fdir[r * w + c];
|
||
if k < 0 {
|
||
continue;
|
||
}
|
||
let (dr, dc) = d8_offset(k as u8);
|
||
let nr = r as i32 + dr;
|
||
let nc = (c as i32 + dc).rem_euclid(w as i32) as usize;
|
||
if nr >= 0 && nr < h as i32 {
|
||
in_degree[nr as usize * w + nc] += 1;
|
||
}
|
||
}
|
||
}
|
||
let mut queue = VecDeque::new();
|
||
for (i, °) in in_degree.iter().enumerate().take(n) {
|
||
if deg == 0 {
|
||
queue.push_back(i);
|
||
}
|
||
}
|
||
let mut accum = vec![1i32; n];
|
||
while let Some(idx) = queue.pop_front() {
|
||
let r = idx / w;
|
||
let c = idx % w;
|
||
let k = fdir[idx];
|
||
if k < 0 {
|
||
continue;
|
||
}
|
||
let (dr, dc) = d8_offset(k as u8);
|
||
let nr = r as i32 + dr;
|
||
let nc = (c as i32 + dc).rem_euclid(w as i32) as usize;
|
||
if nr >= 0 && nr < h as i32 {
|
||
let ni = nr as usize * w + nc;
|
||
accum[ni] += accum[idx];
|
||
in_degree[ni] -= 1;
|
||
if in_degree[ni] == 0 {
|
||
queue.push_back(ni);
|
||
}
|
||
}
|
||
}
|
||
accum
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Tests
|
||
// ---------------------------------------------------------------------------
|
||
|
||
#[cfg(test)]
|
||
mod tests {
|
||
use super::*;
|
||
|
||
fn flat_grid(w: u32, h: u32, val: f32) -> Vec<f32> {
|
||
vec![val; (w * h) as usize]
|
||
}
|
||
|
||
fn slope_grid(w: u32, h: u32) -> Vec<f32> {
|
||
let n = (w * h) as usize;
|
||
(0..n)
|
||
.map(|i| {
|
||
let r = i / w as usize;
|
||
let c = i % w as usize;
|
||
1.0 - (r as f32 / h as f32 * 0.5 + c as f32 / w as f32 * 0.5)
|
||
})
|
||
.collect()
|
||
}
|
||
|
||
/// A grid with a single bowl-shaped depression in the middle — an
|
||
/// unambiguous single-basin fixture.
|
||
fn bowl_grid(w: u32, h: u32) -> Vec<f32> {
|
||
let n = (w * h) as usize;
|
||
let cx = w as f32 / 2.0;
|
||
let cy = h as f32 / 2.0;
|
||
let max_r = (cx.min(cy)).max(1.0);
|
||
(0..n)
|
||
.map(|i| {
|
||
let r = (i / w as usize) as f32;
|
||
let c = (i % w as usize) as f32;
|
||
let d = (((c - cx).powi(2) + (r - cy).powi(2)).sqrt() / max_r).min(1.0);
|
||
// High rim (0.9), low center (0.1) — a genuine bowl.
|
||
0.1 + d * 0.8
|
||
})
|
||
.collect()
|
||
}
|
||
|
||
fn default_climate() -> ClimateInputs {
|
||
ClimateInputs { moisture_q: 50 }
|
||
}
|
||
|
||
#[test]
|
||
fn flat_grid_produces_no_lakes() {
|
||
let elev = flat_grid(32, 16, 0.5);
|
||
let result = solve(&elev, 32, 16, 0.3, default_climate());
|
||
assert!(
|
||
result.basins.iter().all(|b| b.cells.is_empty()),
|
||
"a perfectly flat grid has no depressions to fill"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn bowl_grid_produces_at_least_one_lake_basin() {
|
||
let elev = bowl_grid(64, 32);
|
||
let result = solve(&elev, 64, 32, 0.0, default_climate());
|
||
let nonempty: Vec<_> = result
|
||
.basins
|
||
.iter()
|
||
.filter(|b| !b.cells.is_empty())
|
||
.collect();
|
||
assert!(
|
||
!nonempty.is_empty(),
|
||
"a bowl-shaped depression must fill to at least one lake basin"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn every_lake_cell_is_at_or_above_original_elevation() {
|
||
let elev = bowl_grid(64, 32);
|
||
let original: Vec<i64> = elev
|
||
.iter()
|
||
.map(|&e| (e as f64 * ELEV_SCALE) as i64)
|
||
.collect();
|
||
let result = solve(&elev, 64, 32, 0.0, default_climate());
|
||
for basin in &result.basins {
|
||
for &cell in &basin.cells {
|
||
assert!(
|
||
result.filled_scaled[cell] >= original[cell],
|
||
"filled elevation must never be below original terrain"
|
||
);
|
||
}
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn determinism_full_result() {
|
||
let elev = slope_grid(128, 64);
|
||
let r1 = solve(&elev, 128, 64, 0.3, default_climate());
|
||
let r2 = solve(&elev, 128, 64, 0.3, default_climate());
|
||
assert_eq!(
|
||
r1.filled_scaled, r2.filled_scaled,
|
||
"filled surface must be deterministic"
|
||
);
|
||
assert_eq!(
|
||
r1.channel_depth_scaled, r2.channel_depth_scaled,
|
||
"carved channel depth must be deterministic"
|
||
);
|
||
assert_eq!(
|
||
r1.cliff_edge, r2.cliff_edge,
|
||
"cliff-edge flags must be deterministic"
|
||
);
|
||
assert_eq!(
|
||
r1.basins.len(),
|
||
r2.basins.len(),
|
||
"basin count must be deterministic"
|
||
);
|
||
for (a, b) in r1.basins.iter().zip(r2.basins.iter()) {
|
||
assert_eq!(a.basin_id, b.basin_id);
|
||
assert_eq!(a.cells, b.cells);
|
||
assert_eq!(a.spill_level_scaled, b.spill_level_scaled);
|
||
assert_eq!(a.spill_cell, b.spill_cell);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn determinism_bowl_grid() {
|
||
// Second determinism fixture on a grid that actually produces lake
|
||
// basins (slope_grid alone may not) — exercises the overflow-path
|
||
// search + carving determinism, not just the trivial empty case.
|
||
let elev = bowl_grid(96, 48);
|
||
let r1 = solve(&elev, 96, 48, 0.1, default_climate());
|
||
let r2 = solve(&elev, 96, 48, 0.1, default_climate());
|
||
assert_eq!(r1.filled_scaled, r2.filled_scaled);
|
||
assert_eq!(r1.channel_depth_scaled, r2.channel_depth_scaled);
|
||
assert_eq!(r1.cliff_edge, r2.cliff_edge);
|
||
assert_eq!(r1.basins.len(), r2.basins.len());
|
||
for (a, b) in r1.basins.iter().zip(r2.basins.iter()) {
|
||
assert_eq!(format!("{:?}", a.outcome), format!("{:?}", b.outcome));
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn is_endorheic_large_dry_basin_is_endorheic() {
|
||
assert!(is_endorheic(0.01, 20));
|
||
}
|
||
|
||
#[test]
|
||
fn is_endorheic_tiny_basin_always_overflows() {
|
||
assert!(!is_endorheic(0.0001, 10));
|
||
}
|
||
|
||
#[test]
|
||
fn is_endorheic_wet_body_always_overflows() {
|
||
assert!(!is_endorheic(0.01, 90));
|
||
}
|
||
|
||
#[test]
|
||
fn overflowing_basin_has_nonempty_outlet_path() {
|
||
let elev = bowl_grid(64, 32);
|
||
// Force overflow classification via high moisture.
|
||
let result = solve(&elev, 64, 32, 0.0, ClimateInputs { moisture_q: 90 });
|
||
let overflowed: Vec<_> = result
|
||
.basins
|
||
.iter()
|
||
.filter(|b| !b.cells.is_empty())
|
||
.filter(|b| matches!(b.outcome, BasinOutcome::Overflow { .. }))
|
||
.collect();
|
||
assert!(
|
||
!overflowed.is_empty(),
|
||
"a wet-climate bowl basin must overflow, not go endorheic"
|
||
);
|
||
for basin in overflowed {
|
||
if let BasinOutcome::Overflow { outlet_path, .. } = &basin.outcome {
|
||
assert!(
|
||
!outlet_path.is_empty(),
|
||
"overflow basin must have a search path"
|
||
);
|
||
}
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn large_dry_basin_classifies_endorheic() {
|
||
let elev = bowl_grid(96, 48);
|
||
let result = solve(&elev, 96, 48, 0.0, ClimateInputs { moisture_q: 10 });
|
||
let has_endorheic = result
|
||
.basins
|
||
.iter()
|
||
.filter(|b| !b.cells.is_empty())
|
||
.any(|b| matches!(b.outcome, BasinOutcome::Endorheic { .. }));
|
||
assert!(
|
||
has_endorheic,
|
||
"a large, dry-climate bowl basin should classify endorheic under the \
|
||
moisture-governed criterion"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn single_basin_bowl_never_carves_a_gorge() {
|
||
// A single sealed basin's own rim is, by construction of
|
||
// priority-flood, always exactly the basin's spill level — so it
|
||
// never needs carving (see the module's "Gorge carving" test
|
||
// section below for the full structural explanation, and the
|
||
// T-1177 results doc's "Gorge carving — the honest story" section).
|
||
// This bowl fixture is exactly that shape (one basin, no chained
|
||
// second basin), so the correct, non-vacuous assertion here is that
|
||
// NOTHING gets carved — not an `if cliff_edge[i] { assert depth > 0
|
||
// }` loop, which passes vacuously whenever (as here) the condition
|
||
// never fires. The carving arithmetic ITSELF (what happens when a
|
||
// path cell IS above spill level) is covered independently by
|
||
// `carving_arithmetic_matches_original_minus_spill_level`.
|
||
let elev = bowl_grid(64, 32);
|
||
let result = solve(&elev, 64, 32, 0.0, ClimateInputs { moisture_q: 90 });
|
||
assert!(
|
||
result.cliff_edge.iter().all(|&c| !c),
|
||
"a single sealed basin's own rim is always exactly its spill level under this \
|
||
priority-flood algorithm, so it structurally never carves — see the module's \
|
||
\"Gorge carving\" test section for why"
|
||
);
|
||
assert!(
|
||
result.channel_depth_scaled.iter().all(|&d| d == 0),
|
||
"channel_depth_scaled must be all-zero whenever cliff_edge is all-false"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn basin_cells_are_sorted_ascending() {
|
||
let elev = bowl_grid(64, 32);
|
||
let result = solve(&elev, 64, 32, 0.0, default_climate());
|
||
for basin in &result.basins {
|
||
let mut sorted = basin.cells.clone();
|
||
sorted.sort_unstable();
|
||
assert_eq!(
|
||
basin.cells, sorted,
|
||
"basin cells must be in ascending row-major order"
|
||
);
|
||
}
|
||
}
|
||
|
||
// -----------------------------------------------------------------------
|
||
// Gorge carving — direct verification of `cheapest_overflow_path` +
|
||
// the carving-application arithmetic (T-1177).
|
||
// -----------------------------------------------------------------------
|
||
//
|
||
// **Investigation finding, kept here because it is genuinely load-
|
||
// bearing for the workshop's cliff-representation question, not just
|
||
// test-authoring trivia.** Several `solve()`-level fixture attempts (a
|
||
// single sealed basin behind a ridge, a basin chained into a second
|
||
// basin via a connecting corridor) were built to force a nonzero
|
||
// `channel_depth_scaled` end to end, and every one of them produced
|
||
// ZERO carved cells — not because carving is unreachable dead code, but
|
||
// because of a structural property of priority-flood itself:
|
||
// `priority_flood_fill` (this module) finds the TRUE global minimum
|
||
// enclosing rim of a basin, swallowing any monotonically-non-decreasing
|
||
// connected run of cells into the SAME lake and adopting that run's own
|
||
// peak as the basin's real spill level. Since `cheapest_overflow_path`'s
|
||
// termination check compares against that SAME spill level, the first
|
||
// cell it reaches outside the lake is — by construction of how the
|
||
// spill level was computed — never higher than spill level, so nothing
|
||
// gets carved. This holds for chained basins too, unless the connecting
|
||
// corridor is held to a single row/column (no multi-cell monotonic run
|
||
// for the flood to absorb) AND is higher than both basins' own
|
||
// independent rims (else the two basins simply merge into one).
|
||
//
|
||
// That combination is a genuine, narrow geometry — not a common one at
|
||
// continental (512×256 / 768×432 / 3840×2160) working-grid resolution,
|
||
// which is exactly why all three production-scale benches in
|
||
// `hydrology_equilibrium_bench.rs` report zero carved cells on real
|
||
// (GJ1c) and synthetic multi-octave terrain alike. That is an honest
|
||
// measurement result, not a bug — see the T-1177 results doc's "what
|
||
// this prototype does NOT do" section. The tests below verify the
|
||
// carving MECHANISM directly (bypassing `solve()`'s basin/spill
|
||
// detection, which is separately covered by the `determinism_*` and
|
||
// `basin_*` tests above) rather than continuing to hunt for a
|
||
// `solve()`-level fixture that exercises the narrow geometry above.
|
||
|
||
#[test]
|
||
fn cheapest_overflow_path_finds_the_true_minimum_crossing() {
|
||
// 5x5 grid. Basin cell 12 (row2,col2), spill_level 100_000. Every
|
||
// other cell starts as an impassable wall (500_000) EXCEPT: index
|
||
// 13 (row2,col3, the only path toward the far side) at 300_000
|
||
// (above spill level — the gorge cell this test expects the search
|
||
// to carve through), and index 14 (row2,col4) at 50_000 — a
|
||
// genuine below-spill-level terminus (open low ground) beyond it.
|
||
// The cheapest route MUST cross 13 (above spill, carve-worthy) to
|
||
// reach the valid terminus at 14.
|
||
let w = 5usize;
|
||
let h = 5usize;
|
||
let mut original = vec![500_000i64; w * h];
|
||
original[12] = 100_000; // spill cell (row2,col2) itself, basin-side
|
||
original[13] = 300_000; // the only crossing (row2,col3) — above spill level
|
||
original[14] = 50_000; // valid open-low-ground terminus (row2,col4)
|
||
let mut basin_of: Vec<Option<u32>> = vec![None; w * h];
|
||
basin_of[12] = Some(0);
|
||
let (path, target) = cheapest_overflow_path(12, 100_000, &original, &basin_of, 0, -1, w, h);
|
||
assert_eq!(
|
||
target,
|
||
DownstreamTarget::OpenSpillway,
|
||
"cell 14 is open low ground (not sea, not another basin) — a genuine SUCCESS \
|
||
terminus, OpenSpillway, not EdgeUnreachable (which is reserved for search \
|
||
exhaustion)"
|
||
);
|
||
assert_eq!(
|
||
path,
|
||
vec![12, 13, 14],
|
||
"the search must cross the carve-worthy cell 13 to reach terminus 14"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn cheapest_overflow_path_prefers_lower_total_cost_over_shorter_path() {
|
||
// 2-row, 20-column grid: the terminus sits at col 7, only 6 hops
|
||
// east of the spill (col 1) going the short way, but 13+ hops
|
||
// going the wrapped way around cols 0/19 — wide enough that the
|
||
// wrap-around is never cheaper (columns wrap horizontally, D8,
|
||
// matching `drainage.rs`; an earlier draft used a 9-column grid
|
||
// where the wrapped route was SHORTER and the cost-based search
|
||
// correctly, if unhelpfully-for-this-test, took it — not a bug,
|
||
// just an assumption this test needed to control for explicitly).
|
||
// Row 1 is an impassable wall (999_999) at every column, so the
|
||
// row-0 corridor (climbing to 50_000 at every intermediate column)
|
||
// is the only non-wrapped route.
|
||
let w = 20usize;
|
||
let h = 2usize;
|
||
let mut original = vec![999_999i64; w * h];
|
||
for cell in original.iter_mut().take(w) {
|
||
*cell = 50_000; // row 0, the corridor
|
||
}
|
||
original[1] = 0; // spill cell (row0,col1)
|
||
original[7] = -10_000; // sea-level terminus (row0,col7)
|
||
let mut basin_of: Vec<Option<u32>> = vec![None; w * h];
|
||
basin_of[1] = Some(0);
|
||
let (path, target) = cheapest_overflow_path(1, 0, &original, &basin_of, 0, -5_000, w, h);
|
||
assert_eq!(target, DownstreamTarget::Sea);
|
||
assert_eq!(
|
||
path,
|
||
vec![1, 2, 3, 4, 5, 6, 7],
|
||
"row 1 is walled off and the wrapped route is longer, so the direct row-0 \
|
||
corridor is the only cheapest route"
|
||
);
|
||
// Every intermediate cell is at 50_000, strictly above the spill
|
||
// level (0) — this is the exact shape `solve()`'s carving step
|
||
// consumes: `original[cell] - spill_level` for each path cell.
|
||
for &cell in &path[1..path.len() - 1] {
|
||
assert!(
|
||
original[cell] > 0,
|
||
"intermediate cells must be above spill level (spill level is 0 here)"
|
||
);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn carving_arithmetic_matches_original_minus_spill_level() {
|
||
// Direct check of the exact formula `solve()` uses to populate
|
||
// `channel_depth_scaled`/`cliff_edge` — reproduced here standalone
|
||
// so it's verified independent of whether any `solve()`-level
|
||
// fixture manages to trigger it end-to-end (see the module note
|
||
// above on why that geometry is narrow).
|
||
let spill_level: i64 = 500_000;
|
||
let path_elevations: [i64; 4] = [500_000, 620_000, 750_000, 480_000];
|
||
let mut channel_depth = vec![0i64; path_elevations.len()];
|
||
let mut cliff_edge = vec![false; path_elevations.len()];
|
||
for (i, &elev) in path_elevations.iter().enumerate() {
|
||
if elev > spill_level {
|
||
channel_depth[i] = elev - spill_level;
|
||
cliff_edge[i] = true;
|
||
}
|
||
}
|
||
assert_eq!(channel_depth, vec![0, 120_000, 250_000, 0]);
|
||
assert_eq!(cliff_edge, vec![false, true, true, false]);
|
||
}
|
||
|
||
// -----------------------------------------------------------------------
|
||
// T-1185 — basin-outlet -> D8 river-network wiring
|
||
// -----------------------------------------------------------------------
|
||
|
||
/// `bowl_grid` at `moisture_q = 100` (well above `ENDORHEIC_MOISTURE_
|
||
/// CEILING = 60`) — the SAME single-basin fixture `overflowing_basin_
|
||
/// has_nonempty_outlet_path` already relies on to force an `Overflow`
|
||
/// classification.
|
||
fn overflow_bowl_result() -> HydrologyResult {
|
||
let elev = bowl_grid(64, 32);
|
||
solve(&elev, 64, 32, 0.0, ClimateInputs { moisture_q: 100 })
|
||
}
|
||
|
||
/// `bowl_grid` at `moisture_q = 0` (well below the ceiling) — forces
|
||
/// `Endorheic`, same geometry as [`overflow_bowl_result`] so the two
|
||
/// fixtures differ ONLY in classification, not basin shape (mirrors
|
||
/// `run_layer1_with_moisture_changes_endorheic_split_not_lake_extent`'s
|
||
/// own straddle-the-ceiling pattern in `layer1.rs`).
|
||
fn endorheic_bowl_result() -> HydrologyResult {
|
||
let elev = bowl_grid(64, 32);
|
||
solve(&elev, 64, 32, 0.0, ClimateInputs { moisture_q: 0 })
|
||
}
|
||
|
||
fn empty_network() -> crate::atlas::body_world_state::RiverNetwork {
|
||
crate::atlas::body_world_state::RiverNetwork::default()
|
||
}
|
||
|
||
#[test]
|
||
fn overflow_basin_extends_river_network_with_a_new_course() {
|
||
let result = overflow_bowl_result();
|
||
assert!(
|
||
result
|
||
.basins
|
||
.iter()
|
||
.any(|b| matches!(b.outcome, BasinOutcome::Overflow { .. })),
|
||
"fixture sanity: the wet bowl must produce at least one Overflow basin"
|
||
);
|
||
let before = empty_network();
|
||
let after = extend_river_network_with_basin_outlets(&result, 64, 32, before.clone());
|
||
assert!(
|
||
after.river_cells.len() > before.river_cells.len(),
|
||
"an Overflow basin's outlet_path must add at least one new river cell — \
|
||
this IS the map's 'this lake drains' cue (D-227 amendment (4))"
|
||
);
|
||
assert_eq!(
|
||
after.river_cells.len(),
|
||
after.river_class.len(),
|
||
"river_cells/river_class must stay parallel arrays"
|
||
);
|
||
assert_eq!(
|
||
after.river_cells.len(),
|
||
after.river_downstream.len(),
|
||
"river_cells/river_downstream must stay parallel arrays"
|
||
);
|
||
assert_eq!(
|
||
after.river_cells.len(),
|
||
after.river_seaward.len(),
|
||
"river_cells/river_seaward must stay parallel arrays"
|
||
);
|
||
// Every appended downstream pointer must be a real D8 direction (0-7)
|
||
// or a documented terminus sentinel — never left unset / defaulted.
|
||
for &sentinel in &after.river_downstream {
|
||
assert!(
|
||
sentinel <= crate::atlas::body_world_state::RIVER_DOWNSTREAM_TERMINAL,
|
||
"unexpected river_downstream sentinel value: {sentinel}"
|
||
);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn overflow_basin_outlet_survives_build_edges_as_a_real_course() {
|
||
// The end-to-end guarantee the map actually reads: an Overflow
|
||
// basin's outlet must appear as a real `RiverEdge` once
|
||
// `river_course::build_edges` walks the extended network — not just
|
||
// as raw array entries on `RiverNetwork` itself.
|
||
let result = overflow_bowl_result();
|
||
let network = extend_river_network_with_basin_outlets(&result, 64, 32, empty_network());
|
||
let edges = crate::atlas::river_course::build_edges(&network);
|
||
assert!(
|
||
!edges.is_empty(),
|
||
"an Overflow basin's outlet_path must produce at least one buildable \
|
||
RiverEdge — the exit river the map needs to show"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn endorheic_basin_adds_no_course_the_positive_absence_cue() {
|
||
// The other half of the cue, asserted POSITIVELY (per the ticket):
|
||
// an Endorheic basin's lake must end up with NO outgoing course at
|
||
// all, not merely "fewer" — this IS the D-227 amendment (4) signal
|
||
// ("no course leaving this lake").
|
||
let result = endorheic_bowl_result();
|
||
assert!(
|
||
result
|
||
.basins
|
||
.iter()
|
||
.any(|b| matches!(b.outcome, BasinOutcome::Endorheic { .. })),
|
||
"fixture sanity: the dry bowl must produce at least one Endorheic basin"
|
||
);
|
||
let before = empty_network();
|
||
let after = extend_river_network_with_basin_outlets(&result, 64, 32, before.clone());
|
||
assert_eq!(
|
||
after.river_cells, before.river_cells,
|
||
"an Endorheic basin must contribute ZERO new river cells — the network \
|
||
is byte-identical to the unextended input"
|
||
);
|
||
assert_eq!(after.river_class, before.river_class);
|
||
assert_eq!(after.river_downstream, before.river_downstream);
|
||
assert_eq!(after.river_seaward, before.river_seaward);
|
||
// And the end-to-end read: build_edges over the (unchanged) network
|
||
// yields no edges either, since the input network started empty.
|
||
let edges = crate::atlas::river_course::build_edges(&after);
|
||
assert!(
|
||
edges.is_empty(),
|
||
"an Endorheic-only basin set must produce zero courses — 'no course \
|
||
leaves this lake' is the whole cue"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn basin_outlet_wiring_is_deterministic() {
|
||
// Two independent solves + wirings from the same inputs must be
|
||
// byte-identical (D-010) — the mandatory determinism gate every
|
||
// seed-chaining mechanism in this codebase carries.
|
||
let elev = bowl_grid(64, 32);
|
||
let climate = ClimateInputs { moisture_q: 100 };
|
||
let r1 = solve(&elev, 64, 32, 0.0, climate);
|
||
let r2 = solve(&elev, 64, 32, 0.0, climate);
|
||
let n1 = extend_river_network_with_basin_outlets(&r1, 64, 32, empty_network());
|
||
let n2 = extend_river_network_with_basin_outlets(&r2, 64, 32, empty_network());
|
||
assert_eq!(n1.river_cells, n2.river_cells);
|
||
assert_eq!(n1.river_class, n2.river_class);
|
||
assert_eq!(n1.river_downstream, n2.river_downstream);
|
||
assert_eq!(n1.river_seaward, n2.river_seaward);
|
||
}
|
||
|
||
#[test]
|
||
fn endorheic_basin_wiring_is_also_deterministic() {
|
||
let elev = bowl_grid(64, 32);
|
||
let climate = ClimateInputs { moisture_q: 0 };
|
||
let r1 = solve(&elev, 64, 32, 0.0, climate);
|
||
let r2 = solve(&elev, 64, 32, 0.0, climate);
|
||
let n1 = extend_river_network_with_basin_outlets(&r1, 64, 32, empty_network());
|
||
let n2 = extend_river_network_with_basin_outlets(&r2, 64, 32, empty_network());
|
||
assert_eq!(n1.river_cells, n2.river_cells);
|
||
assert_eq!(n1.river_downstream, n2.river_downstream);
|
||
}
|
||
|
||
#[test]
|
||
fn terminal_sentinel_on_a_preexisting_cell_is_still_a_safe_no_op_in_build_edges() {
|
||
// T-1170 Ruling 7b's own guarantee, re-exercised here at the
|
||
// integration seam this ticket wires up: a RIVER_DOWNSTREAM_TERMINAL
|
||
// entry anywhere in `river_downstream` must never break
|
||
// `build_edges`, regardless of whether the entry came from the
|
||
// original D8 extraction or (hypothetically) sat alongside T-1185's
|
||
// appended cells. `extend_river_network_with_basin_outlets` itself
|
||
// never emits TERMINAL (Endorheic basins are skipped outright, never
|
||
// mapped to the sentinel) — this test pins that build_edges' existing
|
||
// no-op handling still holds on a network T-1185's function has
|
||
// touched, not just a hand-built one.
|
||
let result = overflow_bowl_result();
|
||
let mut network = extend_river_network_with_basin_outlets(&result, 64, 32, empty_network());
|
||
assert!(
|
||
!network.river_cells.is_empty(),
|
||
"fixture sanity: the overflow bowl must have appended at least one cell"
|
||
);
|
||
// Manually inject a TERMINAL entry (simulating a future endorheic
|
||
// interior-sink emission from elsewhere in the pipeline, per Ruling
|
||
// 7b's reservation) at the first appended cell.
|
||
network.river_downstream[0] = crate::atlas::body_world_state::RIVER_DOWNSTREAM_TERMINAL;
|
||
let edges = crate::atlas::river_course::build_edges(&network);
|
||
// No panic, and the TERMINAL cell itself produces no edge (build_edges'
|
||
// own documented no-op) — the remaining cells (if any) are unaffected.
|
||
assert!(
|
||
edges.iter().all(|e| e.upstream != network.river_cells[0]),
|
||
"a RIVER_DOWNSTREAM_TERMINAL cell must never produce a RiverEdge"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn extend_river_network_skips_cells_already_present() {
|
||
// A hand-built network that already contains the bowl's expected
|
||
// outlet cell(s) must not gain duplicate parallel-array entries —
|
||
// the "re-entered the pre-existing D8 network" branch.
|
||
let result = overflow_bowl_result();
|
||
let overflow = result
|
||
.basins
|
||
.iter()
|
||
.find(|b| matches!(b.outcome, BasinOutcome::Overflow { .. }))
|
||
.expect("fixture sanity: at least one Overflow basin");
|
||
let BasinOutcome::Overflow { outlet_path, .. } = &overflow.outcome else {
|
||
unreachable!()
|
||
};
|
||
assert!(
|
||
outlet_path.len() >= 2,
|
||
"fixture sanity: outlet_path must have at least one interior cell \
|
||
to pre-seed as already-present"
|
||
);
|
||
let pre_existing_cell = outlet_path[1];
|
||
let (r, c) = (pre_existing_cell / 64, pre_existing_cell % 64);
|
||
let mut network = empty_network();
|
||
network.river_cells.push((r as u16, c as u16));
|
||
network.river_class.push(1);
|
||
network
|
||
.river_downstream
|
||
.push(crate::atlas::body_world_state::RIVER_DOWNSTREAM_EDGE_DRAIN);
|
||
network.river_seaward.push((0, 0));
|
||
|
||
let after = extend_river_network_with_basin_outlets(&result, 64, 32, network);
|
||
// The pre-existing cell must appear exactly once, not twice.
|
||
let occurrences = after
|
||
.river_cells
|
||
.iter()
|
||
.filter(|&&cell| cell == (r as u16, c as u16))
|
||
.count();
|
||
assert_eq!(
|
||
occurrences, 1,
|
||
"a basin-outlet cell already present in the network must not be duplicated"
|
||
);
|
||
assert_eq!(
|
||
after.river_cells.len(),
|
||
after.river_class.len(),
|
||
"parallel arrays must stay in sync even when a duplicate is skipped"
|
||
);
|
||
}
|
||
|
||
/// PR #202 review, Hoshe finding 1 — the i==1 collision regression,
|
||
/// asserted NON-VACUOUSLY on the actual D-227 amendment (4) cue: when
|
||
/// `outlet_path[1]` (the FIRST interior cell after the spill) is already
|
||
/// a pre-existing river cell, the basin's Overflow status must still be
|
||
/// READABLE via `build_edges` — a course must leave the lake. Before the
|
||
/// fix, this exact scenario made the loop `break` on its first
|
||
/// iteration having pushed zero cells for the basin, silently making a
|
||
/// real Overflow basin indistinguishable from Endorheic.
|
||
#[test]
|
||
fn overflow_basin_survives_an_i_equals_one_collision_non_vacuously() {
|
||
let result = overflow_bowl_result();
|
||
let overflow = result
|
||
.basins
|
||
.iter()
|
||
.find(|b| matches!(b.outcome, BasinOutcome::Overflow { .. }))
|
||
.expect("fixture sanity: at least one Overflow basin");
|
||
let BasinOutcome::Overflow { outlet_path, .. } = &overflow.outcome else {
|
||
unreachable!()
|
||
};
|
||
assert!(
|
||
outlet_path.len() >= 2,
|
||
"fixture sanity: outlet_path must have an interior cell at index 1 \
|
||
to pre-seed as the i==1 collision"
|
||
);
|
||
let spill = outlet_path[0];
|
||
let collision_cell = outlet_path[1]; // the exact i==1 collision position
|
||
let (sr, sc) = (spill / 64, spill % 64);
|
||
let (cr, cc) = (collision_cell / 64, collision_cell % 64);
|
||
|
||
// Pre-seed the network with ONLY the i==1 cell as a native river
|
||
// cell (pointing at some unrelated existing direction) — the spill
|
||
// cell itself is NOT pre-existing, isolating this as the exact
|
||
// "first interior cell collides" scenario Hoshe's finding names.
|
||
let mut network = empty_network();
|
||
network.river_cells.push((cr as u16, cc as u16));
|
||
network.river_class.push(1);
|
||
network
|
||
.river_downstream
|
||
.push(crate::atlas::body_world_state::RIVER_DOWNSTREAM_EDGE_DRAIN);
|
||
network.river_seaward.push((0, 0));
|
||
|
||
let after = extend_river_network_with_basin_outlets(&result, 64, 32, network);
|
||
|
||
// The spill cell must have its OWN entry — this is the actual fix:
|
||
// before it, nothing was ever pushed for this basin at all.
|
||
let spill_idx = after
|
||
.river_cells
|
||
.iter()
|
||
.position(|&cell| cell == (sr as u16, sc as u16))
|
||
.expect(
|
||
"the spill cell must have a real river_cells entry even when \
|
||
outlet_path[1] collides with a pre-existing river cell — this \
|
||
is the exact bug PR #202 flagged: a silent break at i==1 must \
|
||
never erase the basin's own outflow anchor",
|
||
);
|
||
// Its downstream pointer must be a real D8 direction (not a sentinel)
|
||
// toward the collision cell, since outlet_path[1] IS the next hop.
|
||
assert!(
|
||
after.river_downstream[spill_idx] < 8,
|
||
"the spill cell's downstream pointer must be a real D8 direction \
|
||
toward outlet_path[1], not a terminus sentinel"
|
||
);
|
||
|
||
// Non-vacuous on the cue itself: build_edges must produce a real
|
||
// edge whose upstream IS the spill cell — an overflow lake's exit
|
||
// river must be readable from the served network, not just present
|
||
// as raw unreachable array data.
|
||
let edges = crate::atlas::river_course::build_edges(&after);
|
||
assert!(
|
||
edges.iter().any(|e| e.upstream == (sr as u16, sc as u16)),
|
||
"an Overflow basin's spill cell must produce a real RiverEdge even \
|
||
when its first interior outlet cell collides with a pre-existing \
|
||
river cell — otherwise this basin is visually indistinguishable \
|
||
from Endorheic, which is exactly the cue this ticket must not break"
|
||
);
|
||
}
|
||
|
||
/// The sibling collision case: the SPILL CELL ITSELF (not an interior
|
||
/// cell) is already a pre-existing river cell — verified to occur on
|
||
/// real GJ1c data (2/51 real Overflow basins). The pre-existing entry's
|
||
/// downstream pointer must be overwritten to point along the outlet
|
||
/// path (the hydrology solve is the more authoritative answer for this
|
||
/// cell's true direction), not left pointing wherever the original flat
|
||
/// D8 extraction guessed — and it must not create a duplicate
|
||
/// `river_cells` entry at the same position.
|
||
#[test]
|
||
fn overflow_basin_overwrites_a_pre_existing_spill_cell_in_place() {
|
||
let result = overflow_bowl_result();
|
||
let overflow = result
|
||
.basins
|
||
.iter()
|
||
.find(|b| matches!(b.outcome, BasinOutcome::Overflow { .. }))
|
||
.expect("fixture sanity: at least one Overflow basin");
|
||
let BasinOutcome::Overflow { outlet_path, .. } = &overflow.outcome else {
|
||
unreachable!()
|
||
};
|
||
assert!(
|
||
outlet_path.len() >= 2,
|
||
"fixture sanity: needs an interior cell"
|
||
);
|
||
let spill = outlet_path[0];
|
||
let (sr, sc) = (spill / 64, spill % 64);
|
||
|
||
// Pre-seed the network with the SPILL cell as a native river cell,
|
||
// pointing in an unrelated direction (simulating the original flat
|
||
// D8 extraction having already claimed this cell for its own,
|
||
// wrong-for-the-lake reasons).
|
||
let mut network = empty_network();
|
||
network.river_cells.push((sr as u16, sc as u16));
|
||
network.river_class.push(2); // deliberately a different class
|
||
network.river_downstream.push(0); // deliberately a real but likely-wrong direction
|
||
network.river_seaward.push((0, 0));
|
||
|
||
let after = extend_river_network_with_basin_outlets(&result, 64, 32, network);
|
||
|
||
// Exactly one entry at the spill position — overwritten, not duplicated.
|
||
let occurrences = after
|
||
.river_cells
|
||
.iter()
|
||
.filter(|&&cell| cell == (sr as u16, sc as u16))
|
||
.count();
|
||
assert_eq!(
|
||
occurrences, 1,
|
||
"the spill cell must be overwritten in place, never duplicated"
|
||
);
|
||
let spill_idx = after
|
||
.river_cells
|
||
.iter()
|
||
.position(|&cell| cell == (sr as u16, sc as u16))
|
||
.unwrap();
|
||
assert!(
|
||
after.river_downstream[spill_idx] < 8,
|
||
"the overwritten spill cell must point along the real outlet path, \
|
||
not retain whatever direction the pre-existing entry had"
|
||
);
|
||
|
||
// Non-vacuous on the cue: the overwritten spill cell must still
|
||
// build into a real, readable RiverEdge.
|
||
let edges = crate::atlas::river_course::build_edges(&after);
|
||
assert!(
|
||
edges.iter().any(|e| e.upstream == (sr as u16, sc as u16)),
|
||
"an overwritten spill cell must still produce a real RiverEdge — \
|
||
the lake's outflow must remain readable"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn d8_direction_between_matches_the_offset_table() {
|
||
// Round-trip sanity: for every direction in the D8 table, stepping
|
||
// from an interior cell and asking for the direction back must
|
||
// recover the same index.
|
||
let w = 16usize;
|
||
let from = (8usize, 8usize);
|
||
for k in 0..D8_LEN {
|
||
let (dr, dc) = d8_offset(k);
|
||
let to = (
|
||
(from.0 as i32 + dr) as usize,
|
||
(from.1 as i32 + dc).rem_euclid(w as i32) as usize,
|
||
);
|
||
assert_eq!(
|
||
d8_direction_between(from, to, w),
|
||
k,
|
||
"round-trip mismatch for D8 direction {k}"
|
||
);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn d8_direction_between_handles_horizontal_wraparound() {
|
||
// Column 0 -> column (w-1) is a valid D8 West step under horizontal
|
||
// wraparound (the equirectangular grid convention every D8 walk in
|
||
// this module already uses) — must resolve to West (k=3), not panic.
|
||
let w = 16usize;
|
||
let from = (5usize, 0usize);
|
||
let to = (5usize, w - 1);
|
||
assert_eq!(
|
||
d8_direction_between(from, to, w),
|
||
3,
|
||
"expected West (wrapped)"
|
||
);
|
||
}
|
||
}
|