//! Geographic feature tag extraction — Layer 1 (D-209).
//!
//! After D8 drainage analysis (D-208), this module extracts the 7
//! `AttractorType` tags from the heightmap + river network. Each attractor has
//! a pixel position and a `strength` (integer 0–100) derived from local terrain
//! quality — computed in floating point, then quantized to an integer at the
//! extraction boundary so all downstream decisions stay integer-deterministic
//! (D-010).
//!
//! **D-209 / D-223 reconciliation:** D-209 reads ocean/lake polygons from
//! `markers.json`, but D-223 reduced markers to a names-only pool — those
//! polygons no longer exist. Coast and lake cells are therefore derived from
//! the heightmap itself: ocean = the largest connected below-sea-level water
//! body; lakes = smaller enclosed below-sea-level bodies.
//!
//! **Determinism (D-010 #4):** all collections iterate in sorted/row-major
//! order; the final attractor list is sorted by `(attractor_type, row, col)`.
//! No `HashMap`/`HashSet` iteration. `strength` (integer 0–100) is not part of
//! that ordering here, so float quantization can't perturb the sort; consumers
//! that rank by strength (e.g. `layer1::attach_feature_names`) do so on the
//! integer value.

use std::collections::{BTreeMap, VecDeque};

use crate::atlas::drainage::DrainageResult;
use crate::atlas::heightmap::BodyHeightmap;
use crate::simulation::generator::AttractorType;

/// 8-neighbor offsets (dr, dc). Columns wrap horizontally (equirectangular
/// globe); rows are bounds-clamped at the poles. Matches `drainage::D8`.
const NB8: [(i32, i32); 8] = [
    (-1, 0),
    (1, 0),
    (0, 1),
    (0, -1),
    (-1, 1),
    (-1, -1),
    (1, 1),
    (1, -1),
];

/// Minimum spacing (cells) between attractors of an areal type, so coastlines /
/// valleys / plains yield a sparse, placement-friendly set rather than one
/// attractor per pixel.
const MIN_SPACING: i32 = 12;

/// Hard cap on attractors per body (keeps the #955 matching tractable).
const MAX_ATTRACTORS: usize = 256;

/// A feature before sub-biome classification: position, type, strength.
/// `layer1` enriches these into `GeographicAttractor`s.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct RawAttractor {
    pub row: u16,
    pub col: u16,
    pub attractor_type: AttractorType,
    /// Strength on a 0–100 integer scale (100 = strongest). Quantized here from
    /// the f32 flow-accumulation ratio — the single f32→integer boundary, after
    /// which every ranking/scoring decision is integer (D-010, #955).
    pub strength: i32,
}

/// Precomputed per-cell terrain fields, shared by feature extraction (D-209)
/// and sub-biome classification (D-210) so neither recomputes them.
#[derive(Debug, Clone)]
pub struct TerrainAnalysis {
    pub w: usize,
    pub h: usize,
    /// `elev < sea_level` (any submerged cell).
    pub ocean_mask: Vec<bool>,
    /// Submerged cells not part of the largest water body (enclosed lakes/seas).
    pub lake_mask: Vec<bool>,
    /// Chebyshev distance (cells) to the nearest ocean cell or river mouth,
    /// capped at `WATER_DIST_CAP`. Moisture proxy for habitability/sub-biome.
    pub water_dist: Vec<u16>,
    /// Local slope proxy in degrees: `atan(max |Δelev| over 8 neighbors)`.
    /// Elevation is normalized [0,1]; this is a relative steepness measure.
    pub slope_deg: Vec<f32>,
    /// Elevation percentile [0,1] among land cells (ocean cells = 0.0).
    pub elev_pct: Vec<f32>,
    /// Settled-equilibrium hydrology sourcing (T-1184, D-227 amendment (4) /
    /// D-255(f) seed-chaining mechanism B). `None` when hydrology hasn't been
    /// solved for this analysis (e.g. every pre-T-1184 call site still using
    /// bare [`TerrainAnalysis::analyze`], and every unit test that constructs
    /// a `TerrainAnalysis` directly without going through the hydrology-aware
    /// entry point) — callers MUST treat `None` as "fall through to the
    /// `ocean_fraction_q` heuristic", never as an error. `Some` when
    /// [`TerrainAnalysis::with_hydrology`] populated it from a real
    /// [`crate::atlas::hydrology_equilibrium::HydrologyResult`].
    pub hydrology: Option<HydrologySample>,
}

/// The two continuous working-grid fields `derive_morphology_zone`'s lake
/// sourcing bilinearly samples (T-1184) — never a discrete basin-membership
/// lookup (that gives blocky, non-refining lake edges, the exact D-166
/// magnified-composite artifact this design avoids; see D-227 amendment (4)
/// / D-255(f) mechanism B). Both fields are row-major, `w × h`, in the SAME
/// `[0.0, 1.0]` normalized domain the raw heightmap and `sea_level` already
/// share — so a bilinear sample of one is directly comparable to a bilinear
/// sample of the other, no rescaling at the call site.
///
/// **Size + clone cost (PR #200 review, Hoshe finding 1):** two `Vec<f32>` at
/// the real 512×256 working grid = ~1.05 MB/entry, added on top of
/// `TerrainAnalysis`'s pre-existing ~1.57 MB of dense fields (~2.62 MB total,
/// ×1.67 growth, not quite a doubling) — see the corrected sizing comment on
/// `GenWorkItem::DeriveWindow` (`gen_queue.rs`) for the full accounting and
/// the `TerrainAnalysisCache` cache-HIT clone-cost note (every hit
/// deep-copies both these `Vec`s, not just the first miss/insert).
///
/// **`elevation` is a deliberate, provably-necessary redundant copy, not an
/// oversight.** `TerrainAnalysis` has no OTHER field that retains the raw
/// `[0,1]` heightmap: `elev_pct` is a RANK PERCENTILE (`rank(elev[i]) /
/// (land_cell_count - 1)`, `compute_elev_percentile`'s own doc/impl) —
/// mathematically a different quantity from absolute elevation, and NOT
/// safe to compare against `filled` (two cells at different true elevations
/// can share adjacent ranks; ocean cells are forced to `0.0` regardless of
/// their real depth). `HydrologyResult` itself carries no elevation field
/// either (`hydrology_equilibrium.rs`: `basins`, `filled_scaled`,
/// `channel_depth_scaled`, `cliff_edge` — no `original`/`elevation` member).
/// So there is no existing bit-identical grid this field could point at
/// instead — carrying its own copy is the only byte-safe option today.
#[derive(Debug, Clone)]
pub struct HydrologySample {
    /// The original (unfilled) heightmap elevation, `[0.0, 1.0]`. Not stored
    /// anywhere else on `TerrainAnalysis` (`elev_pct` is a land-cell RANK
    /// percentile, a different quantity — see its own doc) — this is the
    /// literal `hm.data` the solver's `original` array was built from,
    /// carried alongside `filled` so both halves of the lake comparison
    /// sample from the identical grid at the identical resolution.
    pub elevation: Vec<f32>,
    /// `HydrologyResult.filled_scaled`, rescaled back from the solver's
    /// `i64`-scaled integer domain to `[0.0, 1.0]` (dividing by the same
    /// `ELEV_SCALE` the solver used to go the other way) — the settled
    /// water-surface height at every working-grid cell (equal to
    /// `elevation` wherever no lake exists).
    pub filled: Vec<f32>,
    /// `HydrologyResult.basin_max_depth_scaled`, rescaled back to `[0.0,
    /// 1.0]` fraction units (T-1188): the MAXIMUM settled depth anywhere in
    /// this cell's basin, broadcast to every cell in that basin, `0.0` for
    /// non-lake cells. Used to normalize `lake_margin_q` per-basin instead
    /// of against a single fixed absolute ceiling — see
    /// `district_profile::lake_from_hydrology_at`'s doc for the full
    /// rationale (the PR #206 eyeball finding that motivated this field).
    pub basin_max_depth: Vec<f32>,
}

const WATER_DIST_CAP: u16 = 255;

#[inline]
fn idx(r: usize, c: usize, w: usize) -> usize {
    r * w + c
}

#[inline]
fn wrap_col(c: i32, w: i32) -> usize {
    c.rem_euclid(w) as usize
}

impl TerrainAnalysis {
    /// Compute all shared terrain fields for a body. O(w·h).
    pub fn analyze(hm: &BodyHeightmap, drainage: &DrainageResult) -> TerrainAnalysis {
        let w = hm.width as usize;
        let h = hm.height as usize;
        let n = w * h;
        let elev = &hm.data;
        let sea = hm.sea_level;

        let ocean_mask: Vec<bool> = (0..n).map(|i| elev[i] < sea).collect();
        let lake_mask = compute_lake_mask(&ocean_mask, w, h);
        let water_dist = compute_water_dist(&ocean_mask, &drainage.river_network.mouths, w, h);
        let slope_deg = compute_slope(elev, w, h);
        let elev_pct = compute_elev_percentile(elev, &ocean_mask, w, h);

        TerrainAnalysis {
            w,
            h,
            ocean_mask,
            lake_mask,
            water_dist,
            slope_deg,
            elev_pct,
            hydrology: None,
        }
    }

    /// Populate the settled-hydrology sourcing fields (T-1184, D-227
    /// amendment (4) / D-255(f) mechanism B) from a solved
    /// [`crate::atlas::hydrology_equilibrium::HydrologyResult`].
    ///
    /// Builder-style (consumes and returns `self`) rather than a constructor
    /// parameter on [`TerrainAnalysis::analyze`] — `analyze` has ~20 call
    /// sites across production code and tests that have no hydrology input
    /// (and, per D-227, don't need one: hydrology sourcing is a lake-specific
    /// refinement, not a precondition for every other terrain field this
    /// struct carries). Keeping `analyze`'s signature untouched means every
    /// existing caller keeps working byte-identically; only the two
    /// production sites that actually solve hydrology
    /// (`layer1::run_layer1`, `gen_queue::TerrainAnalysisCache::get_or_derive`)
    /// opt in by chaining this call.
    ///
    /// Panics if `result`'s grids aren't `self.w * self.h` cells — a
    /// programmer error (mismatched working-grid resolution between the
    /// heightmap this `TerrainAnalysis` was built from and the elevation grid
    /// `solve()` was called on), never a legitimate runtime state.
    pub fn with_hydrology(
        mut self,
        elevation: &[f32],
        result: &crate::atlas::hydrology_equilibrium::HydrologyResult,
    ) -> TerrainAnalysis {
        let n = self.w * self.h;
        assert_eq!(
            elevation.len(),
            n,
            "with_hydrology: elevation grid size does not match TerrainAnalysis dims"
        );
        assert_eq!(
            result.filled_scaled.len(),
            n,
            "with_hydrology: HydrologyResult grid size does not match TerrainAnalysis dims"
        );
        let filled: Vec<f32> = result
            .filled_scaled
            .iter()
            .map(|&s| crate::atlas::hydrology_equilibrium::scaled_to_fraction(s))
            .collect();
        let basin_max_depth: Vec<f32> = result
            .basin_max_depth_scaled
            .iter()
            .map(|&s| crate::atlas::hydrology_equilibrium::scaled_to_fraction(s))
            .collect();
        self.hydrology = Some(HydrologySample {
            elevation: elevation.to_vec(),
            filled,
            basin_max_depth,
        });
        self
    }

    #[inline]
    pub fn is_ocean(&self, r: usize, c: usize) -> bool {
        self.ocean_mask[idx(r, c, self.w)]
    }

    /// Compass bearing toward the nearest water from cell `(r, c)`, quantized to
    /// 8 octants (0=N, 45=NE … 315=NW); `360` = "no water in range" (#957, D-234).
    ///
    /// Reads the `water_dist` field's local gradient — the 8-neighbour with the
    /// smallest distance-to-water points toward water. Integer-only (no `atan2`)
    /// for D-010 determinism. Returns `360` when the cell is itself water or no
    /// neighbour is closer to water (flat/inland).
    pub fn water_bearing(&self, r: usize, c: usize) -> u16 {
        let here = self.water_dist[idx(r, c, self.w)];
        if here == 0 || here >= WATER_DIST_CAP {
            return NO_WATER_BEARING; // on water, or no water within range
        }
        let mut best = here;
        let mut bdir = (0i32, 0i32);
        for &(dr, dc) in &NB8 {
            let nr = r as i32 + dr;
            if nr < 0 || nr >= self.h as i32 {
                continue;
            }
            let nc = wrap_col(c as i32 + dc, self.w as i32);
            let nd = self.water_dist[idx(nr as usize, nc, self.w)];
            if nd < best {
                best = nd;
                bdir = (dr, dc);
            }
        }
        if bdir == (0, 0) {
            NO_WATER_BEARING
        } else {
            octant_bearing(bdir.0, bdir.1)
        }
    }
}

/// Sentinel for [`TerrainAnalysis::water_bearing`] meaning "no water direction".
pub const NO_WATER_BEARING: u16 = 360;

/// Quantize a (Δrow, Δcol) step to a compass octant bearing (0=N … 315=NW).
/// `Δrow < 0` is north (rows increase downward). Integer-only (D-010).
fn octant_bearing(drow: i32, dcol: i32) -> u16 {
    let (ar, ac) = (drow.abs(), dcol.abs());
    let north = drow < 0;
    let east = dcol > 0;
    if ar >= ac * 2 {
        if north {
            0
        } else {
            180
        }
    } else if ac >= ar * 2 {
        if east {
            90
        } else {
            270
        }
    } else {
        match (north, east) {
            (true, true) => 45,
            (true, false) => 315,
            (false, true) => 135,
            (false, false) => 225,
        }
    }
}

/// Largest connected below-sea-level component = ocean; all others = lakes.
/// Deterministic: BFS seeds scanned row-major; ties broken by lowest cell index.
fn compute_lake_mask(ocean_mask: &[bool], w: usize, h: usize) -> Vec<bool> {
    let n = w * h;
    let mut comp = vec![-1i32; n];
    let mut comp_sizes: Vec<usize> = Vec::new();
    let mut next_comp = 0i32;

    for start in 0..n {
        if !ocean_mask[start] || comp[start] >= 0 {
            continue;
        }
        // Flood fill this component (row-major BFS = deterministic).
        let mut size = 0usize;
        let mut q = VecDeque::new();
        comp[start] = next_comp;
        q.push_back(start);
        while let Some(cur) = q.pop_front() {
            size += 1;
            let (r, c) = (cur / w, cur % w);
            for &(dr, dc) in &NB8 {
                let nr = r as i32 + dr;
                if nr < 0 || nr >= h as i32 {
                    continue;
                }
                let nc = wrap_col(c as i32 + dc, w as i32);
                let ni = idx(nr as usize, nc, w);
                if ocean_mask[ni] && comp[ni] < 0 {
                    comp[ni] = next_comp;
                    q.push_back(ni);
                }
            }
        }
        comp_sizes.push(size);
        next_comp += 1;
    }

    if comp_sizes.is_empty() {
        return vec![false; n]; // no water at all
    }

    // Largest component (tie → lowest comp id, which is the earliest row-major).
    let mut ocean_comp = 0i32;
    let mut best = 0usize;
    for (cid, &sz) in comp_sizes.iter().enumerate() {
        if sz > best {
            best = sz;
            ocean_comp = cid as i32;
        }
    }

    // Lakes = submerged cells in any non-ocean component.
    (0..n)
        .map(|i| comp[i] >= 0 && comp[i] != ocean_comp)
        .collect()
}

/// Multi-source BFS Chebyshev distance to nearest ocean cell or river mouth.
fn compute_water_dist(ocean_mask: &[bool], mouths: &[(u16, u16)], w: usize, h: usize) -> Vec<u16> {
    let n = w * h;
    let mut dist = vec![WATER_DIST_CAP; n];
    let mut q = VecDeque::new();
    // Seeds in row-major order for determinism.
    for i in 0..n {
        if ocean_mask[i] {
            dist[i] = 0;
            q.push_back(i);
        }
    }
    for &(mr, mc) in mouths {
        let i = idx(mr as usize, mc as usize, w);
        if dist[i] != 0 {
            dist[i] = 0;
            q.push_back(i);
        }
    }
    while let Some(cur) = q.pop_front() {
        let d = dist[cur];
        if d >= WATER_DIST_CAP {
            continue;
        }
        let (r, c) = (cur / w, cur % w);
        for &(dr, dc) in &NB8 {
            let nr = r as i32 + dr;
            if nr < 0 || nr >= h as i32 {
                continue;
            }
            let nc = wrap_col(c as i32 + dc, w as i32);
            let ni = idx(nr as usize, nc, w);
            if dist[ni] > d + 1 {
                dist[ni] = d + 1;
                q.push_back(ni);
            }
        }
    }
    dist
}

/// Local slope proxy: `atan(max |Δelev| to 8 neighbors)` in degrees.
fn compute_slope(elev: &[f32], w: usize, h: usize) -> Vec<f32> {
    let n = w * h;
    let mut slope = vec![0.0f32; n];
    for r in 0..h {
        for c in 0..w {
            let i = idx(r, c, w);
            let e = elev[i];
            let mut max_grad = 0.0f32;
            for &(dr, dc) in &NB8 {
                let nr = r as i32 + dr;
                if nr < 0 || nr >= h as i32 {
                    continue;
                }
                let nc = wrap_col(c as i32 + dc, w as i32);
                let g = (e - elev[idx(nr as usize, nc, w)]).abs();
                if g > max_grad {
                    max_grad = g;
                }
            }
            slope[i] = max_grad.atan().to_degrees();
        }
    }
    slope
}

/// Elevation percentile [0,1] among land cells; ocean cells get 0.0.
fn compute_elev_percentile(elev: &[f32], ocean_mask: &[bool], w: usize, h: usize) -> Vec<f32> {
    let n = w * h;
    // (scaled_elev, idx) for land cells; integer key for deterministic sort.
    let mut land: Vec<(i64, usize)> = (0..n)
        .filter(|&i| !ocean_mask[i])
        .map(|i| ((elev[i] as f64 * 1_000_000.0) as i64, i))
        .collect();
    land.sort_unstable_by(|a, b| a.0.cmp(&b.0).then(a.1.cmp(&b.1)));
    let mut pct = vec![0.0f32; n];
    let m = land.len();
    if m <= 1 {
        for &(_, i) in &land {
            pct[i] = 0.5;
        }
        return pct;
    }
    for (rank, &(_, i)) in land.iter().enumerate() {
        pct[i] = rank as f32 / (m - 1) as f32;
    }
    pct
}

/// Habitability score [0,1] from elevation band, flatness, and moisture.
/// Used by `ValleyFloor`/`PlainCenter` strengths and the D-210 classifier.
pub fn habitability(elev_pct: f32, slope_deg: f32, water_dist: u16) -> f32 {
    let elev_score = (1.0 - (elev_pct - 0.35).abs() / 0.65).clamp(0.0, 1.0);
    let flat_score = (1.0 - slope_deg / 15.0).clamp(0.0, 1.0);
    let water_score = (1.0 - water_dist as f32 / 40.0).clamp(0.0, 1.0);
    (0.4 * elev_score + 0.3 * flat_score + 0.3 * water_score).clamp(0.0, 1.0)
}

/// Extract the 7 D-209 attractor tags. Returns raw attractors (no sub-biome),
/// sorted by `(attractor_type, row, col)` for determinism. Higher-priority
/// types claim their cells first so a cell is tagged at most once.
pub fn extract_attractors(
    hm: &BodyHeightmap,
    drainage: &DrainageResult,
    ta: &TerrainAnalysis,
) -> Vec<RawAttractor> {
    let w = hm.width as usize;
    let h = hm.height as usize;
    let elev = &hm.data;
    let accum = &drainage.flow_accumulation;
    let max_accum = drainage.max_accumulation.max(1) as f32;

    let mut claimed = vec![false; w * h];
    let mut out: Vec<RawAttractor> = Vec::new();

    // O(1) river-cell membership (avoids a binary_search per valley candidate).
    let mut river_mask = vec![false; w * h];
    for &(r, c) in &drainage.river_network.river_cells {
        river_mask[idx(r as usize, c as usize, w)] = true;
    }

    let claim = |out: &mut Vec<RawAttractor>,
                 claimed: &mut [bool],
                 r: usize,
                 c: usize,
                 at: AttractorType,
                 strength: f32| {
        let i = idx(r, c, w);
        if claimed[i] {
            return;
        }
        claimed[i] = true;
        out.push(RawAttractor {
            row: r as u16,
            col: c as u16,
            attractor_type: at,
            // The single f32→integer boundary: quantize the 0.0–1.0 ratio to 0–100.
            strength: (strength.clamp(0.0, 1.0) * 100.0).round() as i32,
        });
    };

    // 1. RiverMouth — strength = accum / max_accum.
    for &(r, c) in &drainage.river_network.mouths {
        let i = idx(r as usize, c as usize, w);
        let s = accum[i] as f32 / max_accum;
        claim(
            &mut out,
            &mut claimed,
            r as usize,
            c as usize,
            AttractorType::RiverMouth,
            s,
        );
    }

    // 2. RiverCrossing — confluences, strength = accum / max_accum * 0.7.
    for &(r, c) in &drainage.river_network.confluences {
        let i = idx(r as usize, c as usize, w);
        let s = accum[i] as f32 / max_accum * 0.7;
        claim(
            &mut out,
            &mut claimed,
            r as usize,
            c as usize,
            AttractorType::RiverCrossing,
            s,
        );
    }

    // 3. CoastalAccess — land within 3 cells of ocean, thinned by spacing.
    //    strength = 0.6 + coast-density bonus (capped).
    let mut coastal: Vec<(usize, usize, f32)> = Vec::new();
    for r in 0..h {
        for c in 0..w {
            let i = idx(r, c, w);
            // water_dist (precomputed) is a cheap pre-filter: only cells within
            // 3 of water can be coastal, so skip the 49-cell scan for inland.
            if ta.ocean_mask[i] || claimed[i] || ta.water_dist[i] > 3 {
                continue;
            }
            let near = ocean_cells_within(ta, r, c, 3);
            if near > 0 {
                let bonus = (near as f32 / 24.0).min(0.3);
                coastal.push((r, c, 0.6 + bonus));
            }
        }
    }
    for (r, c, s) in thin_by_spacing(coastal, &claimed, w) {
        claim(
            &mut out,
            &mut claimed,
            r,
            c,
            AttractorType::CoastalAccess,
            s,
        );
    }

    // 4. ValleyFloor — gentle slope, mid elevation, positive habitability.
    let mut valleys: Vec<(usize, usize, f32)> = Vec::new();
    for r in 0..h {
        for c in 0..w {
            let i = idx(r, c, w);
            if ta.ocean_mask[i] || claimed[i] {
                continue;
            }
            if river_mask[i] || ta.slope_deg[i] >= 5.0 {
                continue;
            }
            if ta.elev_pct[i] < 0.10 || ta.elev_pct[i] > 0.60 {
                continue;
            }
            let hab = habitability(ta.elev_pct[i], ta.slope_deg[i], ta.water_dist[i]);
            if hab > 0.0 {
                valleys.push((r, c, hab));
            }
        }
    }
    for (r, c, s) in thin_by_spacing(valleys, &claimed, w) {
        claim(&mut out, &mut claimed, r, c, AttractorType::ValleyFloor, s);
    }

    // 5. PassEntrance — morphological saddles in higher terrain.
    //    strength = 1 - elev_pct (lower passes score higher).
    let mut passes: Vec<(usize, usize, f32)> = Vec::new();
    for r in 1..h.saturating_sub(1) {
        for c in 0..w {
            let i = idx(r, c, w);
            if ta.ocean_mask[i] || claimed[i] || ta.elev_pct[i] < 0.5 {
                continue;
            }
            if is_saddle(elev, r, c, w, h) {
                passes.push((r, c, 1.0 - ta.elev_pct[i]));
            }
        }
    }
    for (r, c, s) in thin_by_spacing(passes, &claimed, w) {
        claim(&mut out, &mut claimed, r, c, AttractorType::PassEntrance, s);
    }

    // 6. LakeShore — land adjacent to an enclosed lake. strength = 0.5.
    let mut shores: Vec<(usize, usize, f32)> = Vec::new();
    for r in 0..h {
        for c in 0..w {
            let i = idx(r, c, w);
            if ta.ocean_mask[i] || ta.lake_mask[i] || claimed[i] {
                continue;
            }
            if adjacent_to_lake(ta, r, c) {
                shores.push((r, c, 0.5));
            }
        }
    }
    for (r, c, s) in thin_by_spacing(shores, &claimed, w) {
        claim(&mut out, &mut claimed, r, c, AttractorType::LakeShore, s);
    }

    // 7. PlainCenter — very flat, away from everything. strength = hab * 0.4.
    let mut plains: Vec<(usize, usize, f32)> = Vec::new();
    for r in 0..h {
        for c in 0..w {
            let i = idx(r, c, w);
            if ta.ocean_mask[i] || claimed[i] || ta.slope_deg[i] >= 2.0 {
                continue;
            }
            let hab = habitability(ta.elev_pct[i], ta.slope_deg[i], ta.water_dist[i]);
            plains.push((r, c, hab * 0.4));
        }
    }
    for (r, c, s) in thin_by_spacing(plains, &claimed, w) {
        claim(&mut out, &mut claimed, r, c, AttractorType::PlainCenter, s);
    }

    // Cap to MAX_ATTRACTORS while preserving type diversity. D-209 calls
    // RiverMouth "always high-value", but its *normalized* strength
    // (accum / max_accum) is tiny for all but the largest river, so a pure
    // global-strength cap lets abundant ValleyFloor/CoastalAccess crowd every
    // RiverMouth out. Instead: group by type, sort each group strongest-first,
    // then round-robin across types so every present type keeps representation.
    // Deterministic (BTreeMap type order, integer strength key, fixed rotation).
    if out.len() > MAX_ATTRACTORS {
        let mut by_type: std::collections::BTreeMap<u8, Vec<RawAttractor>> =
            std::collections::BTreeMap::new();
        for a in out.drain(..) {
            by_type.entry(a.attractor_type as u8).or_default().push(a);
        }
        for group in by_type.values_mut() {
            group.sort_by(|a, b| {
                // strength is integer now — rank by it directly (descending).
                b.strength
                    .cmp(&a.strength)
                    .then(a.row.cmp(&b.row))
                    .then(a.col.cmp(&b.col))
            });
        }
        let mut kept: Vec<RawAttractor> = Vec::with_capacity(MAX_ATTRACTORS);
        let mut depth = 0usize;
        'fill: loop {
            let mut progressed = false;
            for group in by_type.values() {
                if let Some(a) = group.get(depth) {
                    kept.push(*a);
                    progressed = true;
                    if kept.len() >= MAX_ATTRACTORS {
                        break 'fill;
                    }
                }
            }
            if !progressed {
                break;
            }
            depth += 1;
        }
        out = kept;
    }
    out.sort_by(|a, b| {
        (a.attractor_type as u8, a.row, a.col).cmp(&(b.attractor_type as u8, b.row, b.col))
    });
    out
}

fn ocean_cells_within(ta: &TerrainAnalysis, r: usize, c: usize, radius: i32) -> usize {
    let mut count = 0;
    for dr in -radius..=radius {
        let nr = r as i32 + dr;
        if nr < 0 || nr >= ta.h as i32 {
            continue;
        }
        for dc in -radius..=radius {
            let nc = wrap_col(c as i32 + dc, ta.w as i32);
            if ta.ocean_mask[idx(nr as usize, nc, ta.w)] {
                count += 1;
            }
        }
    }
    count
}

fn adjacent_to_lake(ta: &TerrainAnalysis, r: usize, c: usize) -> bool {
    for &(dr, dc) in &NB8 {
        let nr = r as i32 + dr;
        if nr < 0 || nr >= ta.h as i32 {
            continue;
        }
        let nc = wrap_col(c as i32 + dc, ta.w as i32);
        if ta.lake_mask[idx(nr as usize, nc, ta.w)] {
            return true;
        }
    }
    false
}

/// Morphological saddle: walking the 8-neighbor ring, the sign of
/// `(neighbor - cell)` alternates at least 4 times (≥2 higher sectors
/// separated by ≥2 lower sectors).
fn is_saddle(elev: &[f32], r: usize, c: usize, w: usize, h: usize) -> bool {
    // Ring order (clockwise) so transitions are meaningful.
    const RING: [(i32, i32); 8] = [
        (-1, 0),
        (-1, 1),
        (0, 1),
        (1, 1),
        (1, 0),
        (1, -1),
        (0, -1),
        (-1, -1),
    ];
    let e = elev[idx(r, c, w)];
    let mut signs = [0i8; 8];
    for (k, &(dr, dc)) in RING.iter().enumerate() {
        let nr = r as i32 + dr;
        if nr < 0 || nr >= h as i32 {
            return false; // poles can't be saddles in this scheme
        }
        let nc = wrap_col(c as i32 + dc, w as i32);
        signs[k] = if elev[idx(nr as usize, nc, w)] > e {
            1
        } else {
            -1
        };
    }
    let mut transitions = 0;
    for k in 0..8 {
        if signs[k] != signs[(k + 1) % 8] {
            transitions += 1;
        }
    }
    transitions >= 4
}

/// Greedy spatial thinning: sort candidates by descending strength (ties by
/// row, col), keep one per `MIN_SPACING` Chebyshev neighborhood.
///
/// Uses a bucket grid (cell size = `MIN_SPACING`) so each candidate only checks
/// the 3×3 neighboring buckets — O(k) amortized rather than O(k²). The kept
/// order is fully determined by the sorted candidate iteration; the bucket map
/// is `BTreeMap` (the project bans `HashMap` for determinism) and is lookup-only
/// regardless.
fn thin_by_spacing(
    mut cands: Vec<(usize, usize, f32)>,
    _claimed: &[bool],
    _w: usize,
) -> Vec<(usize, usize, f32)> {
    // Deterministic order: strength desc, then row, col asc.
    cands.sort_by(|a, b| {
        let sa = (a.2 * 1e6) as i64;
        let sb = (b.2 * 1e6) as i64;
        sb.cmp(&sa).then(a.0.cmp(&b.0)).then(a.1.cmp(&b.1))
    });
    let sp = MIN_SPACING.max(1) as usize;
    let mut buckets: BTreeMap<(usize, usize), Vec<(usize, usize)>> = BTreeMap::new();
    let mut kept: Vec<(usize, usize, f32)> = Vec::new();
    for (r, c, s) in cands {
        let (br, bc) = (r / sp, c / sp);
        let mut ok = true;
        'scan: for nbr in br.saturating_sub(1)..=br + 1 {
            for nbc in bc.saturating_sub(1)..=bc + 1 {
                if let Some(pts) = buckets.get(&(nbr, nbc)) {
                    for &(kr, kc) in pts {
                        let dr = (kr as i32 - r as i32).abs();
                        let dc = (kc as i32 - c as i32).abs();
                        if dr.max(dc) < MIN_SPACING {
                            ok = false;
                            break 'scan;
                        }
                    }
                }
            }
        }
        if ok {
            buckets.entry((br, bc)).or_default().push((r, c));
            kept.push((r, c, s));
        }
    }
    kept
}

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

    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(data: Vec<f32>, w: u32, h: u32, sea: f32) -> BodyHeightmap {
        BodyHeightmap {
            body_id: "T".into(),
            width: w,
            height: h,
            data,
            sea_level: sea,
        }
    }

    #[test]
    fn deterministic_extraction() {
        let h = hm(slope_grid(64, 32), 64, 32, 0.3);
        let dr = drainage::analyze(&h.data, 64, 32, 0.3);
        let ta = TerrainAnalysis::analyze(&h, &dr);
        let a1 = extract_attractors(&h, &dr, &ta);
        let a2 = extract_attractors(&h, &dr, &ta);
        assert_eq!(a1, a2, "attractor extraction must be deterministic");
    }

    #[test]
    fn attractors_sorted_and_bounded() {
        let h = hm(slope_grid(128, 64), 128, 64, 0.3);
        let dr = drainage::analyze(&h.data, 128, 64, 0.3);
        let ta = TerrainAnalysis::analyze(&h, &dr);
        let a = extract_attractors(&h, &dr, &ta);
        assert!(a.len() <= MAX_ATTRACTORS);
        // Sorted by (type as u8, row, col).
        for win in a.windows(2) {
            let ka = (win[0].attractor_type as u8, win[0].row, win[0].col);
            let kb = (win[1].attractor_type as u8, win[1].row, win[1].col);
            assert!(ka <= kb, "attractors must be sorted");
        }
    }

    #[test]
    fn percentile_in_range() {
        let h = hm(slope_grid(32, 16), 32, 16, 0.3);
        let dr = drainage::analyze(&h.data, 32, 16, 0.3);
        let ta = TerrainAnalysis::analyze(&h, &dr);
        assert!(ta.elev_pct.iter().all(|&p| (0.0..=1.0).contains(&p)));
        assert_eq!(ta.slope_deg.len(), 32 * 16);
    }

    /// Multi-octave sine terrain (continents + many small coastal streams) —
    /// produces > MAX_ATTRACTORS candidates with plenty of river mouths.
    fn sine_grid(w: u32, h: u32) -> Vec<f32> {
        use std::f32::consts::{PI, TAU};
        (0..(w * h))
            .map(|i| {
                let r = (i / w) as f32;
                let c = (i % w) as f32;
                let x = c / w as f32 * TAU;
                let y = r / h as f32 * PI;
                (0.5 + 0.25 * (x * 3.0).sin() * (y * 2.0).sin()
                    + 0.15 * (x * 7.0).cos() * (y * 5.0).sin()
                    + 0.08 * (x * 13.0).sin() * (y * 11.0).cos()
                    + 0.05 * (x * 23.0).cos() * (y * 19.0).sin())
                .clamp(0.0, 1.0)
            })
            .collect()
    }

    #[test]
    fn river_mouths_survive_cap() {
        // D-209 + the type-aware cap: even though RiverMouth normalized strength
        // is tiny, a body full of mouths must still keep RiverMouth attractors
        // (a global-strength cap would drop all of them — the bug Hoshe caught).
        let h = hm(sine_grid(512, 256), 512, 256, 0.40);
        let dr = drainage::analyze(&h.data, 512, 256, 0.40);
        assert!(
            !dr.river_network.mouths.is_empty(),
            "fixture must have mouths"
        );
        let ta = TerrainAnalysis::analyze(&h, &dr);
        let a = extract_attractors(&h, &dr, &ta);
        assert!(a.len() <= MAX_ATTRACTORS);
        assert!(
            a.iter()
                .any(|x| x.attractor_type == AttractorType::RiverMouth),
            "RiverMouth attractors must survive the cap when mouths exist"
        );
    }

    // -----------------------------------------------------------------------
    // Per-type attractor reachability (T-964): crafted heightmaps that
    // guarantee at least one attractor of the named type, isolating each
    // extraction branch instead of relying on the slope/sine fixtures above
    // (which reliably exercise RiverMouth/CoastalAccess/ValleyFloor, but never
    // guarantee LakeShore/PassEntrance/PlainCenter/RiverCrossing).
    // -----------------------------------------------------------------------

    #[test]
    fn lake_shore_reachable_via_enclosed_depression() {
        // Two separate below-sea-level components: a wide strip along the
        // west edge (the largest — becomes ocean) and a small isolated pit
        // elsewhere (smaller — becomes an enclosed lake, D-209/compute_lake_mask).
        // Land cells 8-adjacent to the pit must classify LakeShore.
        let (w, h) = (32usize, 16usize);
        let mut data = vec![0.6f32; w * h];
        for r in 0..h {
            for c in 0..4 {
                data[r * w + c] = 0.1; // wide ocean strip
            }
        }
        for r in 6..8 {
            for c in 16..18 {
                data[r * w + c] = 0.1; // small isolated pit, far from the ocean
            }
        }
        let heightmap = hm(data, w as u32, h as u32, 0.3);
        let dr = drainage::analyze(&heightmap.data, w as u32, h as u32, 0.3);
        let ta = TerrainAnalysis::analyze(&heightmap, &dr);
        assert!(
            ta.lake_mask.iter().any(|&x| x),
            "fixture sanity: the isolated pit must register as a lake, not ocean"
        );
        let a = extract_attractors(&heightmap, &dr, &ta);
        assert!(
            a.iter()
                .any(|x| x.attractor_type == AttractorType::LakeShore),
            "land adjacent to an enclosed lake must classify LakeShore"
        );
    }

    #[test]
    fn pass_entrance_reachable_via_morphological_saddle() {
        // Classic saddle: the 8-ring around the center alternates high/low
        // going clockwise (N,NE,E,SE,S,SW,W,NW), giving 8 sign transitions
        // (is_saddle requires >= 4). Whole grid is high-elevation land so the
        // saddle's elev_pct clears the >= 0.5 PassEntrance gate.
        let (w, h) = (32usize, 16usize);
        let (cr, cc) = (h / 2, w / 2);
        let mut data = vec![0.7f32; w * h];
        data[cr * w + cc] = 0.75; // the saddle point itself
        const RING: [(i32, i32); 8] = [
            (-1, 0),
            (-1, 1),
            (0, 1),
            (1, 1),
            (1, 0),
            (1, -1),
            (0, -1),
            (-1, -1),
        ];
        let ring_vals = [0.95, 0.55, 0.95, 0.55, 0.95, 0.55, 0.95, 0.55];
        for (k, &(dr_off, dc_off)) in RING.iter().enumerate() {
            let rr = (cr as i32 + dr_off) as usize;
            let cc_ = (cc as i32 + dc_off) as usize;
            data[rr * w + cc_] = ring_vals[k];
        }
        let heightmap = hm(data, w as u32, h as u32, 0.0);
        let dr = drainage::analyze(&heightmap.data, w as u32, h as u32, 0.0);
        let ta = TerrainAnalysis::analyze(&heightmap, &dr);
        assert!(
            ta.elev_pct[cr * w + cc] >= 0.5,
            "fixture sanity: saddle point must clear the PassEntrance elev_pct gate"
        );
        let a = extract_attractors(&heightmap, &dr, &ta);
        assert!(
            a.iter()
                .any(|x| x.attractor_type == AttractorType::PassEntrance),
            "a genuine morphological saddle at high elevation must classify PassEntrance"
        );
    }

    #[test]
    fn plain_center_reachable_via_flat_uniform_terrain() {
        // A uniformly flat, non-ocean grid: slope_deg is 0 everywhere (well
        // under the < 2.0 PlainCenter gate), so at least one cell survives
        // thin_by_spacing as PlainCenter even where ValleyFloor also
        // competes for the uniform elev_pct=0.5 rank tie.
        let (w, h) = (32usize, 16usize);
        let data = vec![0.9f32; w * h];
        let heightmap = hm(data, w as u32, h as u32, 0.0);
        let dr = drainage::analyze(&heightmap.data, w as u32, h as u32, 0.0);
        let ta = TerrainAnalysis::analyze(&heightmap, &dr);
        let a = extract_attractors(&heightmap, &dr, &ta);
        assert!(
            a.iter()
                .any(|x| x.attractor_type == AttractorType::PlainCenter),
            "flat, non-ocean terrain must produce at least one PlainCenter attractor"
        );
    }

    #[test]
    fn river_crossing_reachable_via_confluence() {
        // Two V-shaped tributary valleys (west + east branches) converge into
        // a single trunk valley at (confluence_row, confluence_col) — the
        // trunk cell has 2+ river-cell inflows, so `drainage::analyze` must
        // report it as a confluence (drainage.rs's own D8 confluence rule),
        // and extract_attractors must tag it RiverCrossing.
        let (w, h) = (64usize, 64usize);
        let confluence_col = (w / 2) as f32;
        let confluence_row = (h / 2) as f32;
        let data: Vec<f32> = (0..(w * h))
            .map(|i| {
                let r = (i / w) as f32;
                let c = (i % w) as f32;
                if r <= confluence_row {
                    // Upstream: two separate branches either side of the
                    // confluence column, each sloping down toward it.
                    let branch_center = if c < confluence_col {
                        confluence_col * 0.5
                    } else {
                        confluence_col * 1.5
                    };
                    let lateral = (c - branch_center).abs() / w as f32;
                    let downstream = (confluence_row - r) / h as f32;
                    (0.3 + lateral * 1.5 - downstream * 0.4).clamp(0.0, 1.0)
                } else {
                    // Downstream: single widening trunk valley.
                    let lateral = (c - confluence_col).abs() / w as f32;
                    let downstream = (r - confluence_row) / h as f32;
                    (0.3 + lateral * 1.5 - downstream * 0.6).clamp(0.0, 1.0)
                }
            })
            .collect();
        let heightmap = hm(data, w as u32, h as u32, 0.0);
        let dr = drainage::analyze(&heightmap.data, w as u32, h as u32, 0.0);
        assert!(
            !dr.river_network.confluences.is_empty(),
            "fixture sanity: the converging-tributary fixture must produce a confluence"
        );
        let ta = TerrainAnalysis::analyze(&heightmap, &dr);
        let a = extract_attractors(&heightmap, &dr, &ta);
        assert!(
            a.iter()
                .any(|x| x.attractor_type == AttractorType::RiverCrossing),
            "a genuine D8 confluence must classify RiverCrossing"
        );
    }

    // -----------------------------------------------------------------------
    // thin_by_spacing behavior (T-964): spacing collisions + equirectangular
    // column wrap.
    // -----------------------------------------------------------------------

    #[test]
    fn thin_by_spacing_drops_close_candidates_keeps_strongest() {
        // Three candidates within MIN_SPACING (12) of each other: only the
        // strongest should survive; a fourth, far-away candidate is
        // independent and must survive alongside it.
        let claimed = vec![false; 64 * 64];
        let cands = vec![
            (10usize, 10usize, 0.5f32),
            (10usize, 15usize, 0.9f32), // strongest, within spacing of the other two
            (15usize, 10usize, 0.3f32),
            (50usize, 50usize, 0.4f32), // far away — independent, must survive
        ];
        let kept = thin_by_spacing(cands, &claimed, 64);
        assert_eq!(
            kept.len(),
            2,
            "expected exactly 2 survivors (the strongest of the clustered trio + the \
             far-away independent point), got {kept:?}"
        );
        assert!(
            kept.contains(&(10, 15, 0.9)),
            "the strongest candidate in the cluster must survive: {kept:?}"
        );
        assert!(
            kept.contains(&(50, 50, 0.4)),
            "the far-away independent candidate must survive: {kept:?}"
        );
    }

    #[test]
    fn thin_by_spacing_respects_exact_spacing_boundary() {
        // Chebyshev distance exactly MIN_SPACING (12) apart must NOT collide
        // (the check is `dr.max(dc) < MIN_SPACING`, a strict less-than) — both
        // survive. One cell short of that (11) must collide — only the
        // stronger survives.
        let claimed = vec![false; 64 * 64];

        let at_boundary = vec![(0usize, 0usize, 0.5f32), (12usize, 0usize, 0.5f32)];
        let kept_boundary = thin_by_spacing(at_boundary, &claimed, 64);
        assert_eq!(
            kept_boundary.len(),
            2,
            "cells exactly MIN_SPACING apart must both survive (strict <): {kept_boundary:?}"
        );

        let inside_spacing = vec![(0usize, 0usize, 0.5f32), (11usize, 0usize, 0.9f32)];
        let kept_inside = thin_by_spacing(inside_spacing, &claimed, 64);
        assert_eq!(
            kept_inside.len(),
            1,
            "cells 1 short of MIN_SPACING must collide, keeping only the stronger: \
             {kept_inside:?}"
        );
        assert_eq!(kept_inside[0], (11, 0, 0.9));
    }

    #[test]
    fn thin_by_spacing_is_not_wrap_aware_pins_current_behavior() {
        // thin_by_spacing itself is a pure Chebyshev-distance thinner over
        // (row, col) pairs — it has NO knowledge of the equirectangular
        // column wrap (unlike NB8-based neighbor walks elsewhere in this
        // file, which wrap explicitly via `wrap_col`). Two candidates at
        // opposite ends of a wide grid (col 0 and col w-1) are geographically
        // adjacent on the globe but numerically far apart in (row, col)
        // space, so thin_by_spacing does NOT treat them as colliding — both
        // survive. This is NOT proof that wrap support exists or is verified
        // — it pins the OPPOSITE: the current lack-of-wrap-awareness, so a
        // future change to make thinning wrap-aware is a deliberate, visible
        // decision (this test would need to be rewritten), not a silent
        // behavior drift.
        let w = 64usize;
        let claimed = vec![false; w * 64];
        let cands = vec![(5usize, 0usize, 0.5f32), (5usize, w - 1, 0.6f32)];
        let kept = thin_by_spacing(cands, &claimed, w);
        assert_eq!(
            kept.len(),
            2,
            "column-wrap-adjacent candidates are numerically far apart in (row, col) \
             space — thin_by_spacing must not collide them: {kept:?}"
        );
    }
}