facett-core 0.1.13

//! **Barnes–Hut N-body repulsion** — the shared, dimension-generic core that makes
//! every facett force-directed layout scale past ~30 000 nodes.
//!
//! The Fruchterman–Reingold repulsion term (every node pushes every other) is the
//! O(n²)-per-iteration curse behind the graph-viewer hang. This module replaces it
//! with the classic **Barnes–Hut** approximation: build a spatial tree (a **quadtree**
//! in 2-D, an **octree** in 3-D — the same code, `const D`), and when a whole cluster
//! of far-away nodes subtends a small angle it collapses to a single **center-of-mass**
//! under an opening angle `θ`. That drops the repulsion to **O(n log n)**; edge
//! attraction stays O(edges) in the caller.
//!
//! ## One routine, every copy (L5 reuse)
//! The exact same all-pairs FR loop was **twinned** across the workspace
//! (`facett-core::positions`, `facett-graph`, `facett-graphnav::layout`,
//! `facett-graph3d` `Layout3D::Force` + `callgraph_layouts::force2d`/`force3d`). This
//! is the ONE scalable repulsion kernel they all consume — [`repulsion_forces`] for a
//! full pass, [`relax_new`] for the incremental *absorb* path — so the fix lands once.
//! `D = 1` serves the 1-D depth relaxation (`force3d`), `D = 2` the plane, `D = 3` the
//! cube.
//!
//! ## Parallel to ~all cores (core-saturation LAW; rayon FORBIDDEN)
//! The tree is built serially (cheap, O(n log n)); the **per-node force query** is the
//! hot part and fans across every core through the znippy **gatling** fork-join engine
//! (`gatling_for_each`, one worker per core). Each node's net force is computed by a
//! single worker traversing the shared immutable tree in a **fixed** order, so the
//! result is **byte-deterministic regardless of the worker count** — the same property
//! the determinism goldens rely on.
//!
//! ## Additive / no-regression
//! [`BH_THRESHOLD`] is the node count at/above which a caller should switch to
//! Barnes–Hut; below it the callers keep their exact O(n²) loop, so every small-graph
//! determinism/golden test is untouched. With `θ = 0` Barnes–Hut degenerates to the
//! exact all-pairs sum (every node is its own leaf) — that is the parity the
//! correctness test pins.

use std::hash::{Hash, Hasher};

/// Node count at/above which a caller should switch from the exact O(n²) all-pairs
/// repulsion to the O(n log n) Barnes–Hut approximation. Below this the exact
/// algorithm runs (determinism goldens + small-graph parity unchanged).
pub const BH_THRESHOLD: usize = 2_000;

/// Barnes–Hut opening angle `θ`. A cell of width `w` seen from distance `d` is
/// approximated by its center-of-mass when `w/d < θ`. `0` ⇒ exact (every node a leaf);
/// `~0.5..1.0` is the usual speed/accuracy knee. `0.7` is a good default.
pub const BH_THETA: f32 = 0.7;

/// Max subdivision depth — a guard against unbounded recursion on (near-)coincident
/// points; beyond it a cell holds several bodies as a merged center-of-mass.
const MAX_DEPTH: usize = 48;

/// Distance floor (mirrors the exact FR loops' `.max(1e-3)`), so coincident/self
/// bodies never blow the `1/d²` term up.
const MIN_DIST: f32 = 1e-3;

/// Widest child fan-out we index (`2^3` for the 3-D octree); `1 << D` are live.
const MAX_CH: usize = 8;

/// One cell of the Barnes–Hut tree (a square in 2-D, a cube in 3-D).
#[derive(Clone)]
struct Cell<const D: usize> {
    /// Cell centre (geometric).
    center: [f32; D],
    /// Half the cell's side length.
    half: f32,
    /// Center of mass of the bodies under this cell.
    com: [f32; D],
    /// Body count under this cell (each body has unit mass).
    mass: f32,
    /// Body index when this is a leaf holding exactly one body; `-1` otherwise
    /// (empty, internal, or a merged multi-body leaf at the depth cap).
    body: i32,
    /// `1 << D` child cell indices into the arena (`-1` = absent). All `-1` ⇒ leaf.
    children: [i32; MAX_CH],
    /// Whether this cell is a leaf (no children).
    leaf: bool,
}

impl<const D: usize> Cell<D> {
    fn new(center: [f32; D], half: f32) -> Self {
        Self { center, half, com: [0.0; D], mass: 0.0, body: -1, children: [-1; MAX_CH], leaf: true }
    }
}

/// The child slot (0..`2^D`) for `pb` within a cell centred at `center`: bit `d` set
/// when `pb` is on the high side of the centre along dim `d`.
#[inline]
fn child_index<const D: usize>(center: [f32; D], pb: [f32; D]) -> usize {
    let mut idx = 0usize;
    for d in 0..D {
        if pb[d] >= center[d] {
            idx |= 1 << d;
        }
    }
    idx
}

/// Allocate child cell `idx` of the cell `(pc, ph)` into `cells`, returning its index
/// (half the side, centre offset per quadrant/octant).
fn arena_new_child<const D: usize>(cells: &mut Vec<Cell<D>>, pc: [f32; D], ph: f32, idx: usize) -> usize {
    let ch = ph * 0.5;
    let mut center = [0.0f32; D];
    for d in 0..D {
        let hi = (idx >> d) & 1 == 1;
        center[d] = pc[d] + if hi { ch } else { -ch };
    }
    cells.push(Cell::new(center, ch));
    cells.len() - 1
}

/// Insert body `b` (at `pb`, global index into `all_pos`) into the subtree rooted at
/// cell `ci` of the arena `cells`. The single, shared insertion routine both the serial
/// and the parallel builds consume (L5): identical center-of-mass folding, leaf
/// subdivision, and depth-cap merging — so every build path yields the same tree shape.
fn arena_insert<const D: usize>(
    cells: &mut Vec<Cell<D>>,
    all_pos: &[[f32; D]],
    ci: usize,
    b: usize,
    pb: [f32; D],
    depth: usize,
) {
    // A body passes through this cell: fold it into the running center of mass.
    let m = cells[ci].mass;
    for d in 0..D {
        cells[ci].com[d] = (cells[ci].com[d] * m + pb[d]) / (m + 1.0);
    }
    cells[ci].mass = m + 1.0;

    // First body to land here → an occupied leaf.
    if cells[ci].mass == 1.0 {
        cells[ci].body = b as i32;
        cells[ci].leaf = true;
        return;
    }
    // Depth cap / degenerate cell: keep as a merged multi-body leaf (approximated by
    // its center of mass). Never recurse forever on coincident points.
    if depth >= MAX_DEPTH || cells[ci].half <= MIN_DIST {
        cells[ci].body = -1;
        cells[ci].leaf = true;
        return;
    }
    // Was a single-body leaf: subdivide and push the resident body down first.
    if cells[ci].leaf && cells[ci].body >= 0 {
        let ob = cells[ci].body as usize;
        let obp = all_pos[ob];
        cells[ci].body = -1;
        cells[ci].leaf = false;
        arena_route(cells, all_pos, ci, ob, obp, depth);
    }
    // Route the new body into its quadrant/octant.
    arena_route(cells, all_pos, ci, b, pb, depth);
}

/// Send body `b` into the correct child of internal cell `ci`, creating it lazily.
fn arena_route<const D: usize>(
    cells: &mut Vec<Cell<D>>,
    all_pos: &[[f32; D]],
    ci: usize,
    b: usize,
    pb: [f32; D],
    depth: usize,
) {
    let idx = child_index(cells[ci].center, pb);
    let existing = cells[ci].children[idx];
    let child = if existing < 0 {
        let (pc, ph) = (cells[ci].center, cells[ci].half);
        let nc = arena_new_child(cells, pc, ph, idx);
        cells[ci].children[idx] = nc as i32;
        nc
    } else {
        existing as usize
    };
    arena_insert(cells, all_pos, child, b, pb, depth + 1);
}

/// A built Barnes–Hut tree over a fixed set of body positions. Immutable after build —
/// safe to share across gatling workers (`Sync`). Held as a **forest of arenas**: each
/// `arenas[b]` is one subtree (its cell `0` the root, child indices local to it). A
/// serial [`build`](NBodyTree::build) is a single-arena forest; the parallel
/// [`build_parallel`](NBodyTree::build_parallel) is one arena per spatial bucket, each
/// built on its own core and kept AS-IS — no serial splice/copy into a shared arena, so
/// the only serial residue of a build is the O(n) bounds + bucketing.
pub struct NBodyTree<const D: usize> {
    arenas: Vec<Vec<Cell<D>>>,
}

impl<const D: usize> NBodyTree<D> {
    /// Number of live children per internal cell (`4` for a quadtree, `8` for an octree).
    #[inline]
    fn nch() -> usize {
        1usize << D
    }

    /// The tight bounding cube over `pos` as `(center, half)` — the root cell every
    /// build starts from. `half` is the max half-span across dims (a hair padded so
    /// boundary points fall strictly inside).
    ///
    /// Kept **serial**: the per-dim min/max is memory-bound and only ~30 µs at 30k, so
    /// fanning it out (`gatling_reduce`) costs more in per-iteration thread-spawn than it
    /// saves at the ≤100k target — a `gatling_reduce` bounds fold was measured net-negative
    /// here and would only pay past millions of bodies. (The reduce primitive *is* the
    /// right tool for a big-accumulator scatter — see [`fr_step_scatter`].)
    fn root_bounds(pos: &[[f32; D]]) -> ([f32; D], f32) {
        let mut mn = [f32::MAX; D];
        let mut mx = [f32::MIN; D];
        for p in pos {
            for d in 0..D {
                mn[d] = mn[d].min(p[d]);
                mx[d] = mx[d].max(p[d]);
            }
        }
        let mut span = 0.0f32;
        let mut center = [0.0f32; D];
        for d in 0..D {
            center[d] = 0.5 * (mn[d] + mx[d]);
            span = span.max(mx[d] - mn[d]);
        }
        (center, (span * 0.5).max(MIN_DIST) * 1.0001)
    }

    /// Build the tree over `pos` (one body per position). O(n log n), **serial** — the
    /// reference build (small graphs + the parity golden). See
    /// [`build_parallel`](Self::build_parallel) for the core-saturating build at scale.
    #[must_use]
    pub fn build(pos: &[[f32; D]]) -> Self {
        if pos.is_empty() {
            return Self { arenas: Vec::new() };
        }
        let (center, half) = Self::root_bounds(pos);
        let mut cells = vec![Cell::new(center, half)];
        for (i, &p) in pos.iter().enumerate() {
            arena_insert(&mut cells, pos, 0, i, p, 0);
        }
        Self { arenas: vec![cells] }
    }

    /// **Parallel Barnes–Hut build** (core-saturating): partition space into a fixed
    /// grid of `(2^D)^pd` buckets, then build every bucket's subtree **on its own core**
    /// via gatling and splice the results into one arena as a **forest** of bucket roots.
    ///
    /// This removes the ~40 % serial `build` stall that otherwise caps core-utilisation:
    /// the only serial work left is the O(n) bounds + bucketing (cheap). Each bucket
    /// subtree is built serially by a single worker over that bucket's bodies **in
    /// ascending global-index order**, and buckets are spliced **in bucket-index order**,
    /// so the whole structure is **deterministic for any worker count**. Global body
    /// indices are preserved (self-exclusion in [`force`](Self::force) still works). The
    /// forest tests each bucket's far-field separately — a strictly finer (never coarser)
    /// opening than the single tree, so it stays within the same accuracy envelope and is
    /// EXACT at `theta = 0` (every body its own leaf), matching the serial build's field
    /// to float tolerance.
    #[must_use]
    pub fn build_parallel(pos: &[[f32; D]]) -> Self {
        let n = pos.len();
        // Small graphs (or a would-be trivial partition): the serial build is cheaper
        // than the fan-out, and keeps tiny inputs on the exact reference path.
        if n < 4_096 {
            return Self::build(pos);
        }
        let (root_center, root_half) = Self::root_bounds(pos);

        // Partition depth: pick the shallowest grid with ≥ 4× workers buckets (good load
        // balance) but never so deep the traversal stack (per-root) could overflow.
        let nch = 1usize << D;
        let want = znippy_zoomies::gatling_forkjoin::default_workers().saturating_mul(4).max(nch);
        let mut pd = 1usize;
        while nch.pow((pd + 1) as u32) <= want && pd < 4 {
            pd += 1;
        }
        let buckets = nch.pow(pd as u32);

        // Assign every body to its depth-`pd` bucket (the descent bit-path with level 0
        // as the most-significant group). This O(n·pd) pass is cheap (pd ≤ 4) — kept
        // serial: fanning such light per-item work out costs more in dispatch than it
        // saves at the 30k target (measured), though it stays a minor residue at 100k.
        let mut bucket_of = vec![0u32; n];
        for (i, &p) in pos.iter().enumerate() {
            let mut c = root_center;
            let mut h = root_half;
            let mut idx = 0usize;
            for _ in 0..pd {
                h *= 0.5;
                let ci = child_index(c, p);
                for d in 0..D {
                    c[d] += if (ci >> d) & 1 == 1 { h } else { -h };
                }
                idx = idx * nch + ci;
            }
            bucket_of[i] = idx as u32;
        }
        // Group body indices per bucket (stable ascending order → deterministic).
        let mut groups: Vec<Vec<u32>> = vec![Vec::new(); buckets];
        for (i, &b) in bucket_of.iter().enumerate() {
            groups[b as usize].push(i as u32);
        }

        // Reconstruct a bucket's cell centre/half from its index (level 0 = MSB group).
        let bucket_cell = |b: usize| -> ([f32; D], f32) {
            let mut c = root_center;
            let mut h = root_half;
            for l in 0..pd {
                h *= 0.5;
                let shift = (pd - 1 - l) * D;
                let ci = (b >> shift) & (nch - 1);
                for d in 0..D {
                    c[d] += if (ci >> d) & 1 == 1 { h } else { -h };
                }
            }
            (c, h)
        };

        // Build every bucket subtree in parallel: each returns its own arena whose cell
        // 0 is the bucket root (children indices LOCAL to that arena). The arenas are
        // kept as the forest AS-IS — no serial splice — so the build's only serial cost
        // is the O(n) bounds + bucketing above (the parallel query walks each arena).
        let arenas: Vec<Vec<Cell<D>>> = znippy_zoomies::gatling_forkjoin::gatling_for_each(buckets, 0, |b| {
            let bodies = &groups[b];
            if bodies.is_empty() {
                return Vec::new();
            }
            let (center, half) = bucket_cell(b);
            let mut cells = vec![Cell::new(center, half)];
            for &gi in bodies {
                let gi = gi as usize;
                arena_insert(&mut cells, pos, 0, gi, pos[gi], pd);
            }
            cells
        });
        Self { arenas }
    }

    /// The **net Barnes–Hut repulsion** on the body `self_idx` at position `p`, with
    /// FR ideal-length `k` and opening angle `theta`. Matches the exact all-pairs sum
    /// `Σ_j (p−p_j)·k²/|p−p_j|²` when `theta == 0` (every leaf a single body). Walks each
    /// forest root with a fixed per-root DFS, so the summation order — and thus the
    /// result — is identical for a serial (single-root) or parallel (multi-root) build.
    #[must_use]
    pub fn force(&self, p: [f32; D], self_idx: usize, k: f32, theta: f32) -> [f32; D] {
        let mut acc = [0.0f32; D];
        let k2 = k * k;
        let theta2 = theta * theta;
        // Explicit stack (no recursion in the hot per-node query), reset per arena. A DFS
        // frontier over one tree of depth ≤ MAX_DEPTH branching ≤ 2^D holds well under
        // this fixed capacity — so no per-query allocation, regardless of forest width.
        const CAP: usize = MAX_DEPTH * MAX_CH + 16;
        let mut stack: [usize; CAP] = [0; CAP];
        for cells in &self.arenas {
            if cells.is_empty() {
                continue;
            }
            stack[0] = 0; // each arena's root is its cell 0
            let mut top = 1usize;
            while top > 0 {
                top -= 1;
                let ci = stack[top];
                let c = &cells[ci];
                if c.mass == 0.0 {
                    continue;
                }
                // A single-body leaf that IS the query body contributes nothing.
                if c.leaf && c.body == self_idx as i32 {
                    continue;
                }
                let mut dv = [0.0f32; D];
                let mut d2 = 0.0f32;
                for d in 0..D {
                    let dd = p[d] - c.com[d];
                    dv[d] = dd;
                    d2 += dd * dd;
                }
                let width = c.half * 2.0;
                // Center-of-mass approximation for a leaf, or when the cell is far enough
                // (`width/dist < theta` ⇔ `width² < theta²·d²`). Otherwise descend.
                if c.leaf || width * width < theta2 * d2 {
                    let dist = d2.sqrt().max(MIN_DIST);
                    let f = k2 / (dist * dist) * c.mass; // (k²/d²)·mass; direction is dv
                    for d in 0..D {
                        acc[d] += dv[d] * f;
                    }
                } else {
                    for ch in 0..Self::nch() {
                        let cc = c.children[ch];
                        if cc >= 0 {
                            debug_assert!(top < CAP, "Barnes–Hut traversal stack overflow");
                            stack[top] = cc as usize;
                            top += 1;
                        }
                    }
                }
            }
        }
        acc
    }
}

/// Build the tree over `pos` and return the **Barnes–Hut repulsion force on every
/// node**, computed in parallel across all cores (gatling fork-join; one worker per
/// core). `k` is the Fruchterman–Reingold ideal edge length; `theta` the opening
/// angle ([`BH_THETA`]). Deterministic: identical output for any worker count.
///
/// This is the O(n log n) drop-in for the `for i { for j>i { … } }` repulsion double
/// loop — the caller keeps its own seed, cooling schedule, edge attraction, and
/// normalisation, and swaps only the repulsion term.
#[must_use]
pub fn repulsion_forces<const D: usize>(pos: &[[f32; D]], k: f32, theta: f32) -> Vec<[f32; D]> {
    let n = pos.len();
    if n == 0 {
        return Vec::new();
    }
    let tree = NBodyTree::<D>::build_parallel(pos);
    // `0` workers ⇒ one per core (core-saturation LAW). Each index is computed by a
    // single worker over the shared immutable tree, so the float sums are reproducible.
    znippy_zoomies::gatling_forkjoin::gatling_for_each(n, 0, |i| tree.force(pos[i], i, k, theta))
}

/// **Single-core** Barnes–Hut repulsion — serial build + serial per-node query, no
/// gatling anywhere. This is the honest one-core baseline the benchmark measures the
/// parallel [`repulsion_forces`] / [`fr_step`] speedup against; it is NOT on the shipped
/// hot path. Same field as the parallel version (within float tolerance).
#[must_use]
pub fn repulsion_forces_serial<const D: usize>(pos: &[[f32; D]], k: f32, theta: f32) -> Vec<[f32; D]> {
    let n = pos.len();
    if n == 0 {
        return Vec::new();
    }
    let tree = NBodyTree::<D>::build(pos);
    (0..n).map(|i| tree.force(pos[i], i, k, theta)).collect()
}

/// **One fully-parallel Fruchterman–Reingold step**, fanned across every core through
/// gatling — the drop-in that turns a single-core layout driver into a core-saturating
/// one. It fuses the three per-iteration phases, each over all cores:
///
/// 1. **Repulsion** — the Barnes–Hut O(n log n) field ([`repulsion_forces`], parallel).
/// 2. **Attraction + integrate** — a SINGLE `gatling_for_each` over the `n` nodes:
///    each worker owns a disjoint set of node indices and, for its node `i`, *gathers*
///    the spring pull from every incident edge (via the pre-built neighbour list
///    `adj[i]`) and writes the temperature-capped new position. Gathering per-node
///    (instead of scattering per-edge into `disp[a]` **and** `disp[b]`) removes the
///    write hazard that blocks a naive edge-loop parallelisation — no locks, no partial
///    buffers, disjoint output. Each node's neighbour sum runs in the FIXED `adj[i]`
///    order, so the result is deterministic for any worker count.
///
/// `adj[i]` must list node `i`'s neighbours with **multiplicity preserved** (each input
/// edge `(a,b)` contributes `b` to `adj[a]` and `a` to `adj[b]`), so the gathered
/// attraction reproduces the serial per-edge scatter term-for-term (only the summation
/// order differs — within float tolerance). `k` is the FR ideal length, `theta` the
/// opening angle, `temp` the current cooling cap. Returns the new positions.
#[must_use]
pub fn fr_step<const D: usize>(
    pos: &[[f32; D]],
    adj: &[Vec<u32>],
    k: f32,
    theta: f32,
    temp: f32,
) -> Vec<[f32; D]> {
    let n = pos.len();
    if n == 0 {
        return Vec::new();
    }
    // Build the Barnes–Hut tree ONCE per step, across all cores (the forest build). The
    // two remaining serial phases of the old driver (attraction scatter, integrate) are
    // fused into the single fan-out below, so a step is exactly TWO parallel passes:
    // (1) the parallel build, (2) this fused per-node kernel — no serial stall between.
    let tree = NBodyTree::<D>::build_parallel(pos);
    // Fused pass — repulsion query + attraction gather + temperature-capped integrate,
    // one worker per core. Node `i` is written by exactly one worker → disjoint,
    // hazard-free, deterministic (fixed neighbour order + fixed per-root DFS order).
    znippy_zoomies::gatling_forkjoin::gatling_for_each(n, 0, |i| {
        let pi = pos[i];
        let mut disp = tree.force(pi, i, k, theta);
        // Spring attraction, gathered from every incident edge (fixed neighbour order).
        for &jn in &adj[i] {
            let j = jn as usize;
            if j >= n {
                continue;
            }
            let pj = pos[j];
            let mut d2 = 0.0f32;
            let mut dv = [0.0f32; D];
            for d in 0..D {
                let dd = pi[d] - pj[d];
                dv[d] = dd;
                d2 += dd * dd;
            }
            let dist = d2.sqrt().max(MIN_DIST);
            // FR attraction magnitude dist²/k along −dv/dist ⇒ pull toward j: −dv·dist/k.
            let f = dist / k;
            for d in 0..D {
                disp[d] -= dv[d] * f;
            }
        }
        // Integrate, capped by the cooling temperature (matches the serial applier).
        let mut dl = 0.0f32;
        for d in 0..D {
            dl += disp[d] * disp[d];
        }
        let dl = dl.sqrt().max(MIN_DIST);
        let cap = dl.min(temp) / dl;
        let mut np = pi;
        for d in 0..D {
            np[d] += disp[d] * cap;
        }
        np
    })
}

/// **Scatter-reduce Fruchterman–Reingold step** — the same physics as [`fr_step`], but
/// the attraction phase is expressed as a genuine **parallel scatter-reduce** over the
/// EDGES via the znippy [`gatling_reduce`](znippy_zoomies::gatling_forkjoin::gatling_reduce)
/// primitive, instead of the per-node gather. Each worker scatters its edge-partition's
/// spring pulls into a **private** displacement buffer (so the `disp[a] += …; disp[b] += …`
/// write hazard is contained to one worker's buffer — no lock, no aliasing), and the
/// partials are summed element-wise in fixed partition order (deterministic). This does
/// the attraction distance work **once per edge** (vs twice — once per endpoint — in the
/// gather), at the cost of `workers` full-size partial buffers + a merge pass.
///
/// Provided as the scatter-reduce lane the write-hazard case calls for; [`fr_step`]'s
/// gather is the shipped default (measured faster: it fuses repulsion+attraction+integrate
/// into ONE pass with no partial buffers, whereas this is four passes). Both are
/// deterministic and match the serial reference within float tolerance.
#[must_use]
pub fn fr_step_scatter<const D: usize>(
    pos: &[[f32; D]],
    edges: &[(usize, usize)],
    k: f32,
    theta: f32,
    temp: f32,
) -> Vec<[f32; D]> {
    let n = pos.len();
    if n == 0 {
        return Vec::new();
    }
    // Pass 1 — parallel forest build + per-node repulsion query.
    let rep = repulsion_forces::<D>(pos, k, theta);
    // Pass 2 — attraction as a scatter-reduce over edges: each worker folds its edge
    // range into a private `disp` buffer; the buffers are merged in fixed order.
    let attr = znippy_zoomies::gatling_forkjoin::gatling_reduce(
        edges.len(),
        0,
        || vec![[0.0f32; D]; n],
        |buf, e| {
            let (a, b) = edges[e];
            if a >= n || b >= n || a == b {
                return;
            }
            let mut d2 = 0.0f32;
            let mut dv = [0.0f32; D];
            for d in 0..D {
                let dd = pos[a][d] - pos[b][d];
                dv[d] = dd;
                d2 += dd * dd;
            }
            let dist = d2.sqrt().max(MIN_DIST);
            let f = dist / k; // FR attraction dist²/k along the unit edge; pull both ends in
            for d in 0..D {
                let push = dv[d] * f;
                buf[a][d] -= push;
                buf[b][d] += push;
            }
        },
        |acc, part| {
            for i in 0..n {
                for d in 0..D {
                    acc[i][d] += part[i][d];
                }
            }
        },
    );
    // Pass 3 — per-node temperature-capped integrate (disjoint output).
    znippy_zoomies::gatling_forkjoin::gatling_for_each(n, 0, |i| {
        let mut disp = rep[i];
        for d in 0..D {
            disp[d] += attr[i][d];
        }
        let mut dl = 0.0f32;
        for d in 0..D {
            dl += disp[d] * disp[d];
        }
        let dl = dl.sqrt().max(MIN_DIST);
        let cap = dl.min(temp) / dl;
        let mut np = pos[i];
        for d in 0..D {
            np[d] += disp[d] * cap;
        }
        np
    })
}

/// Build undirected neighbour lists (**multiplicity preserved**) from a directed edge
/// list over `n` nodes — the `adj` input [`fr_step`] gathers attraction from. Each edge
/// `(a,b)` with `a ≠ b`, both in range, pushes `b` onto `adj[a]` and `a` onto `adj[b]`.
/// Self-loops and out-of-range endpoints are skipped (mirroring the serial edge loop).
#[must_use]
pub fn adjacency(n: usize, edges: &[(usize, usize)]) -> Vec<Vec<u32>> {
    let mut adj = vec![Vec::new(); n];
    for &(a, b) in edges {
        if a >= n || b >= n || a == b {
            continue;
        }
        adj[a].push(b as u32);
        adj[b].push(a as u32);
    }
    adj
}

/// **Incremental absorb relaxation.** Grow a laid-out graph by relaxing ONLY the new
/// nodes (indices `new_lo..pos.len()`) while every existing node stays frozen in place:
/// a few local Fruchterman–Reingold iterations (Barnes–Hut repulsion against the whole
/// cloud + spring attraction along incident edges), each new node's step capped by
/// `step`. Keeps the existing layout stable so a consumer can `absorb` a batch of
/// nodes without a full O(n log n) rebuild. `k` is the FR ideal length, `theta` the
/// opening angle, `iters` the local-relaxation budget (small, e.g. 8..20).
pub fn relax_new<const D: usize>(
    pos: &mut [[f32; D]],
    edges: &[(usize, usize)],
    new_lo: usize,
    k: f32,
    theta: f32,
    iters: usize,
    step: f32,
) {
    let n = pos.len();
    if n == 0 || new_lo >= n {
        return;
    }
    for _ in 0..iters {
        // Repulsion field over the whole cloud (parallel), read only for new nodes.
        let rep = repulsion_forces::<D>(pos, k, theta);
        let mut disp = vec![[0.0f32; D]; n];
        for i in new_lo..n {
            disp[i] = rep[i];
        }
        // Spring attraction along edges touching a new node (the moving endpoint pulls).
        for &(a, b) in edges {
            if a >= n || b >= n || a == b {
                continue;
            }
            let a_new = a >= new_lo;
            let b_new = b >= new_lo;
            if !a_new && !b_new {
                continue;
            }
            let mut dv = [0.0f32; D];
            let mut d2 = 0.0f32;
            for d in 0..D {
                let dd = pos[a][d] - pos[b][d];
                dv[d] = dd;
                d2 += dd * dd;
            }
            let dist = d2.sqrt().max(MIN_DIST);
            let f = dist / k; // FR attraction magnitude dist²/k, per unit direction dv/dist
            for d in 0..D {
                let push = dv[d] * f;
                if a_new {
                    disp[a][d] -= push;
                }
                if b_new {
                    disp[b][d] += push;
                }
            }
        }
        // Apply, moving only the new nodes, capped.
        for i in new_lo..n {
            let mut dl = 0.0f32;
            for d in 0..D {
                dl += disp[i][d] * disp[i][d];
            }
            let dl = dl.sqrt().max(MIN_DIST);
            let cap = dl.min(step) / dl;
            for d in 0..D {
                pos[i][d] += disp[i][d] * cap;
            }
        }
    }
}

/// A cheap **structural signature** of a graph layout request — `(layout key, node
/// count, edge list)` folded into a `u64`. A force layout is a pure function of this,
/// so a stateful consumer freezes its converged positions and only re-runs the
/// expensive layout when the signature changes (the *converge-once-then-freeze*
/// contract). Computing it is O(edges) — negligible next to an O(n log n) relayout.
#[must_use]
pub fn structure_sig(layout_key: &str, n: usize, edges: &[(usize, usize)]) -> u64 {
    let mut h = std::collections::hash_map::DefaultHasher::new();
    layout_key.hash(&mut h);
    n.hash(&mut h);
    edges.len().hash(&mut h);
    for e in edges {
        e.hash(&mut h);
    }
    h.finish()
}

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

    /// A deterministic, well-separated 2-D point cloud (no coincident points, so every
    /// body is its own leaf under `theta = 0`).
    fn cloud2(n: usize) -> Vec<[f32; 2]> {
        (0..n)
            .map(|i| {
                let a = std::f32::consts::TAU * i as f32 / n as f32;
                let r = 1.0 + (i as f32 * 0.37).sin() * 0.5;
                [r * a.cos(), r * a.sin()]
            })
            .collect()
    }

    /// Exact all-pairs repulsion (the reference the tree must match).
    fn naive2(pos: &[[f32; 2]], k: f32) -> Vec<[f32; 2]> {
        let n = pos.len();
        let k2 = k * k;
        let mut out = vec![[0.0f32; 2]; n];
        for i in 0..n {
            for j in 0..n {
                if i == j {
                    continue;
                }
                let dx = pos[i][0] - pos[j][0];
                let dy = pos[i][1] - pos[j][1];
                let dist = (dx * dx + dy * dy).sqrt().max(MIN_DIST);
                let f = k2 / (dist * dist);
                out[i][0] += dx * f;
                out[i][1] += dy * f;
            }
        }
        out
    }

    /// **RED-when-broken correctness:** with `theta = 0` Barnes–Hut is the EXACT
    /// all-pairs sum (every leaf a single body), so it must match the naive reference
    /// to float tolerance. A broken tree (mis-routed bodies, wrong com/mass) diverges
    /// wildly and trips this.
    #[test]
    fn bh_theta_zero_matches_naive_2d() {
        let pos = cloud2(400);
        let k = (4.0 / pos.len() as f32).sqrt();
        let bh = repulsion_forces::<2>(&pos, k, 0.0);
        let naive = naive2(&pos, k);
        let mut max_err = 0.0f32;
        for i in 0..pos.len() {
            for d in 0..2 {
                max_err = max_err.max((bh[i][d] - naive[i][d]).abs());
            }
        }
        assert!(max_err < 1e-3, "theta=0 Barnes–Hut must equal naive all-pairs; max_err={max_err}");
    }

    /// With a real opening angle the approximation stays CLOSE to the exact field
    /// (same direction, bounded relative error) — the accuracy the layout relies on.
    #[test]
    fn bh_theta_moderate_close_to_naive_2d() {
        let pos = cloud2(600);
        let k = (4.0 / pos.len() as f32).sqrt();
        let bh = repulsion_forces::<2>(&pos, k, 0.7);
        let naive = naive2(&pos, k);
        // Compare force magnitudes: the mean relative error must be small.
        let mut sum_rel = 0.0f32;
        let mut cnt = 0.0f32;
        for i in 0..pos.len() {
            let bm = (bh[i][0] * bh[i][0] + bh[i][1] * bh[i][1]).sqrt();
            let nm = (naive[i][0] * naive[i][0] + naive[i][1] * naive[i][1]).sqrt();
            if nm > 1e-4 {
                sum_rel += (bm - nm).abs() / nm;
                cnt += 1.0;
            }
        }
        let mean_rel = sum_rel / cnt.max(1.0);
        assert!(mean_rel < 0.15, "theta=0.7 stays within ~15% of exact; mean_rel={mean_rel}");
    }

    /// Determinism: the parallel per-node force is byte-identical run to run (the
    /// property the layout goldens depend on).
    #[test]
    fn repulsion_is_deterministic() {
        let pos = cloud2(1500);
        let k = (4.0 / pos.len() as f32).sqrt();
        let a = repulsion_forces::<2>(&pos, k, 0.7);
        let b = repulsion_forces::<2>(&pos, k, 0.7);
        assert_eq!(a, b, "parallel Barnes–Hut repulsion is deterministic");
    }

    /// **RED-when-broken (parallel build == serial build):** the core-saturating forest
    /// `build_parallel` must produce the SAME repulsion field as the reference serial
    /// `build`. At `theta = 0` both are the exact all-pairs sum, so they agree to tight
    /// tolerance; a broken parallel build (mis-bucketed body, wrong centre/half
    /// reconstruction, bad child-index remap on splice, dropped bucket) diverges and
    /// trips this. Run well above the 4 096 fan-out threshold so the parallel path is
    /// actually exercised.
    #[test]
    fn parallel_build_matches_serial_build() {
        let pos = cloud2(8_000);
        let k = (4.0 / pos.len() as f32).sqrt();
        let serial = NBodyTree::<2>::build(&pos);
        let parallel = NBodyTree::<2>::build_parallel(&pos);
        let mut max_err = 0.0f32;
        for i in 0..pos.len() {
            let fs = serial.force(pos[i], i, k, 0.0);
            let fp = parallel.force(pos[i], i, k, 0.0);
            for d in 0..2 {
                max_err = max_err.max((fs[d] - fp[d]).abs());
            }
        }
        assert!(max_err < 1e-3, "forest build matches serial build at theta=0; max_err={max_err}");
        // Total mass is conserved across the forest (every body bucketed exactly once):
        // each arena's root (cell 0) carries its bucket's body count.
        let mass: f32 = parallel.arenas.iter().filter(|a| !a.is_empty()).map(|a| a[0].mass).sum();
        assert!((mass - pos.len() as f32).abs() < 0.5, "forest holds every body; mass={mass}");
    }

    /// **RED-when-broken (scatter-reduce == gather):** one `fr_step_scatter` (attraction
    /// as a `gatling_reduce` scatter-reduce over edges) equals one `fr_step` (per-node
    /// gather) to float tolerance — the two parallel attraction schemes compute the same
    /// field. A broken scatter (lost partial, wrong endpoint sign) diverges and trips this.
    #[test]
    fn scatter_reduce_step_matches_gather_step() {
        let n = 6_000usize;
        let mut edges = Vec::new();
        for t in 0..n / 3 {
            let (a, b, c) = (t * 3, t * 3 + 1, t * 3 + 2);
            edges.push((a, b));
            edges.push((b, c));
            edges.push((c, a));
            edges.push((a, (a + 3) % n));
        }
        let pos: Vec<[f32; 2]> = (0..n)
            .map(|i| {
                let ang = std::f32::consts::TAU * i as f32 / n as f32;
                [ang.cos(), ang.sin()]
            })
            .collect();
        let k = (4.0 / n as f32).sqrt();
        let adj = adjacency(n, &edges);
        let gather = fr_step::<2>(&pos, &adj, k, BH_THETA, 0.05);
        let scatter = fr_step_scatter::<2>(&pos, &edges, k, BH_THETA, 0.05);
        let mut max_err = 0.0f32;
        for (g, s) in gather.iter().zip(&scatter) {
            max_err = max_err.max((g[0] - s[0]).abs()).max((g[1] - s[1]).abs());
        }
        assert!(max_err < 1e-4, "scatter-reduce == gather attraction; max_err={max_err}");
    }

    /// The forest build is byte-deterministic run to run (any worker count) — a
    /// per-node force sum is identical across builds.
    #[test]
    fn parallel_build_is_deterministic() {
        let pos = cloud2(6_000);
        let k = (4.0 / pos.len() as f32).sqrt();
        let a = repulsion_forces::<2>(&pos, k, 0.7);
        let t = NBodyTree::<2>::build_parallel(&pos);
        let b: Vec<[f32; 2]> = (0..pos.len()).map(|i| t.force(pos[i], i, k, 0.7)).collect();
        // `repulsion_forces` (which uses build_parallel) equals a fresh forest query.
        assert_eq!(a, b, "forest build + query is deterministic");
    }

    /// The 3-D octree also matches the exact all-pairs sum at `theta = 0`.
    #[test]
    fn bh_theta_zero_matches_naive_3d() {
        let n = 300;
        let golden = std::f32::consts::PI * (3.0 - 5.0_f32.sqrt());
        let pos: Vec<[f32; 3]> = (0..n)
            .map(|i| {
                let y = 1.0 - (i as f32 / (n - 1) as f32) * 2.0;
                let r = (1.0 - y * y).max(0.0).sqrt();
                let th = golden * i as f32;
                [th.cos() * r, y, th.sin() * r]
            })
            .collect();
        let k = 0.3;
        let k2 = k * k;
        let bh = repulsion_forces::<3>(&pos, k, 0.0);
        let mut max_err = 0.0f32;
        for i in 0..n {
            let mut naive = [0.0f32; 3];
            for j in 0..n {
                if i == j {
                    continue;
                }
                let mut dv = [0.0f32; 3];
                let mut d2 = 0.0f32;
                for d in 0..3 {
                    dv[d] = pos[i][d] - pos[j][d];
                    d2 += dv[d] * dv[d];
                }
                let dist = d2.sqrt().max(MIN_DIST);
                let f = k2 / (dist * dist);
                for d in 0..3 {
                    naive[d] += dv[d] * f;
                }
            }
            for d in 0..3 {
                max_err = max_err.max((bh[i][d] - naive[d]).abs());
            }
        }
        assert!(max_err < 1e-3, "3-D theta=0 octree must equal naive; max_err={max_err}");
    }

    /// `relax_new` moves only the appended nodes; every existing position is frozen.
    #[test]
    fn relax_new_freezes_existing_nodes() {
        let mut pos = cloud2(50);
        let frozen = pos.clone();
        let old_n = pos.len();
        // Append 5 new nodes near the origin, wired to a couple of existing ones.
        for i in 0..5 {
            pos.push([0.01 * i as f32, 0.02 * i as f32]);
        }
        let edges: Vec<(usize, usize)> = (0..5).map(|i| (old_n + i, i)).collect();
        let k = (4.0 / pos.len() as f32).sqrt();
        relax_new::<2>(&mut pos, &edges, old_n, k, 0.7, 10, 0.05);
        for i in 0..old_n {
            assert_eq!(pos[i], frozen[i], "existing node {i} stays frozen during absorb");
        }
        for i in old_n..pos.len() {
            assert!(pos[i][0].is_finite() && pos[i][1].is_finite(), "new node {i} is finite");
        }
    }

    /// The signature is stable for the same graph and changes when structure changes.
    #[test]
    fn structure_sig_tracks_changes() {
        let e1 = vec![(0usize, 1usize), (1, 2)];
        let e2 = vec![(0usize, 1usize), (1, 2), (2, 3)];
        assert_eq!(structure_sig("force", 3, &e1), structure_sig("force", 3, &e1));
        assert_ne!(structure_sig("force", 3, &e1), structure_sig("force", 4, &e2));
        assert_ne!(structure_sig("force", 3, &e1), structure_sig("circular", 3, &e1));
    }
}