facett-core/src/barnes_hut.rs
//! **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 (com, p) in cells[ci].com.iter_mut().zip(pb) {
*com = (*com * m + p) / (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, cd) in c.iter_mut().enumerate() {
*cd += 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, cd) in c.iter_mut().enumerate() {
*cd += 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 v in &disp {
dl += v * v;
}
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 v in &disp {
dl += v * v;
}
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
}
/// **One serial Fruchterman–Reingold step** — the single-core mirror of [`fr_step`]:
/// identical physics (repulsion + per-node attraction *gather* + temperature-capped
/// integrate), no gatling anywhere. Below [`BH_THRESHOLD`] the repulsion is the exact
/// O(n²) all-pairs loop (so small-graph layouts stay on the exact reference); at/above it
/// the serial Barnes–Hut tree ([`repulsion_forces_serial`]) takes over. Attraction is
/// gathered per node from `adj` (multiplicity preserved, fixed order — build it with
/// [`adjacency`]), so the result matches [`fr_step`] within float tolerance. `k` is the
/// FR ideal length, `theta` the opening angle (used only on the BH branch), `temp` the
/// current cooling cap.
///
/// This is the ONE serial FR step the small-graph / single-core layout drivers share
/// (`facett-graphnav`'s serial lane and `facett-graph3d`'s `force2d` both call it, each
/// supplying its own `k`/`temp`), so the loop body lives here once instead of being
/// re-inlined per leaf. Coincident nodes carry a zero direction vector, so — as in every
/// FR — they are separated by the seed, not this step (both callers seed on a ring).
#[must_use]
pub fn fr_step_serial<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();
}
// Repulsion — exact all-pairs below the BH threshold (keeps the small-graph exact
// reference byte-for-byte), serial Barnes–Hut above it (O(n log n)).
let mut disp = if n >= BH_THRESHOLD {
repulsion_forces_serial::<D>(pos, k, theta)
} else {
let mut d = vec![[0.0f32; D]; n];
for i in 0..n {
for j in (i + 1)..n {
let mut dv = [0.0f32; D];
let mut d2 = 0.0f32;
for c in 0..D {
let dd = pos[i][c] - pos[j][c];
dv[c] = dd;
d2 += dd * dd;
}
let dist = d2.sqrt().max(MIN_DIST);
// FR repulsion k²/dist along the unit vector dv/dist ⇒ dv·k²/dist².
let f = k * k / (dist * dist);
for c in 0..D {
let push = dv[c] * f;
d[i][c] += push;
d[j][c] -= push;
}
}
}
d
};
// Attraction — per-node gather over the neighbour list (fixed order ⇒ deterministic,
// matches [`fr_step`] term-for-term, only the summation order differs).
for i in 0..n {
let pi = pos[i];
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 c in 0..D {
let dd = pi[c] - pj[c];
dv[c] = dd;
d2 += dd * dd;
}
let dist = d2.sqrt().max(MIN_DIST);
// FR attraction dist²/k toward j ⇒ −dv·dist/k.
let f = dist / k;
for c in 0..D {
disp[i][c] -= dv[c] * f;
}
}
}
// Integrate — capped by the cooling temperature (identical applier to [`fr_step`]).
let mut out = pos.to_vec();
for i in 0..n {
let mut dl = 0.0f32;
for v in &disp[i] {
dl += v * v;
}
let dl = dl.sqrt().max(MIN_DIST);
let cap = dl.min(temp) / dl;
for c in 0..D {
out[i][c] += disp[i][c] * cap;
}
}
out
}
/// Rescale a 2-D position set to fit the `[-0.9, 0.9]²` box (aspect preserved): centre on
/// the bounding-box midpoint, divide by the larger span, scale by `1.8`. The single home
/// for the force-layout box-fit that `facett-graphnav` and `facett-graphview` share, so
/// their world lanes occupy the same world space.
pub fn normalize_box(pos: &mut [[f32; 2]]) {
let (mut x0, mut y0, mut x1, mut y1) = (f32::MAX, f32::MAX, f32::MIN, f32::MIN);
for p in pos.iter() {
x0 = x0.min(p[0]);
y0 = y0.min(p[1]);
x1 = x1.max(p[0]);
y1 = y1.max(p[1]);
}
let span = (x1 - x0).max(y1 - y0).max(1e-3);
let (cx, cy) = ((x0 + x1) * 0.5, (y0 + y1) * 0.5);
for p in pos.iter_mut() {
p[0] = (p[0] - cx) / span * 1.8;
p[1] = (p[1] - cy) / span * 1.8;
}
}
/// **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];
disp[new_lo..n].copy_from_slice(&rep[new_lo..n]);
// 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 v in &disp[i] {
dl += v * v;
}
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, &p) in pos.iter().enumerate() {
let fs = serial.force(p, i, k, 0.0);
let fp = parallel.force(p, 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, p) in pos.iter().enumerate().skip(old_n) {
assert!(p[0].is_finite() && p[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));
}
}