settled-reach/server/src/atlas/layer1.rs

//! Layer 1 orchestrator — empty-world topography (#953).
//!
//! Runs the full Layer-1 pipeline for one body, in order:
//!   1. D8 priority-flood drainage (D-208) → river network + basins
//!   2. shared terrain analysis (ocean/lake masks, water distance, slope,
//!      elevation percentile) — D-209/D-210 inputs
//!   3. 7-tag geographic feature extraction (D-209)
//!   4. sub-biome + terrain_modification_cost classification (D-210)
//!   5. per-district dominant D8 basin direction (T-1047, D-239 §8)
//!
//! Output is the in-memory `Layer1Output`, which maps directly onto
//! `BodyWorldState` (D-203). Name attachment (D-223) is a separate, cheap step
//! (`attach_feature_names`) so the compute can be benchmarked in isolation and
//! names sourced from the DB pool independently.
//!
//! `run_layer1` returns `(Layer1Output, TerrainAnalysis)` so `cascade.rs` can
//! reuse the `TerrainAnalysis` held on `CascadeSnapshot.terrain_analysis`
//! (transient — dropped after DistrictProfile + RoadGraph consume it; D-203 /
//! T-1044) without re-running the ~45 ms drainage pass per body.
//!
//! **Determinism (D-010 #4):** every stage is deterministic; the same heightmap
//! yields bit-identical attractors and river networks.

use std::collections::BTreeMap;

use crate::atlas::body_world_state::{DrainageBasin, RiverNetwork};
use crate::atlas::drainage::{self, DrainageResult};
use crate::atlas::features::{self, TerrainAnalysis};
use crate::atlas::heightmap::BodyHeightmap;
use crate::atlas::scale::{BasinDirection, SurveyCellPos, HEIGHTMAP_CELLS_PER_DISTRICT};
use crate::atlas::subbiome;
use crate::simulation::generator::{AttractorType, GeographicAttractor};
use serde::{Deserialize, Serialize};

/// Full Layer-1 result for one body.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct Layer1Output {
    pub body_id: String,
    pub river_network: RiverNetwork,
    pub drainage_basins: Vec<DrainageBasin>,
    /// Geographic attractors (D-209) with sub-biome + cost (D-210), sorted by
    /// `(attractor_type, row, col)`.
    pub attractors: Vec<GeographicAttractor>,
    /// Working-grid dimensions every position in this output (river cells, basin
    /// boundaries, attractor positions) is expressed in — equals the downsampled
    /// heightmap size. The client maps these onto the displayed heightmap (#960),
    /// so the overlay scale stays correct for any source resolution (mod-safe).
    pub grid_w: u32,
    pub grid_h: u32,
    /// Dominant D8 thalweg direction per SURVEY CELL, aggregated from the `fdir`
    /// grid during the Layer-1 drainage pass (T-1047, D-239 §8). Each entry
    /// holds the cardinal direction with the most votes among non-ocean cells in
    /// that cell's covering 8×8 working-grid pixel block. Keyed by
    /// [`SurveyCellPos`] (D-256(b)) using `HEIGHTMAP_CELLS_PER_DISTRICT` as the
    /// grid-to-cell mapping — the SAME survey raster `derive_all_districts`
    /// builds `DistrictProfile`s over, not the true D-243 `DistrictPos` grid
    /// (D-256: this field predates the newtype and escaped the initial sweep;
    /// the aggregate is honestly a survey-cell aggregate — it votes over
    /// exactly the pixel block one `DistrictProfile` summarizes — so
    /// `SurveyCellPos` is its correct, not just convenient, key).
    ///
    /// The VALUES are the **true D8-computed direction** — not a seed-bit
    /// proxy — so `DistrictProfile.basin_direction` (and downstream
    /// `ChunkContext`) respect drainage monotonicity (D-239 §8: respect the D8
    /// thalweg). Only the KEY space is the coarse survey raster.
    ///
    /// **Transient:** skipped in serialization (`#[serde(skip)]`) — this field is
    /// a cascade-internal transport from `run_layer1` to `derive_all_districts`
    /// and is re-derived on each `run_layer1` call. The per-survey-cell direction
    /// is persisted on `DistrictProfile.basin_direction` (`BodyWorldState.districts`)
    /// after the cascade consumes it.
    #[serde(skip)]
    pub survey_basin_dirs: BTreeMap<SurveyCellPos, BasinDirection>,
}

/// Body-wide moisture ceiling fallback for [`run_layer1`]'s hydrology solve
/// when no [`crate::atlas::district_profile::BodyParams`] is available to
/// derive a real one from (`run_layer1`'s signature is heightmap-only, matching
/// `drainage::analyze`'s own "same way it already runs once per body today"
/// shape per T-1177's scope). Matches the T-1177 prototype's own
/// `ClimateInputs { moisture_q: 55 }` population-survey default (moderate
/// hydrosphere, breathable atmosphere) — a reasonable body-agnostic guess,
/// used ONLY by [`run_layer1`]'s two-arg form; [`run_layer1_with_moisture`]
/// (called by every production site that has real `BodyParams` in scope) never
/// reaches this constant.
const DEFAULT_HYDROLOGY_MOISTURE_Q: i32 = 55;

/// Run the Layer-1 topography pipeline for a single body.
///
/// Returns `(Layer1Output, TerrainAnalysis)`. The `TerrainAnalysis` is carried
/// transiently on `CascadeSnapshot.terrain_analysis` so `cascade.rs` can pass
/// it to `derive_all_districts` and `build_road_graph` without re-running the
/// full D8 drainage pass (T-1044 — eliminates the PERF/TODO re-run).
///
/// Solves settled-equilibrium hydrology (T-1177/T-1184, D-227 amendment (4))
/// once per body as part of this same pass, using
/// [`DEFAULT_HYDROLOGY_MOISTURE_Q`] as the body-wide moisture ceiling — this
/// two-arg form has no `BodyParams` to derive a real one from. Every
/// production call site that DOES have `BodyParams` in scope
/// (`cascade::run_cascade_from_heightmap`,
/// `gen_queue::TerrainAnalysisCache::get_or_derive`) calls
/// [`run_layer1_with_moisture`] instead, so this fallback is only ever
/// exercised by call sites (mostly tests) that never had body params to begin
/// with — never a silent downgrade of a real value.
pub fn run_layer1(hm: &BodyHeightmap) -> (Layer1Output, TerrainAnalysis) {
    run_layer1_with_moisture(hm, DEFAULT_HYDROLOGY_MOISTURE_Q)
}

/// [`run_layer1`], with the body-wide hydrology moisture ceiling
/// (`ClimateInputs::moisture_q`, T-1177) supplied explicitly rather than
/// defaulted. Callers with a real [`crate::atlas::district_profile::BodyParams`]
/// in scope should derive it via
/// [`crate::atlas::district_profile::derive_moisture_ceiling_q`] and pass the
/// result here, so the endorheic-vs-overflow basin split reflects the body's
/// actual hydrosphere/atmosphere instead of the fallback constant.
///
/// Determinism (D-010): pure function of `(hm, moisture_q)` — same inputs,
/// byte-identical `TerrainAnalysis.hydrology` every time (inherits
/// `hydrology_equilibrium::solve`'s own determinism guarantee).
pub fn run_layer1_with_moisture(
    hm: &BodyHeightmap,
    moisture_q: i32,
) -> (Layer1Output, TerrainAnalysis) {
    let drainage: DrainageResult = drainage::analyze(&hm.data, hm.width, hm.height, hm.sea_level);
    let mut ta: TerrainAnalysis = TerrainAnalysis::analyze(hm, &drainage);

    // T-1184: solve settled-equilibrium hydrology once per body (the
    // AnalyzeBody cascade populate point, D-227 amendment (4)) and fold the
    // continuous filled-surface field into this TerrainAnalysis so every
    // derive-core caller downstream (`derive_at_metres_with_riparian`) can
    // bilinearly sample it for lake sourcing — mechanism B, D-255(f): a
    // coarse continuous primitive computed once, sampled fresh at every rung,
    // never re-solved. `hm.data` (the raw [0,1] elevation this analysis was
    // built from) is the SAME grid the solver runs on, so the two fields
    // `with_hydrology` stores are always the correct pairing.
    let hydrology = crate::atlas::hydrology_equilibrium::solve(
        &hm.data,
        hm.width,
        hm.height,
        hm.sea_level,
        crate::atlas::hydrology_equilibrium::ClimateInputs { moisture_q },
    );
    ta = ta.with_hydrology(&hm.data, &hydrology);

    let raw = features::extract_attractors(hm, &drainage, &ta);
    let attractors: Vec<GeographicAttractor> = raw
        .iter()
        .map(|r| {
            let (sub_biome, terrain_modification_cost) =
                subbiome::classify(&ta, r.row as usize, r.col as usize);
            GeographicAttractor {
                position: (r.row, r.col),
                attractor_type: r.attractor_type,
                strength: r.strength,
                sub_biome,
                terrain_modification_cost,
                // Layer-1 water-direction extraction (#957, D-234) — feeds D-213
                // founding orientation + the D-234 waterfront rule.
                water_bearing: ta.water_bearing(r.row as usize, r.col as usize),
            }
        })
        .collect();

    // Aggregate per-survey-cell dominant D8 direction from the fdir grid
    // (T-1047, D-239 §8; D-256(b) survey raster — NOT the true D-243 district
    // grid). fdir is available here before it is discarded — do NOT expose
    // the full grid on DrainageResult externally. The compact per-cell map
    // (~6 000 entries) is what propagates into Layer1Output and DistrictProfile.
    //
    // Mapping fdir index → 4-way cardinal (D-010 integer; matches D8 table):
    //   0 N, 1 S, 2 E, 3 W  (pure cardinals)
    //   4 NE → N (|dr|=|dc|=1; row component wins per D8 priority order)
    //   5 NW → N
    //   6 SE → S
    //   7 SW → S
    //  -1 → skip (no outflow: edge, flat peak, ocean)
    let survey_basin_dirs =
        aggregate_survey_basin_dirs(&drainage.fdir, hm.width, hm.height, &ta.ocean_mask);

    let l1 = Layer1Output {
        body_id: hm.body_id.clone(),
        river_network: drainage.river_network,
        drainage_basins: drainage.drainage_basins,
        attractors,
        grid_w: hm.width,
        grid_h: hm.height,
        survey_basin_dirs,
    };
    (l1, ta)
}

/// Aggregate a per-survey-cell dominant D8 flow direction from the full-grid
/// `fdir` (index into the D8 table, -1 = no outflow). Ocean-masked cells are
/// excluded from voting so coastal cells do not skew toward the ocean sink
/// direction.
///
/// Each non-ocean, non-sink cell casts one vote for its cardinal direction
/// (diagonals NE/NW fold to N, SE/SW fold to S). Ties broken by cardinal
/// precedence (N > S > E > W). Survey cells with no valid votes default to
/// `North`.
///
/// Integer arithmetic throughout (D-010).
fn aggregate_survey_basin_dirs(
    fdir: &[i8],
    width: u32,
    height: u32,
    ocean_mask: &[bool],
) -> BTreeMap<SurveyCellPos, BasinDirection> {
    let w = width as usize;
    let h = height as usize;
    let gcpd = HEIGHTMAP_CELLS_PER_DISTRICT;

    // Per-survey-cell vote counts: [N, S, E, W].
    let mut votes: BTreeMap<SurveyCellPos, [i32; 4]> = BTreeMap::new();

    for r in 0..h {
        for c in 0..w {
            let i = r * w + c;
            let k = fdir[i];
            if k < 0 || ocean_mask[i] {
                continue; // no-outflow or ocean — skip
            }
            // Map D8 index to 4-way cardinal vote index: [N=0, S=1, E=2, W=3].
            let vote = match k {
                0 => 0, // N
                1 => 1, // S
                2 => 2, // E
                3 => 3, // W
                4 => 0, // NE → N (row component wins; |dr|=|dc|=1)
                5 => 0, // NW → N
                6 => 1, // SE → S
                7 => 1, // SW → S
                _ => continue,
            };
            let cell_pos = SurveyCellPos((c / gcpd) as i32, (r / gcpd) as i32);
            votes.entry(cell_pos).or_insert([0i32; 4])[vote] += 1;
        }
    }

    // For each survey cell, pick the cardinal with the most votes.
    // Tie-breaking order: N > S > E > W (matches D8 priority).
    let survey_cols = w.div_ceil(gcpd) as i32;
    let survey_rows = h.div_ceil(gcpd) as i32;
    let mut out = BTreeMap::new();
    for dy in 0..survey_rows {
        for dx in 0..survey_cols {
            let pos = SurveyCellPos(dx, dy);
            let dir = if let Some(v) = votes.get(&pos) {
                // N=0, S=1, E=2, W=3 in descending priority for tie-breaking.
                let mut best_votes = -1i32;
                let mut best_dir = BasinDirection::North;
                for (cardinal_idx, &count) in v.iter().enumerate() {
                    // Strictly greater-than preserves the first (highest-priority)
                    // cardinal in case of tie.
                    if count > best_votes {
                        best_votes = count;
                        best_dir = match cardinal_idx {
                            0 => BasinDirection::North,
                            1 => BasinDirection::South,
                            2 => BasinDirection::East,
                            _ => BasinDirection::West,
                        };
                    }
                }
                best_dir
            } else {
                BasinDirection::North // ocean-only or empty cell: default
            };
            out.insert(pos, dir);
        }
    }
    out
}

/// Attach pool names (D-223) to the largest computed rivers and mountains.
///
/// Rivers are ranked by mouth strength (a proxy for catchment size) descending;
/// `RiverMouth` attractors take names from `river_names` in that order. Mountain
/// names attach to the highest-elevation `Alpine`/`PassEntrance` attractors.
/// Returns `(river_assignments, mountain_assignments)` as `(position, name)`
/// pairs; positions that outrun the pool get no name (the pool is finite).
pub fn attach_feature_names(
    output: &Layer1Output,
    river_names: &[String],
    mountain_names: &[String],
) -> (Vec<((u16, u16), String)>, Vec<((u16, u16), String)>) {
    // Rivers: RiverMouth attractors, strongest first (ties by row, col).
    let mut mouths: Vec<&GeographicAttractor> = output
        .attractors
        .iter()
        .filter(|a| a.attractor_type == AttractorType::RiverMouth)
        .collect();
    mouths.sort_by(|a, b| {
        // strength is integer now — rank directly (descending).
        b.strength
            .cmp(&a.strength)
            .then(a.position.0.cmp(&b.position.0))
            .then(a.position.1.cmp(&b.position.1))
    });
    let rivers = mouths
        .iter()
        .zip(river_names.iter())
        .map(|(a, n)| (a.position, n.clone()))
        .collect();

    // Mountains: Alpine attractors, strongest first.
    let mut peaks: Vec<&GeographicAttractor> = output
        .attractors
        .iter()
        .filter(|a| {
            matches!(
                a.sub_biome,
                crate::simulation::generator::SubBiomeVariant::Alpine
            )
        })
        .collect();
    peaks.sort_by(|a, b| {
        // strength is integer now — rank directly (descending).
        b.strength
            .cmp(&a.strength)
            .then(a.position.0.cmp(&b.position.0))
            .then(a.position.1.cmp(&b.position.1))
    });
    let mountains = peaks
        .iter()
        .zip(mountain_names.iter())
        .map(|(a, n)| (a.position, n.clone()))
        .collect();

    (rivers, mountains)
}

#[cfg(test)]
mod tests {
    use super::*;

    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()
    }

    fn hm(w: u32, h: u32) -> BodyHeightmap {
        BodyHeightmap {
            body_id: "TestBody".into(),
            width: w,
            height: h,
            data: slope_grid(w, h),
            sea_level: 0.3,
        }
    }

    #[test]
    fn run_layer1_is_deterministic() {
        let h = hm(128, 64);
        let (o1, _ta1) = run_layer1(&h);
        let (o2, _ta2) = run_layer1(&h);
        assert_eq!(o1.attractors.len(), o2.attractors.len());
        for (a, b) in o1.attractors.iter().zip(o2.attractors.iter()) {
            assert_eq!(a.position, b.position);
            assert_eq!(a.attractor_type, b.attractor_type);
            assert_eq!(a.strength, b.strength);
            assert_eq!(a.sub_biome, b.sub_biome);
            assert_eq!(a.terrain_modification_cost, b.terrain_modification_cost);
        }
        assert_eq!(o1.river_network.river_cells, o2.river_network.river_cells);
        // survey_basin_dirs is deterministic and non-empty on a slope grid.
        assert_eq!(o1.survey_basin_dirs, o2.survey_basin_dirs);
        assert!(
            !o1.survey_basin_dirs.is_empty(),
            "slope grid must produce survey-cell basin directions"
        );
    }

    #[test]
    fn produces_attractors_and_costs() {
        let (o, _ta) = run_layer1(&hm(256, 128));
        assert!(!o.attractors.is_empty(), "expected some attractors");
        assert!(o
            .attractors
            .iter()
            .all(|a| a.terrain_modification_cost >= 100));
        assert!(o.attractors.iter().all(|a| (0..=100).contains(&a.strength)));
    }

    #[test]
    fn name_attachment_respects_pool_size() {
        let (o, _ta) = run_layer1(&hm(256, 128));
        let names = vec!["Aldren".to_string(), "Brook".to_string()];
        let (rivers, _mtn) = attach_feature_names(&o, &names, &[]);
        assert!(rivers.len() <= names.len());
    }

    // -------------------------------------------------------------------
    // T-1184 — hydrology productionization
    // -------------------------------------------------------------------

    /// Bowl-shaped fixture (high rim, low centre) — the same shape
    /// `hydrology_equilibrium.rs`'s own `bowl_grid` test fixture uses,
    /// reproduced here (not imported — that one is `#[cfg(test)]`-private to
    /// its own module) so `run_layer1`'s hydrology wiring can be exercised
    /// end-to-end without depending on solver-internal test helpers.
    /// `sea_level: 0.0` keeps the whole grid land except the filled basin, so
    /// a resulting `MorphologyZone::Lake` can only be hydrology-sourced, never
    /// the `ocean_fraction_q` heuristic fallback.
    fn bowl_hm(w: u32, h: u32, body_id: &str) -> BodyHeightmap {
        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);
        let data = (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);
                0.1 + d * 0.8
            })
            .collect();
        BodyHeightmap {
            body_id: body_id.into(),
            width: w,
            height: h,
            data,
            sea_level: 0.0,
        }
    }

    #[test]
    fn run_layer1_populates_hydrology_on_terrain_analysis() {
        let h = bowl_hm(64, 32, "BowlBody");
        let (_o, ta) = run_layer1(&h);
        let hydro = ta
            .hydrology
            .as_ref()
            .expect("run_layer1 must populate TerrainAnalysis.hydrology (T-1184)");
        assert_eq!(hydro.elevation.len(), (64 * 32) as usize);
        assert_eq!(hydro.filled.len(), (64 * 32) as usize);
        // The bowl centre must be a lake cell: filled strictly exceeds original.
        let centre_idx = (16 * 64 + 32) as usize; // row 16, col 32 — the bowl centre
        assert!(
            hydro.filled[centre_idx] > hydro.elevation[centre_idx],
            "bowl centre must be filled above its original elevation"
        );
    }

    #[test]
    fn run_layer1_hydrology_is_deterministic() {
        let h = bowl_hm(64, 32, "BowlBody");
        let (_o1, ta1) = run_layer1(&h);
        let (_o2, ta2) = run_layer1(&h);
        let h1 = ta1.hydrology.expect("first run must populate hydrology");
        let h2 = ta2.hydrology.expect("second run must populate hydrology");
        assert_eq!(
            h1.elevation, h2.elevation,
            "D-010: identical inputs must produce byte-identical elevation carry"
        );
        assert_eq!(
            h1.filled, h2.filled,
            "D-010: identical inputs must produce byte-identical filled-surface field"
        );
    }

    /// The D-255(f) mandatory determinism gate: whether hydrology is solved
    /// with the fallback default moisture (`run_layer1`) or an explicit
    /// caller-supplied moisture that happens to equal the default
    /// (`run_layer1_with_moisture`), the two code paths must produce
    /// byte-identical `TerrainAnalysis.hydrology` output — the "cache-hit
    /// path == cache-miss path" shape applied to the two entry points that
    /// stand in for it here (both are genuinely fresh `derive()` calls; T-1184
    /// has no separate cached-coarser-canvas to compare against yet, since
    /// that tier is T-1181's rung-0 scope — this test instead pins that
    /// `run_layer1`'s convenience wrapper and its explicit-moisture sibling
    /// never silently diverge, which is the property the next ticket's
    /// resident-cache read will depend on staying true).
    #[test]
    fn run_layer1_default_and_explicit_moisture_agree_at_the_default_value() {
        let h = bowl_hm(64, 32, "BowlBody");
        let (_o1, ta1) = run_layer1(&h);
        let (_o2, ta2) = run_layer1_with_moisture(&h, DEFAULT_HYDROLOGY_MOISTURE_Q);
        let h1 = ta1.hydrology.expect("run_layer1 must populate hydrology");
        let h2 = ta2
            .hydrology
            .expect("run_layer1_with_moisture must populate hydrology");
        assert_eq!(h1.elevation, h2.elevation);
        assert_eq!(h1.filled, h2.filled);
    }

    #[test]
    fn run_layer1_with_moisture_changes_endorheic_split_not_lake_extent() {
        // T-1184 scope note (ticket text): moisture affects the
        // endorheic-vs-overflow split only, never lake EXTENT (filled_scaled
        // is a pure function of elevation/sea_level, moisture-independent).
        let h = bowl_hm(64, 32, "BowlBody");
        let (_o_dry, ta_dry) = run_layer1_with_moisture(&h, 0);
        let (_o_wet, ta_wet) = run_layer1_with_moisture(&h, 100);
        let hydro_dry = ta_dry.hydrology.expect("dry run must populate hydrology");
        let hydro_wet = ta_wet.hydrology.expect("wet run must populate hydrology");
        assert_eq!(
            hydro_dry.filled, hydro_wet.filled,
            "lake extent (filled_scaled) must be moisture-independent — only \
             the endorheic/overflow split may vary with moisture_q"
        );
    }
}