facett-core 0.1.15

//! **The cull index** — Ragnar's S+tree, not a Hilbert R-tree.
//!
//! Gemini's spec proposed a packed Hilbert R-tree for the CPU-side slicing.
//! Rejected: the constellation already owns a **measured** better structure. The
//! S(+)-tree in `znippy_zoomies::stree32` is Ragnar Groot Koerkamp's cache-optimal
//! implicit tree — one 64-byte cache line per node, a branchless AVX2
//! `count_lt`, and a software-pipelined prefetch batch traversal. Its published
//! headline is **40× faster than binary search**, and it is already a dependency
//! of `facett-core` (the `gatling` edge), so adopting it costs zero new deps.
//!
//! ## How a *range* query becomes a *key* query
//! The S+tree answers "where does key `q` sit in this sorted `u32` array?"
//! ([`lower_bound`](znippy_zoomies::stree32::STree32::lower_bound)). A viewport
//! query is turned into that shape the way every tile pyramid already does it:
//!
//! 1. Quantise each feature's centre into a `2^k × 2^k` lattice over the dataset
//!    bbox and interleave the bits — a **Morton key**.
//! 2. Sort features by Morton key once, at build time. Now spatial locality is
//!    array locality.
//! 3. A viewport covers `rows × cols` lattice cells. For **each row**, the cells
//!    `x0..=x1` lie between `morton(x0, y)` and `morton(x1, y)` — a contiguous
//!    key interval. Two `lower_bound`s per row give a contiguous slice.
//! 4. Those slices are a **superset** (Morton rows interleave), so each survivor
//!    gets one exact AABB test. `candidates` records exactly how many — the
//!    anti-lie field, the same discipline as
//!    [`NearestHit::candidates`](crate::legibility::NearestHit::candidates).
//!
//! The lattice level `k` is chosen **per query** so the viewport spans at most
//! [`MAX_ROWS`] rows. That is not a tuning constant dressed up as a bound: it is
//! the same "pick the pyramid level that matches the zoom" decision the
//! [`Hierarchy`](super::source::Hierarchy) makes, which is the point — *drill-down
//! and LOD are the same operation*.
//!
//! ## What this does NOT yet do
//! Features are indexed by their **centre**. A feature far wider than a lattice
//! cell (a motorway, a country boundary) can therefore be missed by a viewport
//! that clips it without containing its centre. [`CullIndex::build`] detects that
//! case and records it in [`CullIndex::oversize`]; the engine currently draws
//! oversize features unconditionally, which is correct but not yet cheap. The
//! fix is a per-level index (each feature entered at the level where it fits) —
//! the same `Bands` structure, applied to the culler.

use super::source::{WorldAabb, WorldPos};
use znippy_zoomies::stree32::STree32;

/// The most lattice rows one viewport query will decompose into. Above this the
/// lattice coarsens (fewer, larger cells) so the query count stays bounded and
/// the superset grows instead — the classic index trade, made explicit.
pub const MAX_ROWS: u32 = 64;

/// What a [`CullIndex::query`] found, and **how much work it did**.
#[derive(Clone, Debug, Default, PartialEq)]
pub struct CullReport {
    /// Features returned.
    pub visible: usize,
    /// Features that survived the Morton slice and got an exact AABB test. A
    /// linear scan reports `candidates == len`; an honest index reports far less.
    pub candidates: usize,
    /// Contiguous key intervals the viewport decomposed into (2 `lower_bound`s each).
    pub ranges: usize,
    /// Oversize features drawn without a spatial test (see the module note).
    pub oversize_forced: usize,
    /// The lattice level the query chose — the LOD/zoom band, made observable.
    pub level: u32,
}

impl CullReport {
    /// `state_json` fragment.
    #[must_use]
    pub fn state_json(&self) -> serde_json::Value {
        serde_json::json!({
            "visible": self.visible,
            "candidates": self.candidates,
            "ranges": self.ranges,
            "oversize_forced": self.oversize_forced,
            "level": self.level,
        })
    }
}

/// A static, Morton-ordered, S+tree-backed spatial index over world-space AABBs.
pub struct CullIndex {
    /// The dataset bbox the lattice is normalised against.
    bbox: WorldAabb,
    /// Morton keys, ascending. Index `i` here is `order[i]` in caller id space.
    keys: Vec<u32>,
    /// Caller feature ids, in Morton order.
    order: Vec<u32>,
    /// AABBs in Morton order, for the exact re-test.
    boxes: Vec<WorldAabb>,
    /// `None` for an empty dataset — `STree32::new` asserts on an empty slice,
    /// and an empty layer is a normal frame (a pane before its data arrives),
    /// not a crash.
    tree: Option<STree32>,
    /// Features wider than one lattice cell at the finest level — see module doc.
    oversize: Vec<u32>,
    /// Bits per axis in the stored keys (`LEVEL_BITS`).
    bits: u32,
}

// `STree32` is a packed SIMD array with no `Debug`; the index reports its shape
// rather than its bytes (which are never useful in a failure message anyway).
impl std::fmt::Debug for CullIndex {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.debug_struct("CullIndex")
            .field("len", &self.order.len())
            .field("oversize", &self.oversize.len())
            .field("bits", &self.bits)
            .finish()
    }
}

/// Bits per axis. `16 + 16 = 32` — exactly a `u32` Morton key, exactly the
/// S+tree's native width.
const LEVEL_BITS: u32 = 16;

#[inline]
fn part1by1(mut v: u32) -> u32 {
    v &= 0x0000_ffff;
    v = (v ^ (v << 8)) & 0x00ff_00ff;
    v = (v ^ (v << 4)) & 0x0f0f_0f0f;
    v = (v ^ (v << 2)) & 0x3333_3333;
    v = (v ^ (v << 1)) & 0x5555_5555;
    v
}

/// Interleave two 16-bit lattice coordinates into a 32-bit Morton key.
#[inline]
#[must_use]
pub fn morton(x: u32, y: u32) -> u32 {
    part1by1(x) | (part1by1(y) << 1)
}

impl CullIndex {
    /// Build over `(id, aabb)` pairs. `O(n log n)` once; queries are then
    /// `O(rows · log n + candidates)`.
    #[must_use]
    pub fn build(items: &[(u32, WorldAabb)]) -> Self {
        let mut bbox = WorldAabb::empty();
        for (_, b) in items {
            if b.min.x.is_finite() && b.min.y.is_finite() {
                bbox = bbox.union(*b);
            }
        }
        let n = 1u32 << LEVEL_BITS;
        let (sx, sy) = Self::scales(&bbox, n);

        let mut rows: Vec<(u32, u32, WorldAabb)> = Vec::with_capacity(items.len());
        let mut oversize = Vec::new();
        // A cell edge at the FINEST level, in world units.
        let cell_w = 1.0 / sx.max(f64::MIN_POSITIVE);
        let cell_h = 1.0 / sy.max(f64::MIN_POSITIVE);
        for &(id, b) in items {
            if !b.min.x.is_finite() || !b.min.y.is_finite() {
                continue;
            }
            let c = b.centre();
            let gx = (((c.x - bbox.min.x) * sx) as i64).clamp(0, (n - 1) as i64) as u32;
            let gy = (((c.y - bbox.min.y) * sy) as i64).clamp(0, (n - 1) as i64) as u32;
            rows.push((morton(gx, gy), id, b));
            if (b.max.x - b.min.x) > cell_w * (MAX_ROWS as f64) || (b.max.y - b.min.y) > cell_h * (MAX_ROWS as f64) {
                oversize.push(id);
            }
        }
        rows.sort_unstable_by_key(|r| r.0);

        let keys: Vec<u32> = rows.iter().map(|r| r.0).collect();
        let order: Vec<u32> = rows.iter().map(|r| r.1).collect();
        let boxes: Vec<WorldAabb> = rows.iter().map(|r| r.2).collect();
        let tree = (!keys.is_empty()).then(|| STree32::new(&keys));
        Self { bbox, keys, order, boxes, tree, oversize, bits: LEVEL_BITS }
    }

    fn scales(bbox: &WorldAabb, n: u32) -> (f64, f64) {
        let span_x = (bbox.max.x - bbox.min.x).max(f64::MIN_POSITIVE);
        let span_y = (bbox.max.y - bbox.min.y).max(f64::MIN_POSITIVE);
        (n as f64 / span_x, n as f64 / span_y)
    }

    #[must_use]
    pub fn len(&self) -> usize {
        self.order.len()
    }
    #[must_use]
    pub fn is_empty(&self) -> bool {
        self.order.is_empty()
    }
    /// Features too wide for the lattice — drawn unconditionally (see module doc).
    #[must_use]
    pub fn oversize(&self) -> &[u32] {
        &self.oversize
    }
    /// The dataset bbox the lattice normalises against.
    #[must_use]
    pub fn bbox(&self) -> WorldAabb {
        self.bbox
    }

    /// **Every feature whose AABB overlaps `view`**, ascending by Morton order,
    /// plus the work report. Exact: no false positives, no false negatives (for
    /// non-oversize features).
    pub fn query(&self, view: WorldAabb, out: &mut Vec<u32>) -> CullReport {
        out.clear();
        let mut rep = CullReport { level: self.bits, ..Default::default() };
        if self.order.is_empty() || view.is_empty() {
            return rep;
        }
        let full = 1u32 << self.bits;
        let (sx, sy) = Self::scales(&self.bbox, full);

        // Fine lattice cell range covered by the viewport.
        let clampi = |v: i64| v.clamp(0, (full - 1) as i64) as u32;
        let fx0 = clampi(((view.min.x - self.bbox.min.x) * sx).floor() as i64);
        let fx1 = clampi(((view.max.x - self.bbox.min.x) * sx).ceil() as i64);
        let fy0 = clampi(((view.min.y - self.bbox.min.y) * sy).floor() as i64);
        let fy1 = clampi(((view.max.y - self.bbox.min.y) * sy).ceil() as i64);

        // Coarsen until the row count fits the budget. Shifting BOTH axes by
        // `shift` keeps the Morton interleave valid: a coarse cell (gx>>s, gy>>s)
        // is exactly the contiguous key block `[m<<2s, ((m+1)<<2s)-1]`.
        let mut shift = 0u32;
        while ((fy1 >> shift) - (fy0 >> shift) + 1) > MAX_ROWS && shift < self.bits {
            shift += 1;
        }
        rep.level = self.bits - shift;

        let (cx0, cx1) = (fx0 >> shift, fx1 >> shift);
        let (cy0, cy1) = (fy0 >> shift, fy1 >> shift);
        let block = 2 * shift; // key bits consumed by one coarse cell

        let Some(tree) = self.tree.as_ref() else { return rep };

        // Per-row key intervals. Morton is NOT row-major, so consecutive rows'
        // intervals OVERLAP (row 0 of a 4-wide span is keys [0..5], row 1 is
        // [2..7]). Scanning them unmerged emits the same feature twice — a
        // duplicated draw, a duplicated pick, and a `visible` count that lies.
        // Merge first; the merge also shrinks the candidate set.
        let mut spans: Vec<(usize, usize)> = Vec::with_capacity((cy1 - cy0 + 1) as usize);
        for cy in cy0..=cy1 {
            let lo = (morton(cx0, cy) as u64) << block;
            let hi = (((morton(cx1, cy) as u64) + 1) << block) - 1;
            let a = tree.lower_bound(lo.min(u32::MAX as u64) as u32);
            // `hi + 1` as a u64 so the top-of-range case does not wrap.
            let hi_next = hi + 1;
            let b = if hi_next > u32::MAX as u64 {
                self.keys.len()
            } else {
                tree.lower_bound(hi_next as u32)
            };
            let b = b.min(self.keys.len());
            if a < b {
                spans.push((a, b));
            }
        }
        spans.sort_unstable();
        let mut merged: Vec<(usize, usize)> = Vec::with_capacity(spans.len());
        for (a, b) in spans {
            match merged.last_mut() {
                Some(last) if a <= last.1 => last.1 = last.1.max(b),
                _ => merged.push((a, b)),
            }
        }
        rep.ranges = merged.len();

        let mut candidates = 0usize;
        for (a, b) in merged {
            for i in a..b {
                candidates += 1;
                if self.boxes[i].overlaps_2d(&view) {
                    out.push(self.order[i]);
                }
            }
        }
        rep.candidates = candidates;
        rep.visible = out.len();
        rep
    }

    /// The reference answer — a full linear scan. Public so the S+tree path can
    /// be **red-proven** against it on random data rather than trusted.
    #[must_use]
    pub fn query_linear(&self, view: WorldAabb) -> Vec<u32> {
        let mut v: Vec<u32> =
            self.boxes.iter().enumerate().filter(|(_, b)| b.overlaps_2d(&view)).map(|(i, _)| self.order[i]).collect();
        v.sort_unstable();
        v
    }
}

/// The viewport, in world space, as an AABB — what the engine hands the culler.
#[must_use]
pub fn view_aabb(centre: WorldPos, half_w: f64, half_h: f64) -> WorldAabb {
    WorldAabb {
        min: WorldPos::flat(centre.x - half_w, centre.y - half_h),
        max: WorldPos::flat(centre.x + half_w, centre.y + half_h),
    }
}

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

    fn pt(x: f64, y: f64) -> WorldAabb {
        WorldAabb::point(WorldPos::flat(x, y))
    }

    /// Deterministic LCG — no dev-dependency for a random test.
    fn lcg(seed: &mut u64) -> f64 {
        *seed = seed.wrapping_mul(6_364_136_223_846_793_005).wrapping_add(1_442_695_040_888_963_407);
        ((*seed >> 11) as f64) / ((1u64 << 53) as f64)
    }

    #[test]
    fn stree_query_equals_linear_scan_on_random_points() {
        let mut seed = 0x5eed_1234_u64;
        let items: Vec<(u32, WorldAabb)> =
            (0..20_000u32).map(|i| (i, pt(lcg(&mut seed) * 1000.0, lcg(&mut seed) * 1000.0))).collect();
        let idx = CullIndex::build(&items);

        let mut out = Vec::new();
        let mut total_cand = 0usize;
        let mut total_vis = 0usize;
        for _ in 0..200 {
            let (x, y) = (lcg(&mut seed) * 1000.0, lcg(&mut seed) * 1000.0);
            let (w, h) = (lcg(&mut seed) * 60.0 + 1.0, lcg(&mut seed) * 60.0 + 1.0);
            let view = view_aabb(WorldPos::flat(x, y), w, h);
            let rep = idx.query(view, &mut out);
            let mut got = out.clone();
            got.sort_unstable();
            assert_eq!(got, idx.query_linear(view), "index disagreed with linear scan at {x},{y}");
            assert_eq!(rep.visible, got.len());
            total_cand += rep.candidates;
            total_vis += rep.visible;
        }
        // The anti-lie assertion: the index must actually be an index. A linear
        // scan would report 200 * 20_000 = 4_000_000 candidates.
        assert!(total_cand < 400_000, "index tested {total_cand} candidates — that is a linear scan in disguise");
        assert!(total_vis > 0, "the probe never hit anything — the test proves nothing");
    }

    #[test]
    fn empty_and_degenerate_inputs_do_not_panic() {
        let idx = CullIndex::build(&[]);
        let mut out = Vec::new();
        let rep = idx.query(view_aabb(WorldPos::flat(0.0, 0.0), 1.0, 1.0), &mut out);
        assert_eq!(rep.visible, 0);

        let same: Vec<(u32, WorldAabb)> = (0..100u32).map(|i| (i, pt(7.0, 7.0))).collect();
        let idx = CullIndex::build(&same);
        let rep = idx.query(view_aabb(WorldPos::flat(7.0, 7.0), 1.0, 1.0), &mut out);
        assert_eq!(rep.visible, 100, "every coincident point must survive a viewport containing them");
    }

    #[test]
    fn row_budget_bounds_the_range_count() {
        let mut seed = 99u64;
        let items: Vec<(u32, WorldAabb)> =
            (0..5_000u32).map(|i| (i, pt(lcg(&mut seed) * 1000.0, lcg(&mut seed) * 1000.0))).collect();
        let idx = CullIndex::build(&items);
        let mut out = Vec::new();
        // The whole world — the worst case for row decomposition.
        let rep = idx.query(view_aabb(WorldPos::flat(500.0, 500.0), 600.0, 600.0), &mut out);
        assert!(rep.ranges <= MAX_ROWS as usize + 1, "ranges={} exceeded the budget", rep.ranges);
        assert_eq!(rep.visible, 5_000, "a world-sized viewport must return everything");
    }

    #[test]
    fn oversize_features_are_flagged_not_silently_dropped() {
        let mut items: Vec<(u32, WorldAabb)> = (0..100u32).map(|i| (i, pt(i as f64, i as f64))).collect();
        items.push((
            999,
            WorldAabb { min: WorldPos::flat(0.0, 0.0), max: WorldPos::flat(100.0, 100.0) },
        ));
        let idx = CullIndex::build(&items);
        assert!(idx.oversize().contains(&999), "a dataset-spanning feature must be flagged oversize");
    }
}