graph_explorer_layout/

force_sim.rs

//! A persistent force simulation: node positions and velocities that survive
//! across frames, so a graph change nudges the layout instead of recomputing it.
//!
//! Parameter names mirror `d3-force` deliberately — anyone who has tuned a
//! force graph already knows what `alpha_decay` and `velocity_decay` do.
//!
//! Pure: no wasm, no GPU, and no time source. Ticks are driven by the caller,
//! never self-scheduled, which is what makes the whole thing host-testable.

use crate::QuadTree;
use graph_explorer_core::{Graph, NodeId};
use serde::Deserialize;
use std::collections::HashMap;

#[derive(Debug, Clone, Deserialize)]
#[serde(deny_unknown_fields)]
pub struct SimulationSpec {
    /// When false the simulation never ANIMATES — but a Force-source
    /// load/relayout still runs a synchronous settle so disabled hosts get a
    /// finished static layout (see `settle_forced`).
    #[serde(default = "yes")]
    pub enabled: bool,
    /// Per-tick multiplicative decay of `alpha`. The default reaches
    /// `alpha_min` in ~300 ticks: `1 - 0.001^(1/300)`.
    #[serde(default = "d_alpha_decay")]
    pub alpha_decay: f32,
    /// Below this the sim is considered settled and stops ticking, so a graph
    /// at rest costs nothing.
    #[serde(default = "d_alpha_min")]
    pub alpha_min: f32,
    /// Friction: velocity is multiplied by `1 - velocity_decay` each tick.
    #[serde(default = "d_velocity_decay")]
    pub velocity_decay: f32,
    /// `alpha` restored on a graph change. Lower than 1 so a reveal nudges
    /// rather than re-scrambling.
    #[serde(default = "d_reheat_alpha")]
    pub reheat_alpha: f32,
    /// Barnes-Hut accuracy. 0 is exact and `O(n²)`; larger approximates more.
    #[serde(default = "d_theta")]
    pub theta: f32,
    /// Ideal edge length, in world units.
    #[serde(default = "d_link_distance")]
    pub link_distance: f32,
    /// Many-body strength. Negative repels.
    #[serde(default = "d_charge")]
    pub charge_strength: f32,
    /// Pull toward the centroid. Without it, disconnected components drift
    /// apart forever — the old (now-deleted) one-shot layout had no such force.
    #[serde(default = "d_center")]
    pub center_strength: f32,
}

fn yes() -> bool { true }
fn d_alpha_decay() -> f32 { 0.0228 }
fn d_alpha_min() -> f32 { 0.001 }
fn d_velocity_decay() -> f32 { 0.4 }
fn d_reheat_alpha() -> f32 { 0.3 }
fn d_theta() -> f32 { 0.9 }
fn d_link_distance() -> f32 { 120.0 }
fn d_charge() -> f32 { -800.0 }
fn d_center() -> f32 { 0.02 }

impl Default for SimulationSpec {
    fn default() -> Self {
        Self {
            enabled: yes(),
            alpha_decay: d_alpha_decay(),
            alpha_min: d_alpha_min(),
            velocity_decay: d_velocity_decay(),
            reheat_alpha: d_reheat_alpha(),
            theta: d_theta(),
            link_distance: d_link_distance(),
            charge_strength: d_charge(),
            center_strength: d_center(),
        }
    }
}

pub struct ForceSimulation {
    spec: SimulationSpec,
    index: HashMap<NodeId, usize>,
    ids: Vec<NodeId>,
    pos: Vec<[f32; 2]>,
    vel: Vec<[f32; 2]>,
    links: Vec<(usize, usize)>,
    /// At most one node is held in place — the one the user is standing on.
    fixed: Option<usize>,
    /// Where the pinned node is held. Captured when the pin is set, so the
    /// node returns to exactly that point rather than to wherever the last
    /// integration happened to leave it.
    pinned_at: Option<[f32; 2]>,
    alpha: f32,
    /// Monotonic count of nodes that have entered without a seed. Persisting it
    /// on the struct — rather than restarting per `sync` — is what stops two
    /// arrivals in *different* syncs from being handed the same spiral point
    /// and stacking, which is the exact failure the spiral exists to prevent.
    spawn_seq: usize,
}

/// Upper bound for `theta`. Beyond roughly 1.25 — and as low as 0.75 for a
/// tight cluster with a distant outlier, which is precisely the shape an
/// expanded aggregate makes — a cell that *contains* body `i` can pass the
/// approximation test, and `i` then repels itself through its own cell's mass.
/// Error stays within 21% at 1.0, so that is where the config boundary sits.
/// Clamping here rather than in `walk` keeps the check off the hot path.
const THETA_MAX: f32 = 1.0;

/// Floor for `alpha_min`. `is_active` is `alpha >= alpha_min`, and `alpha`
/// decays multiplicatively toward zero without ever reaching it — so at
/// `alpha_min = 0` the predicate is permanently true and the simulation never
/// settles. That is not a slow layout, it is three broken guarantees: the
/// host's render loop runs a Barnes-Hut tick every frame at rest, the
/// settle-detecting `!tick()` edge never fires (so a caller that rebuilds its
/// scene there never does, and its target goes stale indefinitely), and
/// `settle(max)` burns its entire budget and still reports active, which is
/// what a reduced-motion caller uses to mean "snap, do not animate".
///
/// `1e-4` is the smallest value that still terminates promptly: from a full
/// reheat at the default decay it is reached in ~400 ticks.
const ALPHA_MIN_FLOOR: f32 = 1e-4;

/// Ceiling for `alpha_min`, guarding the opposite failure. `alpha_min` above
/// the alpha a reheat restores means `is_active` is false the moment it is
/// set, so every reheat is silently a no-op and the simulation appears dead.
/// `0.1` leaves headroom under the default `reheat_alpha` of 0.3. A host that
/// also lowers `reheat_alpha` below its `alpha_min` can still arrange that
/// no-op, but it fails safe — no motion, no CPU — rather than wedging.
const ALPHA_MIN_CEIL: f32 = 0.1;

/// Floor for `alpha_decay`. The schedule is `alpha -= alpha * alpha_decay`, so
/// at zero `alpha` never decreases: the same never-settles failure as
/// `alpha_min = 0`, reached from the other side.
///
/// The value is derived, not chosen. `settle(600)` is the reduced-motion
/// contract — settle synchronously, then snap, with no motion afterwards — and
/// it only holds if 600 ticks is actually enough. Worst case is a full reheat
/// (`alpha = 1.0`, the ceiling `reheat_alpha` is clamped to) decaying to the
/// smallest legal `alpha_min`, `ALPHA_MIN_FLOOR`, which takes
/// `ln(1e-4) / ln(1 - alpha_decay)` ticks: 571 at `0.016`, against a budget of
/// 600. Anything slower silently converts reduced motion into permanent
/// motion. The price is that a host cannot ask for a settle longer than about
/// 9.5 seconds at 60fps, which is a fair ceiling on how long a graph may wobble.
const ALPHA_DECAY_MIN: f32 = 0.016;

/// Ceiling for `alpha_decay`. At 1.0 `alpha` drops to zero in a single tick,
/// which is already the fastest meaningful schedule; above it `alpha` goes
/// negative, which is not a faster settle, just a nonsensical state.
const ALPHA_DECAY_MAX: f32 = 1.0;

/// Floor for `velocity_decay`, which enters the integrator as
/// `friction = 1 - velocity_decay`. At zero there is no damping at all: a node
/// keeps whatever velocity it accumulated and coasts in a straight line until
/// `alpha` decays out, freezing it wherever it happened to be rather than at a
/// force fixpoint. Below zero `friction > 1` and this explicit-Euler step
/// *amplifies* velocity every tick — exponential divergence, NaN within the
/// usual few hundred ticks. `0.01` guarantees energy actually leaves the
/// system.
const VELOCITY_DECAY_MIN: f32 = 0.01;

/// Ceiling for `velocity_decay`. 1.0 is full damping — velocity is discarded
/// each tick and the step degenerates to plain gradient descent, which is
/// stable and occasionally what a host wants. Above it `friction` is negative
/// and velocity flips sign every tick, so the layout oscillates instead of
/// converging (and diverges outright past 2.0).
const VELOCITY_DECAY_MAX: f32 = 1.0;

/// Radius scale for an unseeded spawn, in world units.
const SPAWN_RADIUS: f32 = 30.0;

/// How many spiral rings a *late* arrival may use before the radius wraps.
/// `SPAWN_RADIUS * sqrt(11.5) ≈ 102`, so an arrival lands within about one
/// link distance of the centroid instead of at `30·√(global rank)` — which for
/// a 500-node graph was 671 units, 5.6 link distances, and nowhere near the
/// "near the centroid" the doc contract promises. The angle keeps advancing
/// off the unwrapped sequence, so wrapping the radius never repeats a point.
const ARRIVAL_RINGS: usize = 12;

impl ForceSimulation {
    pub fn new(spec: SimulationSpec) -> Self {
        Self {
            spec: Self::clamped(spec),
            index: HashMap::new(),
            ids: Vec::new(),
            pos: Vec::new(),
            vel: Vec::new(),
            links: Vec::new(),
            fixed: None,
            pinned_at: None,
            alpha: 0.0,
            spawn_seq: 0,
        }
    }

    pub fn set_spec(&mut self, spec: SimulationSpec) { self.spec = Self::clamped(spec); }
    pub fn spec(&self) -> &SimulationSpec { &self.spec }

    /// Bring a host-supplied spec inside the range the forces can actually
    /// integrate, and the alpha schedule inside the range that can actually
    /// terminate. Done once, where the spec is adopted, so no per-tick branch
    /// pays for it — and here rather than at the call site because this is the
    /// config boundary: `new` and `set_spec` are the only two ways a spec gets
    /// in, and `set_spec` is reachable from a host in a single call.
    ///
    /// The alpha-schedule bounds are not taste. `is_active` is the sole thing
    /// standing between a settled graph and a render loop that ticks forever,
    /// and it is one comparison against two of these values; see the constants
    /// for what each unbounded end breaks.
    fn clamped(mut spec: SimulationSpec) -> SimulationSpec {
        spec.theta = spec.theta.clamp(0.0, THETA_MAX);
        spec.alpha_min = spec.alpha_min.clamp(ALPHA_MIN_FLOOR, ALPHA_MIN_CEIL);
        spec.alpha_decay = spec.alpha_decay.clamp(ALPHA_DECAY_MIN, ALPHA_DECAY_MAX);
        spec.velocity_decay =
            spec.velocity_decay.clamp(VELOCITY_DECAY_MIN, VELOCITY_DECAY_MAX);
        // Nothing above 1 is reachable by `reheat_full`, and a negative alpha is
        // just an immediately-settled one wearing a disguise. Bounding it keeps
        // `alpha` inside `[0, 1]` for every code path, which is what makes the
        // tick-count arithmetic in the constants above true.
        spec.reheat_alpha = spec.reheat_alpha.clamp(0.0, 1.0);
        spec
    }

    /// Adopt the graph's topology, preserving position and velocity for every
    /// id already present. New ids take `seed` when given — a revealed ghost
    /// enters exactly where it was drawn, so physics continues from the
    /// preview rather than teleporting. Does not reheat; callers decide that.
    pub fn sync(&mut self, graph: &Graph, seed: &HashMap<NodeId, [f32; 2]>) {
        let mut ids: Vec<NodeId> = graph.nodes().map(|n| n.id.clone()).collect();
        ids.sort(); // deterministic index assignment, so layouts reproduce

        let mut pos = Vec::with_capacity(ids.len());
        let mut vel = Vec::with_capacity(ids.len());
        let mut index = HashMap::with_capacity(ids.len());
        let centroid = self.centroid();
        // A cold start is the one case where spiral rank *is* arrival count, so
        // `30·√N` is the right spread and is kept exactly. Every later sync is
        // an arrival into an existing layout, where that radius would be a
        // function of where the id happens to sort rather than of anything
        // physical, so those wrap instead.
        let cold_start = self.index.is_empty();
        let mut spawned_here = 0usize;

        for (i, id) in ids.iter().enumerate() {
            match self.index.get(id) {
                Some(&old) => { pos.push(self.pos[old]); vel.push(self.vel[old]); }
                None => {
                    match seed.get(id) {
                        Some(&at) => pos.push(at),
                        None => {
                            let ring = if cold_start {
                                spawned_here
                            } else {
                                self.spawn_seq % ARRIVAL_RINGS
                            };
                            pos.push(Self::phyllotaxis(self.spawn_seq, ring, centroid));
                            self.spawn_seq += 1;
                            spawned_here += 1;
                        }
                    }
                    vel.push([0.0, 0.0]);
                }
            }
            index.insert(id.clone(), i);
        }

        self.links = graph
            .edges()
            .filter_map(|e| Some((*index.get(&e.source)?, *index.get(&e.target)?)))
            .collect();
        // `Graph::edges()` iterates a HashMap, so its order differs between two
        // independently built copies of the same graph. Float addition is not
        // associative, so an unsorted `links` makes the spring accumulation —
        // and therefore the whole settled layout — vary run to run. Sorting
        // costs one pass per topology change and restores the reproducibility
        // the sorted `ids` above was already reaching for. Equal pairs are
        // interchangeable (identical contributions), so a stable tie is fine.
        self.links.sort_unstable();

        // Re-resolve the pin through the id, since indices just changed.
        let pinned_id = self.fixed.and_then(|f| self.ids.get(f)).cloned();
        self.ids = ids;
        self.pos = pos;
        self.vel = vel;
        self.index = index;
        self.fixed = pinned_id.and_then(|id| self.index.get(&id).copied());
        self.pinned_at = self.fixed.map(|i| self.pos[i]);
    }

    /// Deterministic spread for a node with no seed — a sunflower spiral, so
    /// unseeded arrivals never stack on one point.
    ///
    /// `seq` drives the angle and `ring` the radius. On a cold start the two
    /// are equal and this is the classic `30·√i` sunflower. For a later arrival
    /// `ring` wraps (see [`ARRIVAL_RINGS`]) so the radius stays near the
    /// centroid, while `seq` keeps climbing — the golden angle is an irrational
    /// fraction of a turn, so a never-repeating angle is on its own enough to
    /// guarantee two spawns never coincide.
    fn phyllotaxis(seq: usize, ring: usize, center: [f32; 2]) -> [f32; 2] {
        let r = SPAWN_RADIUS * (ring as f32 + 0.5).sqrt();
        let a = seq as f32 * 2.399_963; // golden angle
        [center[0] + r * a.cos(), center[1] + r * a.sin()]
    }

    fn centroid(&self) -> [f32; 2] {
        if self.pos.is_empty() { return [0.0, 0.0]; }
        let n = self.pos.len() as f32;
        let s = self.pos.iter().fold([0.0f32, 0.0], |a, p| [a[0] + p[0], a[1] + p[1]]);
        [s[0] / n, s[1] / n]
    }

    /// Re-energise after a graph change. Sets `alpha` unconditionally —
    /// `enabled` is enforced solely by `is_active`/`tick`, so a disabled sim
    /// still accumulates the alpha a later `settle_forced` needs to have
    /// anything to act on (reduced-motion and physics-disabled hosts reheat
    /// exactly like a live one, then force a settle instead of ticking).
    pub fn reheat(&mut self) {
        self.alpha = self.spec.reheat_alpha;
    }

    /// Re-energise from scratch (`alpha = 1`), for a cold start. See `reheat`
    /// on why this does not gate on `enabled`.
    pub fn reheat_full(&mut self) {
        self.alpha = 1.0;
    }

    /// Hold `id` in place. `None` releases. Unknown ids are ignored.
    pub fn pin(&mut self, id: Option<&str>) {
        self.fixed = id.and_then(|i| self.index.get(i).copied());
        self.pinned_at = self.fixed.map(|i| self.pos[i]);
    }

    pub fn is_active(&self) -> bool { self.spec.enabled && self.alpha >= self.spec.alpha_min }

    /// Whether `alpha` is still above the settle threshold — `is_active`
    /// without the `enabled` term, so a disabled sim can be asked whether it
    /// was frozen mid-settle (and so needs a `settle_forced` to finish).
    pub fn is_unsettled(&self) -> bool { self.alpha >= self.spec.alpha_min }

    pub fn position(&self, id: &str) -> Option<[f32; 2]> {
        self.index.get(id).map(|&i| self.pos[i])
    }

    pub fn velocity(&self, id: &str) -> Option<[f32; 2]> {
        self.index.get(id).map(|&i| self.vel[i])
    }

    pub fn positions(&self) -> impl Iterator<Item = (&NodeId, [f32; 2])> + '_ {
        self.ids.iter().zip(self.pos.iter().copied())
    }

    /// Tick to convergence, bounded by `max` so a pathological graph cannot
    /// hang the caller. Used by reduced-motion, which snaps rather than
    /// animating.
    pub fn settle(&mut self, max: usize) {
        for _ in 0..max {
            if !self.tick() { return; }
        }
    }

    /// `settle`, but runs even when `spec.enabled` is false — for hosts that
    /// need a finished layout without any animation (reduced motion, or
    /// physics disabled while the Force layout source is active). The flag is
    /// restored unconditionally before returning; `tick` runs no host code,
    /// so a scoped set/restore is panic-safe in practice.
    pub fn settle_forced(&mut self, max: usize) {
        let was = self.spec.enabled;
        self.spec.enabled = true;
        self.settle(max);
        self.spec.enabled = was;
    }

    /// Advance one step. Returns false once settled (`alpha < alpha_min`), so
    /// a caller can stop scheduling frames when the graph is at rest.
    pub fn tick(&mut self) -> bool {
        if !self.is_active() {
            return false;
        }
        let n = self.pos.len();
        if n == 0 {
            self.alpha = 0.0;
            return false;
        }

        let mut force = vec![[0.0f32, 0.0]; n];

        // Many-body repulsion, Barnes-Hut approximated. Rebuilt each tick
        // because every position moved.
        let tree = QuadTree::build(&self.pos);
        for (i, f) in force.iter_mut().enumerate() {
            let r = tree.repulsion(i, &self.pos, self.spec.theta, self.spec.charge_strength);
            f[0] += r[0];
            f[1] += r[1];
        }

        // Link springs: pull each edge toward `link_distance`.
        //
        // Both the stiffness and the split are normalised by degree, exactly as
        // d3-force does, and it is load-bearing rather than cosmetic. A flat
        // 0.5/0.5 split gives a degree-k node a summed stiffness of k/2, and
        // this integrator is explicit Euler: once `friction * k/2 > 2` the node
        // overshoots further every tick and the layout diverges. Measured with
        // a flat split: stable to degree 18, exponential blow-up by 20, NaN at
        // 26 — and the demo's seed graph has a node with 26 neighbours.
        // `1/min(deg)` bounds each node's total stiffness at ~1 instead, which
        // is stable for any degree.
        let mut deg = vec![0u32; n];
        for &(a, b) in &self.links {
            deg[a] += 1;
            deg[b] += 1;
        }
        for &(a, b) in &self.links {
            let dx = self.pos[b][0] - self.pos[a][0];
            let dy = self.pos[b][1] - self.pos[a][1];
            let d = (dx * dx + dy * dy).sqrt().max(1e-3);
            let strength = 1.0 / deg[a].min(deg[b]).max(1) as f32;
            // Positive when too far apart, so both ends move toward each other.
            let pull = (d - self.spec.link_distance) / d * strength;
            // Share of the correction taken by `b`. The busier end moves less,
            // so a hub stays put and its leaves swing around it.
            let (da, db) = (deg[a] as f32, deg[b] as f32);
            let share_b = da / (da + db);
            force[a][0] += dx * pull * (1.0 - share_b);
            force[a][1] += dy * pull * (1.0 - share_b);
            force[b][0] -= dx * pull * share_b;
            force[b][1] -= dy * pull * share_b;
        }

        // Weak centering, so disconnected components stay in frame.
        let c = self.centroid();
        for (i, f) in force.iter_mut().enumerate() {
            f[0] += (c[0] - self.pos[i][0]) * self.spec.center_strength;
            f[1] += (c[1] - self.pos[i][1]) * self.spec.center_strength;
        }

        // Integrate. Force is scaled by alpha so motion damps out with the
        // schedule, and velocity is damped by friction — d3's Verlet-ish step.
        let friction = 1.0 - self.spec.velocity_decay;
        let alpha = self.alpha;
        for ((p, v), f) in self.pos.iter_mut().zip(self.vel.iter_mut()).zip(force.iter()) {
            v[0] = (v[0] + f[0] * alpha) * friction;
            v[1] = (v[1] + f[1] * alpha) * friction;
            p[0] += v[0];
            p[1] += v[1];
        }

        // The pinned node exerts force but never moves. Restoring position (not
        // just zeroing velocity) is required: `sync` recaptures `pinned_at`
        // from the current position, so drift would ratchet the hold point.
        if let Some(f) = self.fixed {
            if let Some(at) = self.pinned_at {
                self.pos[f] = at;
            }
            self.vel[f] = [0.0, 0.0];
        }

        self.alpha -= self.alpha * self.spec.alpha_decay;
        self.is_active()
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use graph_explorer_core::{Edge, Graph, Node};
    use std::collections::HashMap;

    fn graph_of(nodes: &[&str], edges: &[(&str, &str, &str)]) -> Graph {
        let mut g = Graph::default();
        for n in nodes { g.add_node(Node::new(*n)); }
        for (id, a, b) in edges { g.add_edge(Edge::new(*id, *a, *b)); }
        g
    }

    fn seeded(pairs: &[(&str, [f32; 2])]) -> HashMap<String, [f32; 2]> {
        pairs.iter().map(|(k, v)| ((*k).to_string(), *v)).collect()
    }

    #[test]
    fn is_unsettled_tracks_alpha_not_enabled() {
        let mut sim = ForceSimulation::new(SimulationSpec::default());
        sim.sync(&graph_of(&["a", "b"], &[("e", "a", "b")]), &seeded(&[]));
        sim.reheat_full();
        assert!(sim.is_unsettled());
        sim.settle(10_000);
        assert!(!sim.is_unsettled());
    }

    #[test]
    fn defaults_match_the_documented_schedule() {
        let s = SimulationSpec::default();
        assert!(s.enabled);
        assert_eq!(s.alpha_min, 0.001);
        assert_eq!(s.reheat_alpha, 0.3);
        assert_eq!(s.theta, 0.9);
        // ~300 ticks from alpha 1 to alpha_min, which is d3's schedule.
        assert!((s.alpha_decay - 0.0228).abs() < 1e-4);
    }

    #[test]
    fn spec_rejects_unknown_fields() {
        // `deny_unknown_fields` is what turns a host's typo into an error
        // instead of a silently ignored setting.
        let r: Result<SimulationSpec, _> = serde_json::from_str(r#"{"charge": -30}"#);
        assert!(r.is_err(), "unknown field must be rejected");
    }

    #[test]
    fn alpha_decays_to_min_in_about_three_hundred_ticks_then_stops() {
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&["a", "b"], &[("e", "a", "b")]), &HashMap::new());
        s.reheat_full();
        let mut ticks = 0;
        while s.tick() { ticks += 1; assert!(ticks < 1000, "runaway"); }
        assert!((250..=350).contains(&ticks), "settled in {ticks} ticks, expected ~300");
        assert!(!s.is_active(), "a settled sim must report inactive");
        assert!(!s.tick(), "ticking a settled sim stays settled");
    }

    #[test]
    fn reheat_reenergises_a_settled_sim() {
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&["a", "b"], &[("e", "a", "b")]), &HashMap::new());
        while s.tick() {}
        assert!(!s.is_active());
        s.reheat();
        assert!(s.is_active(), "reheat must restart the schedule");
    }

    #[test]
    fn sync_preserves_position_and_velocity_for_known_nodes() {
        // Warm start is the whole reason a one-shot layout could not be reused:
        // without it, every reveal teleports the graph.
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&["a", "b"], &[("e", "a", "b")]), &seeded(&[("a", [100.0, 0.0]), ("b", [-100.0, 0.0])]));
        s.reheat();
        s.tick();
        let a_before = s.position("a").unwrap();
        let v_before = s.velocity("a").unwrap();
        assert!(v_before[0] != 0.0 || v_before[1] != 0.0, "expected motion to compare against");

        // A third node arrives — a reveal.
        s.sync(&graph_of(&["a", "b", "c"], &[("e", "a", "b")]), &seeded(&[("c", [0.0, 50.0])]));
        assert_eq!(s.position("a").unwrap(), a_before, "existing position preserved");
        assert_eq!(s.velocity("a").unwrap(), v_before, "existing velocity preserved");
        assert_eq!(s.position("c").unwrap(), [0.0, 50.0], "new node enters at its seed");
    }

    #[test]
    fn sync_drops_nodes_that_left_the_graph() {
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&["a", "b"], &[]), &HashMap::new());
        assert!(s.position("b").is_some());
        s.sync(&graph_of(&["a"], &[]), &HashMap::new());
        assert!(s.position("b").is_none(), "a removed node must leave the sim");
    }

    #[test]
    fn an_unseeded_new_node_still_gets_a_finite_distinct_position() {
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&["a", "b", "c"], &[]), &HashMap::new());
        let ps: Vec<[f32; 2]> = ["a", "b", "c"].iter().map(|i| s.position(i).unwrap()).collect();
        for p in &ps { assert!(p[0].is_finite() && p[1].is_finite()); }
        assert!(ps[0] != ps[1] && ps[1] != ps[2], "unseeded nodes must not stack");
    }

    #[test]
    fn a_pinned_node_holds_still_while_its_neighbours_move() {
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&["a", "b"], &[("e", "a", "b")]), &seeded(&[("a", [0.0, 0.0]), ("b", [400.0, 0.0])]));
        s.pin(Some("a"));
        s.reheat_full();
        for _ in 0..30 { s.tick(); }
        assert_eq!(s.position("a").unwrap(), [0.0, 0.0], "pinned node must not drift");
        assert!(s.position("b").unwrap()[0] < 400.0, "its neighbour should have been pulled in");
    }

    #[test]
    fn a_pinned_node_still_exerts_force() {
        // Pinning fixes position, not influence — otherwise the graph would
        // collapse through the node the user is standing on.
        //
        // `b` starts *off* both axes so this is a test about force rather than
        // about motion. Asserting only that `b` moved, or that one coordinate
        // grew, is satisfied by a model that slides the whole graph a constant
        // step each tick with no repulsion at all. Repulsion from a fixed `a`
        // has a signature that drift does not: it is purely radial, so `b` must
        // travel straight out along the a->b ray and stay on it.
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&["a", "b"], &[]), &seeded(&[("a", [0.0, 0.0]), ("b", [3.0, 4.0])]));
        s.pin(Some("a"));
        s.reheat_full();
        for _ in 0..20 { s.tick(); }

        let (a, b) = (s.position("a").unwrap(), s.position("b").unwrap());
        assert_eq!(a, [0.0, 0.0], "the pinned node must not have moved");

        let d = ((b[0] - a[0]).powi(2) + (b[1] - a[1]).powi(2)).sqrt();
        assert!(d > 5.0, "b sits {d} from a, no further than it started; nothing repelled it");

        // `a` is at the origin, so the ray is just `b` normalised: (0.6, 0.8).
        let unit = [(b[0] - a[0]) / d, (b[1] - a[1]) / d];
        assert!(
            (unit[0] - 0.6).abs() < 0.01 && (unit[1] - 0.8).abs() < 0.01,
            "b reached {b:?}, off the a->b ray (direction {unit:?}); that is drift, not repulsion"
        );
    }

    /// `sync` re-resolves the pin through the node id, because sorting the ids
    /// reshuffles every index as soon as one is added or removed. That is the
    /// subtlest logic in the warm start and it fails *silently* — a wrong index
    /// holds the wrong node rather than panicking — so it needs its own test.
    #[test]
    fn a_pin_survives_the_index_shift_from_a_new_node() {
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&["m", "z"], &[]), &seeded(&[("m", [10.0, 0.0]), ("z", [200.0, 0.0])]));
        s.pin(Some("m")); // m is index 0 of [m, z]

        // "a" sorts before "m", so m slides from index 0 to index 1.
        s.sync(&graph_of(&["a", "m", "z"], &[]), &seeded(&[("a", [0.0, 300.0])]));
        s.reheat_full();
        for _ in 0..60 { s.tick(); }

        assert_eq!(
            s.position("m").unwrap(), [10.0, 0.0],
            "the pin must still hold m, at m's own hold point, after the reindex"
        );
        assert!(s.position("z").unwrap() != [200.0, 0.0], "z was not pinned and should have moved");
        assert!(s.position("a").unwrap() != [0.0, 300.0], "a was not pinned and should have moved");
    }

    /// A pinned node that leaves the graph must drop the pin, not hand it to
    /// whichever survivor inherits its index — that would freeze an innocent
    /// node the user never pinned.
    #[test]
    fn a_pinned_node_leaving_the_graph_drops_the_pin() {
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(
            &graph_of(&["m", "q", "z"], &[]),
            &seeded(&[("m", [0.0, 0.0]), ("q", [40.0, 0.0]), ("z", [0.0, 40.0])]),
        );
        s.pin(Some("m")); // index 0 of [m, q, z]

        // m leaves; "q" inherits index 0. A pin tracked by index would hold q.
        s.sync(&graph_of(&["q", "z"], &[]), &HashMap::new());
        s.reheat_full();
        for _ in 0..60 { s.tick(); }

        assert!(s.position("m").is_none(), "m left the graph");
        assert!(
            s.position("q").unwrap() != [40.0, 0.0],
            "q inherited m's index and is frozen; the pin transferred instead of dropping"
        );
        assert!(s.position("z").unwrap() != [0.0, 40.0], "z should be free to move");
    }

    #[test]
    fn unpinning_releases_the_node() {
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&["a", "b"], &[]), &seeded(&[("a", [0.0, 0.0]), ("b", [5.0, 0.0])]));
        s.pin(Some("a"));
        s.reheat_full();
        for _ in 0..10 { s.tick(); }
        s.pin(None);
        s.reheat_full();
        for _ in 0..20 { s.tick(); }
        assert!(s.position("a").unwrap() != [0.0, 0.0], "released node should move");
    }

    #[test]
    fn pinning_an_unknown_id_is_a_no_op_not_a_panic() {
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&["a"], &[]), &HashMap::new());
        s.pin(Some("ghost-that-is-not-here"));
        s.reheat_full();
        s.tick();
        assert!(s.position("a").unwrap()[0].is_finite());
    }

    #[test]
    fn disabled_spec_never_activates() {
        let spec = SimulationSpec { enabled: false, ..SimulationSpec::default() };
        let mut s = ForceSimulation::new(spec);
        s.sync(&graph_of(&["a", "b"], &[]), &HashMap::new());
        s.reheat();
        assert!(!s.is_active(), "a disabled sim must never run");
        assert!(!s.tick());
    }

    #[test]
    fn settle_converges_and_respects_its_bound() {
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&["a", "b", "c"], &[("e1", "a", "b"), ("e2", "b", "c")]), &HashMap::new());
        s.reheat_full();
        s.settle(1000);
        assert!(!s.is_active(), "settle should reach convergence");

        s.reheat_full();
        s.settle(5);
        assert!(s.is_active(), "a bounded settle must stop early rather than spin");
    }

    #[test]
    fn settle_forced_runs_while_disabled_and_restores_the_flag() {
        let spec = SimulationSpec { enabled: false, ..SimulationSpec::default() };
        let mut s = ForceSimulation::new(spec);
        s.sync(&graph_of(&["a", "b", "c"], &[("e1", "a", "b"), ("e2", "b", "c")]), &HashMap::new());
        s.reheat_full();
        let before: Vec<[f32; 2]> = s.positions().map(|(_, p)| p).collect();
        s.settle_forced(600);
        assert!(!s.spec().enabled, "enabled flag must be restored");
        assert!(!s.is_active(), "disabled sim reports inactive after the forced settle");
        let after: Vec<[f32; 2]> = s.positions().map(|(_, p)| p).collect();
        assert_ne!(before, after, "forced settle must actually move nodes");
    }

    #[test]
    fn settle_forced_on_an_enabled_sim_behaves_like_settle() {
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&["a", "b"], &[("e1", "a", "b")]), &HashMap::new());
        s.reheat_full();
        s.settle_forced(1000);
        assert!(s.spec().enabled);
        assert!(!s.is_active(), "converged");
    }

    #[test]
    fn a_triangle_settles_to_roughly_equilateral() {
        let mut s = ForceSimulation::new(SimulationSpec::default());
        let g = graph_of(&["a", "b", "c"], &[("e1", "a", "b"), ("e2", "b", "c"), ("e3", "c", "a")]);
        s.sync(&g, &seeded(&[("a", [0.0, 0.0]), ("b", [10.0, 0.0]), ("c", [0.0, 10.0])]));
        s.reheat_full();
        s.settle(2000);

        let d = |x: &str, y: &str| {
            let (p, q) = (s.position(x).unwrap(), s.position(y).unwrap());
            ((p[0] - q[0]).powi(2) + (p[1] - q[1]).powi(2)).sqrt()
        };
        let (ab, bc, ca) = (d("a", "b"), d("b", "c"), d("c", "a"));
        let mean = (ab + bc + ca) / 3.0;
        for (name, side) in [("ab", ab), ("bc", bc), ("ca", ca)] {
            assert!((side - mean).abs() / mean < 0.15, "{name}={side} vs mean {mean}: not equilateral");
        }
        assert!(mean > 20.0, "sides collapsed to {mean}; repulsion is not holding them apart");
    }

    #[test]
    fn linked_nodes_are_pulled_toward_link_distance() {
        let mut s = ForceSimulation::new(SimulationSpec::default());
        // Start far apart; the spring should reel them in.
        s.sync(&graph_of(&["a", "b"], &[("e", "a", "b")]), &seeded(&[("a", [0.0, 0.0]), ("b", [900.0, 0.0])]));
        s.reheat_full();
        s.settle(2000);
        let (p, q) = (s.position("a").unwrap(), s.position("b").unwrap());
        let d = ((p[0] - q[0]).powi(2) + (p[1] - q[1]).powi(2)).sqrt();
        assert!(d < 400.0, "linked nodes stayed {d} apart; the spring is not pulling");
    }

    #[test]
    fn a_disconnected_component_does_not_drift_to_infinity() {
        // Pure repulsion pushes unlinked components apart forever. The
        // centering force is what bounds them, and nothing else tests it.
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&["a", "b", "c", "d"], &[("e", "a", "b")]), &HashMap::new());
        s.reheat_full();
        s.settle(3000);
        for id in ["a", "b", "c", "d"] {
            let p = s.position(id).unwrap();
            assert!(p[0].is_finite() && p[1].is_finite(), "{id} went non-finite");
            assert!(p[0].abs() < 5000.0 && p[1].abs() < 5000.0, "{id} drifted to {p:?}");
        }
    }

    /// `a_disconnected_component_does_not_drift_to_infinity` does not actually
    /// test the centering force: with `center_strength = 0` its worst
    /// coordinate is 495, well inside the 5000 bound, so it passes either way.
    /// One alpha schedule is simply too short for the drift to show. Repeated
    /// reheats are what separate the two: unbounded growth without centering
    /// (measured 504 -> 972 -> 1269 -> 1504 over ten reheats, still climbing)
    /// versus a fixed point with it (291 -> 269 -> 265 -> 265).
    #[test]
    fn centering_bounds_spread_across_repeated_reheats() {
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&["a", "b", "c", "d"], &[("e", "a", "b")]), &HashMap::new());
        for _ in 0..10 {
            s.reheat_full();
            s.settle(3000);
        }
        let c = s.centroid();
        let radius = s
            .positions()
            .map(|(_, p)| ((p[0] - c[0]).powi(2) + (p[1] - c[1]).powi(2)).sqrt())
            .fold(0.0f32, f32::max);
        assert!(
            radius < 600.0,
            "spread reached {radius} after ten reheats; the centering force is not bounding it"
        );
    }

    #[test]
    fn an_isolated_single_node_does_not_move() {
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&["only"], &[]), &seeded(&[("only", [42.0, -7.0])]));
        s.reheat_full();
        s.settle(500);
        let p = s.position("only").unwrap();
        assert!((p[0] - 42.0).abs() < 0.01 && (p[1] + 7.0).abs() < 0.01,
            "a lone node has nothing to push against and should sit still, got {p:?}");
    }

    #[test]
    fn identical_inputs_produce_identical_layouts() {
        // The old one-shot layout was seeded for reproducibility; the sim must not
        // regress that, or tests and screenshots become non-deterministic.
        let run = || {
            let mut s = ForceSimulation::new(SimulationSpec::default());
            s.sync(&graph_of(&["a", "b", "c"], &[("e1", "a", "b"), ("e2", "b", "c")]), &HashMap::new());
            s.reheat_full();
            s.settle(500);
            ["a", "b", "c"].iter().map(|i| s.position(i).unwrap()).collect::<Vec<_>>()
        };
        assert_eq!(run(), run());
    }

    /// A flat 0.5/0.5 spring split gives a degree-k node a summed stiffness of
    /// k/2, which explicit Euler cannot integrate: measured, it was stable to
    /// degree 18, diverging by 20, and NaN at 26. The demo's seed graph has a
    /// node with 26 neighbours, so this is the shipping configuration, not a
    /// pathological one.
    #[test]
    fn a_high_degree_hub_does_not_explode() {
        let ids: Vec<String> = (0..=26).map(|i| format!("n{i:03}")).collect();
        let refs: Vec<&str> = ids.iter().map(|s| s.as_str()).collect();
        let edges: Vec<(String, String, String)> = (1..=26)
            .map(|i| (format!("e{i}"), "n000".to_string(), format!("n{i:03}")))
            .collect();
        let eref: Vec<(&str, &str, &str)> =
            edges.iter().map(|(i, a, b)| (i.as_str(), a.as_str(), b.as_str())).collect();

        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&refs, &eref), &HashMap::new());
        s.reheat_full();
        s.settle(3000);

        let hub = s.position("n000").unwrap();
        assert!(hub[0].is_finite() && hub[1].is_finite(), "hub went non-finite: {hub:?}");
        for id in &refs {
            let p = s.position(id).unwrap();
            assert!(p[0].is_finite() && p[1].is_finite(), "{id} went non-finite: {p:?}");
            let r = ((p[0] - hub[0]).powi(2) + (p[1] - hub[1]).powi(2)).sqrt();
            assert!(r < 2000.0, "{id} flew to {r} from the hub; the spring is diverging");
        }
    }

    /// Degree normalisation must not simply freeze every spring: a star's
    /// spokes should still be reeled in toward `link_distance`.
    #[test]
    fn a_hub_still_reels_its_spokes_to_roughly_link_distance() {
        let ids: Vec<String> = (0..=12).map(|i| format!("n{i:03}")).collect();
        let refs: Vec<&str> = ids.iter().map(|s| s.as_str()).collect();
        let edges: Vec<(String, String, String)> = (1..=12)
            .map(|i| (format!("e{i}"), "n000".to_string(), format!("n{i:03}")))
            .collect();
        let eref: Vec<(&str, &str, &str)> =
            edges.iter().map(|(i, a, b)| (i.as_str(), a.as_str(), b.as_str())).collect();

        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&refs, &eref), &HashMap::new());
        s.reheat_full();
        s.settle(3000);

        let hub = s.position("n000").unwrap();
        for i in 1..=12 {
            let p = s.position(&format!("n{i:03}")).unwrap();
            let r = ((p[0] - hub[0]).powi(2) + (p[1] - hub[1]).powi(2)).sqrt();
            assert!((30.0..600.0).contains(&r), "spoke n{i:03} settled at {r} from the hub");
        }
    }

    /// `Graph::edges()` iterates a HashMap, so two independently built copies
    /// of one graph can enumerate edges in different orders. Without a sort in
    /// `sync` the non-associative f32 spring sum diverges the layouts — this
    /// reproduced 26 times in 60 runs before `links.sort_unstable()`.
    #[test]
    fn link_order_is_independent_of_edge_hash_order() {
        let build = || {
            let mut s = ForceSimulation::new(SimulationSpec::default());
            s.sync(
                &graph_of(
                    &["a", "b", "c", "d", "e", "f"],
                    &[("e1", "a", "b"), ("e2", "b", "c"), ("e3", "c", "d"),
                      ("e4", "d", "e"), ("e5", "e", "f"), ("e6", "f", "a")],
                ),
                &HashMap::new(),
            );
            s.links.clone()
        };
        for _ in 0..50 {
            assert_eq!(build(), build(), "links order must not depend on hash iteration");
        }
    }

    /// Build a wide ring of seeded nodes, so the centroid is well defined and
    /// far from any one member.
    fn big_seeded_ring(n: usize) -> (Vec<String>, HashMap<String, [f32; 2]>) {
        let ids: Vec<String> = (0..n).map(|k| format!("n{k:04}")).collect();
        let seed = ids
            .iter()
            .enumerate()
            .map(|(k, id)| {
                let a = k as f32 / n as f32 * std::f32::consts::TAU;
                (id.clone(), [400.0 * a.cos(), 400.0 * a.sin()])
            })
            .collect();
        (ids, seed)
    }

    /// The spiral radius used to be `30·√(global rank)`, so a single unseeded
    /// arrival into a 500-node graph landed 671 units out — 5.6 link distances,
    /// nowhere near the "near the centroid" the doc contract promises.
    #[test]
    fn a_late_unseeded_arrival_lands_near_the_centroid() {
        let (ids, seed) = big_seeded_ring(500);
        let refs: Vec<&str> = ids.iter().map(|s| s.as_str()).collect();
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&refs, &[]), &seed);
        let c = s.centroid();

        let mut with_new = refs.clone();
        with_new.push("zzz-arrives-last");
        s.sync(&graph_of(&with_new, &[]), &HashMap::new());

        let p = s.position("zzz-arrives-last").unwrap();
        let r = ((p[0] - c[0]).powi(2) + (p[1] - c[1]).powi(2)).sqrt();
        let link = s.spec().link_distance;
        assert!(
            r < 2.0 * link,
            "arrival landed {r} from the centroid ({:.1} link distances), not near it",
            r / link
        );
    }

    /// Radius came from the id's position in the *sorted* id list, so the same
    /// node entering the same graph spawned at r≈21 or r≈671 purely according
    /// to how its id happened to sort. Nothing physical justifies that.
    #[test]
    fn an_arrival_position_does_not_depend_on_where_its_id_sorts() {
        let radius_for = |new_id: &str| {
            let (ids, seed) = big_seeded_ring(500);
            let refs: Vec<&str> = ids.iter().map(|s| s.as_str()).collect();
            let mut s = ForceSimulation::new(SimulationSpec::default());
            s.sync(&graph_of(&refs, &[]), &seed);
            let c = s.centroid();
            let mut with_new = refs.clone();
            with_new.push(new_id);
            s.sync(&graph_of(&with_new, &[]), &HashMap::new());
            let p = s.position(new_id).unwrap();
            ((p[0] - c[0]).powi(2) + (p[1] - c[1]).powi(2)).sqrt()
        };
        // "aaa-new" sorts before every "nNNNN"; "zzz-new" sorts after them all.
        let (first, last) = (radius_for("aaa-new"), radius_for("zzz-new"));
        assert!(
            (first - last).abs() < 1e-3,
            "spawn radius depends on id ordering: {first} when sorting first, {last} when last"
        );
    }

    /// The spawn counter lives on the struct, not on the `sync` call. A
    /// per-sync counter would hand every solitary arrival `seq = 0` and stack
    /// them all on one point — the precise failure the spiral exists to avoid.
    #[test]
    fn arrivals_in_separate_syncs_do_not_stack() {
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&["base"], &[]), &seeded(&[("base", [0.0, 0.0])]));

        let mut live = vec!["base".to_string()];
        let mut arrivals = Vec::new();
        for k in 0..40 {
            let id = format!("late{k:03}");
            live.push(id.clone());
            let refs: Vec<&str> = live.iter().map(|x| x.as_str()).collect();
            s.sync(&graph_of(&refs, &[]), &HashMap::new()); // one arrival per sync
            arrivals.push(s.position(&id).unwrap());
        }

        for i in 0..arrivals.len() {
            for j in (i + 1)..arrivals.len() {
                let (p, q) = (arrivals[i], arrivals[j]);
                let d = ((p[0] - q[0]).powi(2) + (p[1] - q[1]).powi(2)).sqrt();
                assert!(d > 1.0, "arrivals {i} and {j} spawned {d} apart, effectively stacked");
            }
        }
        // And they must all still be near the centroid, not spiralling outward.
        let link = s.spec().link_distance;
        for (k, p) in arrivals.iter().enumerate() {
            let r = (p[0] * p[0] + p[1] * p[1]).sqrt();
            assert!(r < 2.0 * link, "arrival {k} spawned {r} out; the radius is not bounded");
        }
    }

    /// A cold start is the one case where spiral rank really is arrival count,
    /// and its `30·√N` spread is correct today. Bounding the *arrival* radius
    /// must not quietly shrink it.
    #[test]
    fn a_cold_start_keeps_its_full_spiral_spread() {
        let ids: Vec<String> = (0..500).map(|k| format!("n{k:04}")).collect();
        let refs: Vec<&str> = ids.iter().map(|s| s.as_str()).collect();
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&refs, &[]), &HashMap::new());

        let c = s.centroid();
        let max = s
            .positions()
            .map(|(_, p)| ((p[0] - c[0]).powi(2) + (p[1] - c[1]).powi(2)).sqrt())
            .fold(0.0f32, f32::max);
        let expected = 30.0 * 499.5f32.sqrt(); // the outermost seat, ~670
        assert!(
            (max - expected).abs() < 10.0,
            "cold-start spread is {max}, expected ~{expected}; the sunflower was truncated"
        );
    }

    /// Above theta ≈ 1.25 — and ≈ 0.75 for a tight cluster with a distant
    /// outlier, which is exactly the shape an expanded aggregate makes — a cell
    /// that contains body `i` can pass the approximation test, and `i` starts
    /// repelling itself through its own cell's mass. The spec boundary is where
    /// that gets excluded, so no per-tick branch has to.
    #[test]
    fn theta_is_clamped_when_the_spec_is_adopted() {
        let s = ForceSimulation::new(SimulationSpec { theta: 4.0, ..SimulationSpec::default() });
        assert_eq!(s.spec().theta, 1.0, "new() must clamp an out-of-range theta");

        let mut t = ForceSimulation::new(SimulationSpec::default());
        t.set_spec(SimulationSpec { theta: -3.0, ..SimulationSpec::default() });
        assert_eq!(t.spec().theta, 0.0, "set_spec must clamp a negative theta to exact O(n^2)");

        let u = ForceSimulation::new(SimulationSpec { theta: 0.9, ..SimulationSpec::default() });
        assert_eq!(u.spec().theta, 0.9, "an in-range theta must be left exactly alone");
    }

    #[test]
    fn the_alpha_schedule_is_clamped_when_the_spec_is_adopted() {
        // Both entry points, since `set_spec` is the one a host reaches in a
        // single call and is therefore the one that has to be safe.
        let s = ForceSimulation::new(SimulationSpec {
            alpha_min: 0.0,
            alpha_decay: 0.0,
            velocity_decay: -1.0,
            reheat_alpha: 7.0,
            ..SimulationSpec::default()
        });
        assert_eq!(s.spec().alpha_min, ALPHA_MIN_FLOOR);
        assert_eq!(s.spec().alpha_decay, ALPHA_DECAY_MIN);
        assert_eq!(s.spec().velocity_decay, VELOCITY_DECAY_MIN);
        assert_eq!(s.spec().reheat_alpha, 1.0);

        let mut t = ForceSimulation::new(SimulationSpec::default());
        t.set_spec(SimulationSpec {
            alpha_min: 5.0,
            alpha_decay: 9.0,
            velocity_decay: 3.0,
            reheat_alpha: -2.0,
            ..SimulationSpec::default()
        });
        assert_eq!(t.spec().alpha_min, ALPHA_MIN_CEIL);
        assert_eq!(t.spec().alpha_decay, ALPHA_DECAY_MAX);
        assert_eq!(t.spec().velocity_decay, VELOCITY_DECAY_MAX);
        assert_eq!(t.spec().reheat_alpha, 0.0);

        // The defaults sit inside every bound and must survive untouched, or
        // the clamps would be silently retuning the shipped schedule.
        let d = ForceSimulation::new(SimulationSpec::default());
        let want = SimulationSpec::default();
        assert_eq!(d.spec().alpha_min, want.alpha_min);
        assert_eq!(d.spec().alpha_decay, want.alpha_decay);
        assert_eq!(d.spec().velocity_decay, want.velocity_decay);
        assert_eq!(d.spec().reheat_alpha, want.reheat_alpha);
    }

    #[test]
    fn a_zero_alpha_min_cannot_wedge_the_sim_active_forever() {
        // `is_active` is `alpha >= alpha_min` and `alpha` decays multiplicatively
        // toward zero without reaching it, so an unclamped `alpha_min: 0` is
        // permanently true: the host's render loop would tick Barnes-Hut every
        // frame at rest, and the settle edge that callers rebuild their scene on
        // would never fire.
        let mut s = ForceSimulation::new(SimulationSpec {
            alpha_min: 0.0,
            ..SimulationSpec::default()
        });
        s.sync(&graph_of(&["a", "b"], &[("e", "a", "b")]), &HashMap::new());
        s.reheat_full();
        s.settle(5000);
        assert!(!s.is_active(), "a zero alpha_min must not leave the sim active forever");
        assert!(!s.tick(), "and the next tick must still report settled");
    }

    #[test]
    fn every_legal_schedule_settles_within_the_reduced_motion_budget() {
        // The guarantee `ALPHA_DECAY_MIN` is derived from. A reduced-motion
        // caller settles synchronously with `settle(600)` and then snaps; if 600
        // ticks is not enough the sim is still active afterwards and animates
        // every frame — "no visible motion" inverted into permanent motion, plus
        // 600 wasted Barnes-Hut ticks per reheat. Each spec below pushes a
        // different schedule value to the end that used to break it.
        let hostile = [
            ("zero alpha_min", SimulationSpec { alpha_min: 0.0, ..Default::default() }),
            ("zero alpha_decay", SimulationSpec { alpha_decay: 0.0, ..Default::default() }),
            ("negative alpha_decay", SimulationSpec { alpha_decay: -1.0, ..Default::default() }),
            ("zero velocity_decay", SimulationSpec { velocity_decay: 0.0, ..Default::default() }),
            ("oversized reheat_alpha", SimulationSpec { reheat_alpha: 9.0, ..Default::default() }),
            (
                "all of them at once",
                SimulationSpec {
                    alpha_min: 0.0,
                    alpha_decay: 0.0,
                    velocity_decay: -3.0,
                    reheat_alpha: 9.0,
                    ..Default::default()
                },
            ),
        ];
        for (what, spec) in hostile {
            let mut s = ForceSimulation::new(spec);
            s.sync(
                &graph_of(&["a", "b", "c"], &[("e1", "a", "b"), ("e2", "b", "c")]),
                &HashMap::new(),
            );
            // `reheat_full`, not `reheat`: alpha 1.0 is the worst case the
            // budget has to cover.
            s.reheat_full();
            s.settle(600);
            assert!(!s.is_active(), "{what}: still active after settle(600)");
            assert!(!s.tick(), "{what}: the next tick must still report settled");
        }
    }

    #[test]
    fn a_negative_velocity_decay_cannot_diverge_the_layout() {
        // `friction = 1 - velocity_decay`, so a negative decay makes friction
        // exceed 1 and this explicit-Euler step amplifies velocity every tick
        // instead of damping it — unclamped, the positions reach NaN.
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.set_spec(SimulationSpec { velocity_decay: -5.0, ..SimulationSpec::default() });
        s.sync(
            &graph_of(&["a", "b", "c"], &[("e1", "a", "b"), ("e2", "b", "c")]),
            &seeded(&[("a", [0.0, 0.0]), ("b", [10.0, 0.0]), ("c", [0.0, 10.0])]),
        );
        s.reheat_full();
        s.settle(2000);
        for id in ["a", "b", "c"] {
            let p = s.position(id).expect("node present");
            assert!(p[0].is_finite() && p[1].is_finite(), "{id} diverged to {p:?}");
        }
    }

    #[test]
    fn the_alpha_decay_floor_leaves_headroom_under_the_settle_budget() {
        // Pins the derivation in `ALPHA_DECAY_MIN`'s doc, so that lowering the
        // floor (or raising `ALPHA_MIN_FLOOR`) without redoing the arithmetic
        // fails here rather than silently in reduced motion.
        let ticks = (ALPHA_MIN_FLOOR.ln() / (1.0 - ALPHA_DECAY_MIN).ln()).ceil();
        assert!(ticks <= 600.0, "worst-case settle is {ticks} ticks, over the 600 budget");
    }

    #[test]
    fn coincident_start_positions_do_not_produce_nan() {
        let mut s = ForceSimulation::new(SimulationSpec::default());
        s.sync(&graph_of(&["a", "b", "c"], &[]),
               &seeded(&[("a", [0.0, 0.0]), ("b", [0.0, 0.0]), ("c", [0.0, 0.0])]));
        s.reheat_full();
        s.settle(200);
        for id in ["a", "b", "c"] {
            let p = s.position(id).unwrap();
            assert!(p[0].is_finite() && p[1].is_finite(), "{id} = {p:?}");
        }
    }
}