graph_core/

force.rs

//! Barnes-Hut 2D force simulation.
//!
//! Port of the d3-force semantics onto a quadtree-accelerated
//! n-body repulsion kernel, plus link springs and a center-of-mass
//! gravity. The result is deterministic given a seed so the WASM and
//! native call sites (Wave 2+) produce identical layouts for the
//! same graph.
//!
//! ## Single-threaded for Wave 1
//!
//! The naive per-tick traversal on 50k nodes is well under 5 ms on
//! a modern laptop core so we leave it single-threaded here. A
//! rayon-parallel variant would live behind `#[cfg(feature = "native")]`
//! once we have profile data to justify it. Wave 1's priority is
//! correctness + WASM-compatibility of the algorithm, not raw peak
//! throughput on a dense 100k-node graph.
//!
//! ## SIMD choice
//!
//! Positions are 2-float vectors, so per-element SIMD (`f32x2`)
//! buys nothing — the dot product over two floats is one FMA already.
//! We use `glam::Vec2` scalar math everywhere in this file and keep
//! `graph_core::util::dot_simd` for the wide-float kernels (cosine,
//! BM25, bloom) that need it. A future revision that batches the
//! center-gravity force over all positions at once and runs it as a
//! fused reduction could reach for `dot_simd`; the per-node force
//! loop wouldn't benefit.

use glam::{Vec2, Vec3};
use rand::{Rng, SeedableRng};
use rand_chacha::ChaCha8Rng;

use crate::force_3d::{OctNode, accumulate_repulsion_3d, build_octree};

/// Minimum quadtree cell size, in world units. Below this the tree
/// stops splitting and treats coincident bodies as a single
/// accumulated point mass. 1.0 is plenty given initial positions
/// are scattered on a radius-10+ disc.
const MIN_CELL_SIZE: f32 = 1.0;

/// Force simulation configuration. The defaults match d3-force's
/// defaults so graphs that used to render there look the same here.
#[derive(Debug, Clone, Copy)]
pub struct SimulationConfig {
    /// Negative values repel, positive values attract. d3 default -30.
    pub repulsion_strength: f32,
    /// Spring stiffness scaling on the link force (0..1 range is typical).
    pub link_strength: f32,
    /// Natural length of each link edge in world units.
    pub link_distance: f32,
    /// Center-of-mass gravity pulling nodes toward the origin.
    pub center_strength: f32,
    /// Per-tick velocity damping. 0.4 matches d3's default velocity decay.
    pub velocity_decay: f32,
    /// Barnes-Hut approximation threshold. Larger = faster but less
    /// accurate. d3 uses 0.9 under the name `theta`.
    pub theta: f32,
    /// RNG seed used to scatter initial positions deterministically.
    pub seed: u64,
    /// Dimensionality of the simulation: 2 (default, quadtree) or 3
    /// (octree). Phase 224 Wave 2 toggle — anything else is clamped
    /// to 2. Positions keep the 2D layout when `dimensions == 2` and
    /// switch to 3-vector positions (accessible via `positions_3d()`)
    /// when `dimensions == 3`.
    pub dimensions: u8,
}

impl Default for SimulationConfig {
    fn default() -> Self {
        Self {
            repulsion_strength: -30.0,
            link_strength: 0.1,
            link_distance: 30.0,
            center_strength: 0.05,
            velocity_decay: 0.4,
            theta: 0.9,
            seed: 0,
            dimensions: 2,
        }
    }
}

/// Single edge in the force graph. `weight` scales the link spring
/// force so heavier edges pull endpoints together faster.
#[derive(Debug, Clone, Copy)]
struct Edge {
    src: u32,
    dst: u32,
    weight: f32,
}

/// Public force simulation handle. Owns positions, velocities, edges,
/// and a scratch quadtree that gets rebuilt every tick.
///
/// In 2D mode (`config.dimensions == 2`, the default) only the
/// `positions`/`velocities`/`forces` vectors and the `tree` quadtree
/// are used. In 3D mode the parallel `positions_3d`/`velocities_3d`/
/// `forces_3d` vectors are populated alongside an `oct_tree` octree,
/// and `positions` mirrors the x/y pair of each 3D position so the
/// `positions()` accessor stays zero-copy for the common 2D readers.
pub struct Simulation {
    config: SimulationConfig,
    positions: Vec<Vec2>,
    velocities: Vec<Vec2>,
    // Per-tick scratch force accumulator, kept to avoid reallocating.
    forces: Vec<Vec2>,
    edges: Vec<Edge>,
    // Scratch arena for the Barnes-Hut quadtree. Cleared and refilled
    // every tick; retained here to amortise allocation.
    tree: Vec<QuadNode>,
    // 3D-only state. Empty vectors when `config.dimensions == 2`.
    positions_3d: Vec<Vec3>,
    velocities_3d: Vec<Vec3>,
    forces_3d: Vec<Vec3>,
    oct_tree: Vec<OctNode>,
}

impl Simulation {
    /// Build a simulation with `n_nodes` nodes. Initial positions are
    /// scattered uniformly inside a circle of radius `sqrt(n) * 10`
    /// using a deterministic ChaCha8 RNG seeded from `config.seed`.
    ///
    /// When `config.dimensions == 3`, the 3D `positions_3d` field is
    /// also populated with points scattered inside a ball of the same
    /// radius; the 2D `positions` field mirrors the x/y slice of each
    /// 3D point so consumers of the 2D accessor continue to see a
    /// sensible layout.
    pub fn new(n_nodes: u32, config: SimulationConfig) -> Self {
        let n = n_nodes as usize;
        let mut rng = ChaCha8Rng::seed_from_u64(config.seed);
        let radius = (n_nodes as f32).max(1.0).sqrt() * 10.0;

        // Clamp dimensions; anything outside {2, 3} degrades gracefully to 2D.
        let dims = if config.dimensions == 3 { 3u8 } else { 2u8 };
        let mut cfg = config;
        cfg.dimensions = dims;

        let mut positions = Vec::with_capacity(n);
        let mut positions_3d: Vec<Vec3> = if dims == 3 {
            Vec::with_capacity(n)
        } else {
            Vec::new()
        };

        for _ in 0..n {
            if dims == 3 {
                // Rejection-sample a point in the unit ball, then scale.
                // Same pattern as the 2D case; keeps the disc/ball
                // analogue visually consistent.
                loop {
                    let x: f32 = rng.gen_range(-1.0..=1.0);
                    let y: f32 = rng.gen_range(-1.0..=1.0);
                    let z: f32 = rng.gen_range(-1.0..=1.0);
                    if x * x + y * y + z * z <= 1.0 {
                        let p = Vec3::new(x * radius, y * radius, z * radius);
                        positions_3d.push(p);
                        positions.push(Vec2::new(p.x, p.y));
                        break;
                    }
                }
            } else {
                // Rejection-sample a point in the unit disc, then scale.
                // Gives a visibly uniform cloud rather than the polar
                // artefacts you get from sampling (r, θ) independently.
                loop {
                    let x: f32 = rng.gen_range(-1.0..=1.0);
                    let y: f32 = rng.gen_range(-1.0..=1.0);
                    if x * x + y * y <= 1.0 {
                        positions.push(Vec2::new(x * radius, y * radius));
                        break;
                    }
                }
            }
        }

        let (velocities_3d, forces_3d) = if dims == 3 {
            (vec![Vec3::ZERO; n], vec![Vec3::ZERO; n])
        } else {
            (Vec::new(), Vec::new())
        };

        Self {
            config: cfg,
            positions,
            velocities: vec![Vec2::ZERO; n],
            forces: vec![Vec2::ZERO; n],
            edges: Vec::new(),
            tree: Vec::new(),
            positions_3d,
            velocities_3d,
            forces_3d,
            oct_tree: Vec::new(),
        }
    }

    /// Replace the edge set. Each edge is `(src, dst, weight)`.
    /// Out-of-range or self-loop indices are silently dropped — callers
    /// who care should validate upstream.
    pub fn set_edges(&mut self, edges: &[(u32, u32, f32)]) {
        let n = self.positions.len() as u32;
        self.edges.clear();
        self.edges.reserve(edges.len());
        for &(src, dst, weight) in edges {
            if src < n && dst < n && src != dst {
                self.edges.push(Edge { src, dst, weight });
            }
        }
    }

    /// Number of nodes in the simulation.
    #[inline]
    pub fn node_count(&self) -> usize {
        self.positions.len()
    }

    /// Borrow positions as a flat `&[Vec2]` for zero-copy consumption
    /// by downstream bindings.
    #[inline]
    pub fn positions(&self) -> &[Vec2] {
        &self.positions
    }

    /// Borrow 3D positions as a flat `&[Vec3]`. Returns an empty slice
    /// when the simulation was constructed in 2D mode. Phase 224 Wave 2.
    #[inline]
    pub fn positions_3d(&self) -> &[Vec3] {
        &self.positions_3d
    }

    /// Advance the simulation by one tick of size `dt`. Routes to the
    /// 3D octree pass when `config.dimensions == 3`, otherwise runs
    /// the 2D quadtree path.
    pub fn tick(&mut self, dt: f32) {
        if self.config.dimensions == 3 {
            self.tick_3d(dt);
            return;
        }
        self.tick_2d(dt);
    }

    fn tick_2d(&mut self, dt: f32) {
        let n = self.positions.len();
        if n == 0 {
            return;
        }

        // Reset force accumulator.
        for f in self.forces.iter_mut() {
            *f = Vec2::ZERO;
        }

        // --- Repulsion via Barnes-Hut quadtree ----------------------
        self.build_tree();
        let theta2 = self.config.theta * self.config.theta;
        let repulsion = self.config.repulsion_strength;
        if !self.tree.is_empty() {
            for i in 0..n {
                let pos = self.positions[i];
                let f = accumulate_repulsion(&self.tree, 0, pos, i as u32, theta2, repulsion);
                self.forces[i] += f;
            }
        }

        // --- Link springs ------------------------------------------
        let link_k = self.config.link_strength;
        let link_d = self.config.link_distance;
        for e in &self.edges {
            let a = self.positions[e.src as usize];
            let b = self.positions[e.dst as usize];
            let delta = b - a;
            let dist = delta.length().max(1e-6);
            // Hooke's law: force magnitude is (dist - natural) * k * w,
            // along the edge direction. Positive (dist - link_d) pulls
            // endpoints together.
            let mag = (dist - link_d) * link_k * e.weight;
            let dir = delta / dist;
            self.forces[e.src as usize] += dir * mag;
            self.forces[e.dst as usize] -= dir * mag;
        }

        // --- Center gravity ----------------------------------------
        let c = self.config.center_strength;
        if c != 0.0 {
            for i in 0..n {
                self.forces[i] -= self.positions[i] * c;
            }
        }

        // --- Semi-implicit Euler integration ------------------------
        let decay = self.config.velocity_decay;
        for i in 0..n {
            self.velocities[i] *= decay;
            self.velocities[i] += self.forces[i] * dt;
            self.positions[i] += self.velocities[i] * dt;
        }
    }

    /// 3D counterpart of [`Self::tick_2d`]. Uses a Barnes-Hut octree
    /// (`force_3d.rs`) for repulsion and otherwise mirrors the 2D
    /// semi-implicit Euler integration, link springs, and center
    /// gravity. Also mirrors the x/y slice of each 3D position back
    /// into the 2D `positions` array so zero-copy consumers that
    /// only care about the projection still get updated values.
    fn tick_3d(&mut self, dt: f32) {
        let n = self.positions_3d.len();
        if n == 0 {
            return;
        }

        for f in self.forces_3d.iter_mut() {
            *f = Vec3::ZERO;
        }

        // --- Repulsion via Barnes-Hut octree -----------------------
        build_octree(&mut self.oct_tree, &self.positions_3d, MIN_CELL_SIZE);
        let theta2 = self.config.theta * self.config.theta;
        let repulsion = self.config.repulsion_strength;
        if !self.oct_tree.is_empty() {
            for i in 0..n {
                let pos = self.positions_3d[i];
                let f = accumulate_repulsion_3d(
                    &self.oct_tree,
                    0,
                    pos,
                    i as u32,
                    theta2,
                    repulsion,
                );
                self.forces_3d[i] += f;
            }
        }

        // --- Link springs ------------------------------------------
        let link_k = self.config.link_strength;
        let link_d = self.config.link_distance;
        for e in &self.edges {
            let a = self.positions_3d[e.src as usize];
            let b = self.positions_3d[e.dst as usize];
            let delta = b - a;
            let dist = delta.length().max(1e-6);
            let mag = (dist - link_d) * link_k * e.weight;
            let dir = delta / dist;
            self.forces_3d[e.src as usize] += dir * mag;
            self.forces_3d[e.dst as usize] -= dir * mag;
        }

        // --- Center gravity ----------------------------------------
        let c = self.config.center_strength;
        if c != 0.0 {
            for i in 0..n {
                self.forces_3d[i] -= self.positions_3d[i] * c;
            }
        }

        // --- Semi-implicit Euler integration ------------------------
        let decay = self.config.velocity_decay;
        for i in 0..n {
            self.velocities_3d[i] *= decay;
            self.velocities_3d[i] += self.forces_3d[i] * dt;
            self.positions_3d[i] += self.velocities_3d[i] * dt;
            // Mirror x/y into the 2D projection so callers going
            // through positions() still see live values.
            self.positions[i] = Vec2::new(self.positions_3d[i].x, self.positions_3d[i].y);
        }
    }

    /// Build a flat-arena Barnes-Hut quadtree covering all current
    /// positions. The tree is rebuilt every tick (positions drift by
    /// more than a cell between ticks so incremental maintenance
    /// wouldn't help).
    fn build_tree(&mut self) {
        self.tree.clear();
        let n = self.positions.len();
        if n == 0 {
            return;
        }

        let mut min = self.positions[0];
        let mut max = self.positions[0];
        for p in &self.positions[1..] {
            min = min.min(*p);
            max = max.max(*p);
        }
        // Square up the bounding box so quadrants stay square. Also
        // guard against degenerate (all-coincident) inputs.
        let extent = (max - min).max_element().max(MIN_CELL_SIZE);
        let center = (min + max) * 0.5;
        let half = Vec2::splat(extent * 0.5);
        let min = center - half;
        let max = center + half;

        self.tree.push(QuadNode::empty(min, max));
        for i in 0..n {
            let p = self.positions[i];
            insert(&mut self.tree, 0, i as u32, p);
        }

        // Finalise aggregate mass / centre-of-mass on internal nodes.
        finalise(&mut self.tree, 0);
    }
}

// ------------------------------------------------------------------
// Quadtree internals
// ------------------------------------------------------------------

/// A single quadtree node. Empty leaves have `body == None` and
/// `has_children == false`; one-body leaves have `body == Some(_)`
/// and `has_children == false`; internal nodes have `has_children
/// == true` and any mix of child slots populated. `mass` is the
/// aggregate body count of the subtree (unit mass per body) and
/// `com` the unweighted centroid, both filled in by `finalise`.
#[derive(Debug, Clone, Copy)]
struct QuadNode {
    min: Vec2,
    max: Vec2,
    mass: f32,
    com: Vec2,
    body: Option<u32>,
    children: [u32; 4], // u32::MAX sentinel for "no child"
    has_children: bool,
}

impl QuadNode {
    fn empty(min: Vec2, max: Vec2) -> Self {
        Self {
            min,
            max,
            mass: 0.0,
            com: Vec2::ZERO,
            body: None,
            children: [u32::MAX; 4],
            has_children: false,
        }
    }

    #[inline]
    fn size(&self) -> f32 {
        (self.max - self.min).max_element()
    }

    #[inline]
    fn center(&self) -> Vec2 {
        (self.min + self.max) * 0.5
    }

    /// Return 0..4 based on which quadrant `p` falls into relative to
    /// this node's centre. Bit 0 = east (x >= cx), bit 1 = north (y >= cy).
    #[inline]
    fn quadrant_of(&self, p: Vec2) -> usize {
        let c = self.center();
        let east = (p.x >= c.x) as usize;
        let north = (p.y >= c.y) as usize;
        (north << 1) | east
    }

    #[inline]
    fn child_bounds(&self, q: usize) -> (Vec2, Vec2) {
        let c = self.center();
        let (x_min, x_max) = if q & 1 == 1 {
            (c.x, self.max.x)
        } else {
            (self.min.x, c.x)
        };
        let (y_min, y_max) = if q & 2 == 2 {
            (c.y, self.max.y)
        } else {
            (self.min.y, c.y)
        };
        (Vec2::new(x_min, y_min), Vec2::new(x_max, y_max))
    }
}

/// Insert a body into the subtree rooted at `tree[idx]`.
///
/// Three cases:
///
/// 1. Empty leaf → drop the body in, record mass/com directly.
/// 2. One-body leaf with room to split → split into children and
///    reinsert both bodies.
/// 3. Internal node → recurse into the matching quadrant.
///
/// The "room to split" test is `size() > MIN_CELL_SIZE`. Below that
/// threshold we accumulate the new body into the leaf's mass/com
/// without splitting further — this avoids infinite recursion on
/// coincident positions. The accumulated leaf is still treated as
/// a single point in the force loop.
fn insert(tree: &mut Vec<QuadNode>, idx: usize, body: u32, pos: Vec2) {
    // Case 1: empty leaf.
    if !tree[idx].has_children && tree[idx].body.is_none() && tree[idx].mass == 0.0 {
        tree[idx].body = Some(body);
        tree[idx].com = pos;
        tree[idx].mass = 1.0;
        return;
    }

    // Case 2: one-body (or multi-body-accumulated) leaf. Split if we
    // still have room.
    if !tree[idx].has_children {
        if tree[idx].size() <= MIN_CELL_SIZE {
            // Too small to split further — accumulate. Weighted
            // average of existing COM and the new body.
            let new_mass = tree[idx].mass + 1.0;
            tree[idx].com = (tree[idx].com * tree[idx].mass + pos) / new_mass;
            tree[idx].mass = new_mass;
            // `body` stays as the first one inserted; the self-skip
            // in `accumulate_repulsion` will still skip that body
            // and the rest of the pile contributes force as part of
            // the aggregated leaf mass (close enough for physics).
            return;
        }
        let existing = tree[idx].body.take().expect("leaf without a body");
        let existing_pos = tree[idx].com;
        tree[idx].mass = 0.0;
        tree[idx].com = Vec2::ZERO;
        tree[idx].has_children = true;

        // Reinsert the two bodies into their quadrants. `quadrant_of`
        // on the same two distinct positions must eventually route
        // them into distinct cells, because we bail at MIN_CELL_SIZE
        // (case 2 above on the next level down).
        let q_existing = tree[idx].quadrant_of(existing_pos);
        let c_existing = create_or_get_child(tree, idx, q_existing);
        insert(tree, c_existing, existing, existing_pos);

        let q_new = tree[idx].quadrant_of(pos);
        let c_new = create_or_get_child(tree, idx, q_new);
        insert(tree, c_new, body, pos);
        return;
    }

    // Case 3: internal node — route into the matching quadrant.
    let q = tree[idx].quadrant_of(pos);
    let c = create_or_get_child(tree, idx, q);
    insert(tree, c, body, pos);
}

/// Return the existing child index at `quadrant` or allocate a new
/// child cell with the right sub-bounds.
fn create_or_get_child(tree: &mut Vec<QuadNode>, parent: usize, quadrant: usize) -> usize {
    let existing = tree[parent].children[quadrant];
    if existing != u32::MAX {
        return existing as usize;
    }
    let (cmin, cmax) = tree[parent].child_bounds(quadrant);
    let idx = tree.len();
    tree.push(QuadNode::empty(cmin, cmax));
    tree[parent].children[quadrant] = idx as u32;
    idx
}

/// Post-order pass: sum child masses and centres up into each
/// internal node. Leaves already have mass/com set during insertion.
fn finalise(tree: &mut [QuadNode], idx: usize) {
    if !tree[idx].has_children {
        return;
    }
    let children = tree[idx].children;
    let mut mass = 0.0;
    let mut com = Vec2::ZERO;
    for &c in &children {
        if c == u32::MAX {
            continue;
        }
        finalise(tree, c as usize);
        let child = tree[c as usize];
        mass += child.mass;
        com += child.com * child.mass;
    }
    if mass > 0.0 {
        com /= mass;
    }
    tree[idx].mass = mass;
    tree[idx].com = com;
}

/// Walk the Barnes-Hut tree and accumulate the repulsion force on
/// the node at `target_pos` (index `target_idx`, so we can skip
/// self-interactions at the leaf). `theta2` is the squared ratio
/// threshold for the approximation: when
/// `cell_size^2 < theta^2 * distance^2` the whole subtree is
/// collapsed to its centre-of-mass point charge.
fn accumulate_repulsion(
    tree: &[QuadNode],
    idx: usize,
    target_pos: Vec2,
    target_idx: u32,
    theta2: f32,
    repulsion: f32,
) -> Vec2 {
    let node = &tree[idx];
    if node.mass <= 0.0 {
        return Vec2::ZERO;
    }

    // Leaf handling: direct pairwise force, skip self. Multi-body
    // accumulated leaves still skip the stored `body` identity; the
    // remaining mass at that spot contributes a small residual
    // force at ~zero distance, which is capped by the softening
    // term (`+ 0.01`) in `pair_force`.
    if !node.has_children {
        if let Some(body) = node.body {
            if body == target_idx {
                // Subtract self: mass==1.0 for a single-body leaf,
                // > 1 for a pile-up. Everything beyond self still
                // contributes.
                let residual_mass = node.mass - 1.0;
                if residual_mass <= 0.0 {
                    return Vec2::ZERO;
                }
                return pair_force(target_pos, node.com, residual_mass, repulsion);
            }
        }
        return pair_force(target_pos, node.com, node.mass, repulsion);
    }

    let delta = node.com - target_pos;
    let dist2 = delta.length_squared();
    let size = node.size();
    // Barnes-Hut approximation check: treat subtree as a point when
    // it's far enough away relative to its size.
    if size * size < theta2 * dist2 {
        return pair_force(target_pos, node.com, node.mass, repulsion);
    }

    // Otherwise recurse into whichever children exist.
    let mut total = Vec2::ZERO;
    for &c in &node.children {
        if c == u32::MAX {
            continue;
        }
        total += accumulate_repulsion(tree, c as usize, target_pos, target_idx, theta2, repulsion);
    }
    total
}

/// Coulomb-like pairwise force on `from` due to a point mass `mass`
/// at `to`. With `repulsion < 0` (d3 default -30) the returned vector
/// points away from `to`; with `repulsion > 0` it points toward `to`.
///
/// We stay scalar here: for 2-float deltas there's nothing to
/// vectorise. The `graph_core::util::dot_simd` primitive is reserved
/// for the wide-float kernels and for a future batched center-gravity
/// reduction over all positions at once.
#[inline]
fn pair_force(from: Vec2, to: Vec2, mass: f32, repulsion: f32) -> Vec2 {
    let delta = to - from;
    let dist2 = delta.length_squared() + 0.01;
    let dist = dist2.sqrt();
    let dir = delta / dist;
    // Negative `repulsion` means the force on `from` is in the
    // `-dir` direction — away from `to`, which is the repulsion case.
    let magnitude = repulsion * mass / dist2;
    dir * magnitude
}

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

    #[test]
    fn two_linked_nodes_converge() {
        // `velocity_decay=0.4` (d3 default) is aggressive damping — it
        // multiplies velocity by 0.4 every tick, so with a small dt the
        // system is strongly overdamped and reaches equilibrium slowly.
        // Use a milder decay and bump the tick count so we end well
        // inside the 20% tolerance band around `link_distance`.
        let mut sim = Simulation::new(
            2,
            SimulationConfig {
                seed: 42,
                link_strength: 0.8,
                link_distance: 30.0,
                repulsion_strength: -30.0,
                center_strength: 0.0,
                velocity_decay: 0.9,
                ..Default::default()
            },
        );
        sim.set_edges(&[(0, 1, 1.0)]);
        for _ in 0..2000 {
            sim.tick(0.02);
        }
        let d = (sim.positions()[0] - sim.positions()[1]).length();
        let target = 30.0_f32;
        assert!(
            d > target * 0.8 && d < target * 1.2,
            "expected ~30 got {d}"
        );
    }

    #[test]
    fn star_topology_stabilizes() {
        let n = 9u32;
        let mut sim = Simulation::new(
            n,
            SimulationConfig {
                seed: 1,
                link_strength: 0.2,
                link_distance: 30.0,
                repulsion_strength: -30.0,
                center_strength: 0.1,
                velocity_decay: 0.5,
                ..Default::default()
            },
        );
        let edges: Vec<_> = (1..n).map(|i| (0u32, i, 1.0)).collect();
        sim.set_edges(&edges);
        for _ in 0..1000 {
            sim.tick(0.02);
        }
        let max_v = sim
            .velocities
            .iter()
            .map(|v| v.length())
            .fold(0.0_f32, f32::max);
        assert!(max_v < 0.5, "max_v={max_v}");
    }

    #[test]
    fn deterministic_tick_sequence() {
        let config = SimulationConfig {
            seed: 12345,
            ..Default::default()
        };
        let mut a = Simulation::new(50, config);
        let mut b = Simulation::new(50, config);
        let edges: Vec<_> = (0..49u32).map(|i| (i, i + 1, 1.0)).collect();
        a.set_edges(&edges);
        b.set_edges(&edges);
        for _ in 0..100 {
            a.tick(0.02);
            b.tick(0.02);
        }
        for (pa, pb) in a.positions().iter().zip(b.positions().iter()) {
            assert_eq!(pa, pb, "deterministic run diverged");
        }
    }

    #[test]
    fn isolated_nodes_repel() {
        let n = 10u32;
        let mut sim = Simulation::new(
            n,
            SimulationConfig {
                seed: 7,
                repulsion_strength: -80.0,
                center_strength: 0.0,
                velocity_decay: 0.6,
                ..Default::default()
            },
        );
        // Start clustered on a small 2D disc so repulsion has both
        // axes to work along. A perfectly collinear starting layout
        // stays collinear under symmetric repulsion, which makes the
        // per-neighbour spacing fall off too slowly for a clean
        // threshold test.
        let cluster_rng = &mut rand_chacha::ChaCha8Rng::seed_from_u64(99);
        for p in sim.positions.iter_mut() {
            loop {
                let x: f32 = cluster_rng.gen_range(-1.0..=1.0);
                let y: f32 = cluster_rng.gen_range(-1.0..=1.0);
                if x * x + y * y <= 1.0 {
                    *p = Vec2::new(x, y);
                    break;
                }
            }
        }
        for v in sim.velocities.iter_mut() {
            *v = Vec2::ZERO;
        }

        for _ in 0..500 {
            sim.tick(0.02);
        }

        let mut min_d = f32::INFINITY;
        for i in 0..sim.positions.len() {
            for j in (i + 1)..sim.positions.len() {
                let d = (sim.positions[i] - sim.positions[j]).length();
                if d < min_d {
                    min_d = d;
                }
            }
        }
        assert!(min_d > 5.0, "min pairwise distance {min_d} ≤ 5");
    }
}