manatee 0.4.0

#![allow(dead_code)]
#![allow(clippy::collapsible_if)]
#![allow(clippy::manual_div_ceil)]
#![allow(clippy::needless_option_as_deref)]
#![allow(clippy::needless_range_loop)]
#![allow(clippy::nonminimal_bool)]

use crate::algo::FcoseOptions;
use crate::error::Result;
use crate::graph::{Anchor, BoundsExtras, Graph, LayoutResult, Point};
use indexmap::{IndexMap, IndexSet};
use nalgebra as na;
use rustc_hash::FxHashMap;

mod spectral;

#[derive(Debug, Default, Clone)]
struct FcoseLayoutTimings {
    total: std::time::Duration,
    from_graph: std::time::Duration,
    constraints: std::time::Duration,
    spring: FcoseSpringTimings,
    translate: std::time::Duration,
    output: std::time::Duration,
}

#[derive(Debug, Default, Clone)]
struct FcoseSpringTimings {
    total: std::time::Duration,
    opts_prep: std::time::Duration,
    spectral: std::time::Duration,
    root_compound: std::time::Duration,
    collapse_start_positions: std::time::Duration,
    pre_constraints: std::time::Duration,
    constraint_rt: std::time::Duration,
    iterations: std::time::Duration,
}

#[derive(Debug, Default, Clone, Copy)]
struct SpringStats {
    iterations: usize,
    spectral_applied: bool,
}

pub fn layout(graph: &Graph, opts: &FcoseOptions) -> Result<LayoutResult> {
    graph.validate()?;

    let timing_enabled = std::env::var("MANATEE_FCOSE_TIMING").ok().as_deref() == Some("1");
    let mut timings = FcoseLayoutTimings::default();
    let total_start = timing_enabled.then(std::time::Instant::now);

    let from_graph_start = timing_enabled.then(std::time::Instant::now);
    let mut sim = SimGraph::from_graph(graph);
    if let Some(s) = from_graph_start {
        timings.from_graph = s.elapsed();
    }

    let constraints_start = timing_enabled.then(std::time::Instant::now);
    let constraints = Constraints::from_opts(&sim, opts);
    if let Some(s) = constraints_start {
        timings.constraints = s.elapsed();
    }

    let mut rng = XorShift64Star::new(opts.random_seed);
    // Mermaid upstream SVG baselines (ADR-0055) seed `Math.random()` at document start; a small
    // amount of randomness can be consumed before the first spectral sampler draw. Model that as
    // a single deterministic "seed offset" per layout invocation (not per rerun).
    let _ = rng.next_f64_unit();
    let run_count = if opts.rerun { 2 } else { 1 };
    let mut spring_stats = SpringStats::default();

    let debug_rng_calls = std::env::var("MANATEE_FCOSE_DEBUG_RNG_CALLS")
        .ok()
        .as_deref()
        == Some("1");

    for run_idx in 0..run_count {
        let rng_calls_before = rng.calls();
        // Mirror upstream component center bookkeeping (`eles.boundingBox()` before layout) by
        // ensuring compound rects wrap their current children before we compute `orig_center`.
        let compound_padding = opts.compound_padding.unwrap_or(0.0).max(0.0);
        if compound_padding > 0.0 {
            for n in &mut sim.nodes {
                if n.is_compound {
                    n.padding = compound_padding;
                }
            }
        }
        let _ = sim.update_bounds();

        // Mimic fcose's `aux.relocateComponent(...)`: keep the final component center aligned to
        // the original component center to avoid arbitrary global translations affecting viewBox
        // parity.
        //
        // Upstream uses Cytoscape `eles.boundingBox()` to capture the pre-layout component center,
        // and then relocates the final node rects to that center. Importantly, compounds are part
        // of that bbox, and compound sizes include padding (and may include label sizing depending
        // on style). Using leaves-only centers creates deterministic root drift for group-heavy
        // diagrams (e.g. architecture groups-within-groups).
        //
        // Mermaid Architecture runs Cytoscape FCoSE twice (`layout.run()` in `layoutstop`), so we
        // repeat this per run while keeping the RNG stream continuous.
        // Upstream Cytoscape FCoSE keeps the final component aligned to the *pre-layout*
        // `options.eles.boundingBox()` center (nodes + edges + labels).
        //
        // Important: In proof/default quality, `aux.relocateComponent(...)` computes the
        // "current" bbox from layout-base node rects (`node.getRect()`), which excludes labels
        // when `nodeDimensionsIncludeLabels: false`.
        let orig_center = opts
            .relocate_center
            .or_else(|| sim.bounding_box_center_eles(run_idx))
            .unwrap_or((0.0, 0.0));
        let debug_relocate = std::env::var("MANATEE_FCOSE_DEBUG_RELOCATE")
            .ok()
            .as_deref()
            == Some("1");

        let spring_start = timing_enabled.then(std::time::Instant::now);
        spring_stats = sim.run_spring_embedder(
            &constraints,
            opts,
            &mut rng,
            if timing_enabled {
                Some(&mut timings.spring)
            } else {
                None
            },
        );
        if debug_rng_calls {
            eprintln!(
                "[manatee-fcose-rng] run={} calls_before={} calls_after={} (+{})",
                run_idx,
                rng_calls_before,
                rng.calls(),
                rng.calls().saturating_sub(rng_calls_before)
            );
        }
        if let Some(s) = spring_start {
            timings.spring.total = s.elapsed();
        }

        // Ensure compound node rectangles reflect the final child placements before we compute the
        // "current" component bounding box for relocation (`aux.relocateComponent(...)` parity).
        let _ = sim.update_bounds();

        let new_center = sim.bounding_box_center_rects().unwrap_or((0.0, 0.0));
        let translate_start = timing_enabled.then(std::time::Instant::now);
        let disable_relocate = std::env::var("MANATEE_FCOSE_DISABLE_RELOCATE")
            .ok()
            .as_deref()
            == Some("1");
        if !disable_relocate {
            let dx = orig_center.0 - new_center.0;
            let dy = orig_center.1 - new_center.1;
            if debug_relocate {
                eprintln!(
                    "[manatee-fcose-relocate] run={} orig=({:.6},{:.6}) new=({:.6},{:.6}) d=({:.6},{:.6})",
                    run_idx, orig_center.0, orig_center.1, new_center.0, new_center.1, dx, dy
                );
            }
            sim.translate(dx, dy);
        }
        if let Some(s) = translate_start {
            timings.translate = s.elapsed();
        }

        if run_idx + 1 < run_count {
            let _ = sim.update_bounds();
        }
    }

    let output_start = timing_enabled.then(std::time::Instant::now);
    let node_count = sim.nodes.len();
    let edge_count = sim.edges.len();
    let compound_count = sim.compound_parent.len();

    let mut positions: std::collections::BTreeMap<String, Point> =
        std::collections::BTreeMap::new();
    let nodes = std::mem::take(&mut sim.nodes);
    for n in nodes {
        let x = n.center_x();
        let y = n.center_y();
        positions.insert(n.id, Point { x, y });
    }
    if let Some(s) = output_start {
        timings.output = s.elapsed();
    }

    if let Some(s) = total_start {
        timings.total = s.elapsed();
        eprintln!(
            "[manatee-fcose-timing] total={:?} from_graph={:?} constraints={:?} spring_total={:?} spring_opts_prep={:?} spring_spectral={:?} spring_root_compound={:?} spring_collapse_start={:?} spring_pre_constraints={:?} spring_constraint_rt={:?} spring_iterations={:?} translate={:?} output={:?} nodes={} edges={} compounds={} iterations={} spectral_applied={}",
            timings.total,
            timings.from_graph,
            timings.constraints,
            timings.spring.total,
            timings.spring.opts_prep,
            timings.spring.spectral,
            timings.spring.root_compound,
            timings.spring.collapse_start_positions,
            timings.spring.pre_constraints,
            timings.spring.constraint_rt,
            timings.spring.iterations,
            timings.translate,
            timings.output,
            node_count,
            edge_count,
            compound_count,
            spring_stats.iterations,
            spring_stats.spectral_applied,
        );
    }

    Ok(LayoutResult { positions })
}

#[derive(Debug, Clone)]
struct SimNode {
    id: String,
    parent: Option<String>,
    owner_idx: usize,
    is_compound: bool,
    width: f64,
    height: f64,
    bounds_extras: BoundsExtras,
    // layout-base `LNode.estimatedSize` (stable across updateBounds mutations).
    estimated_size: f64,
    // Top-left anchored rectangle (layout-base `LNode.rect` style).
    left: f64,
    top: f64,

    spring_fx: f64,
    spring_fy: f64,
    repulsion_fx: f64,
    repulsion_fy: f64,
    gravitation_fx: f64,
    gravitation_fy: f64,

    // layout-base `LNode.noOfChildren` weight (leaf descendants count).
    no_of_children: f64,

    // Compound padding (Cytoscape style `padding`, mapped onto layout-base `paddingLeft/...`).
    // Used as a margin when computing child graph bounds.
    padding: f64,

    // layout-base FR-grid repulsion caches a per-node "surrounding" list, refreshed periodically.
    surrounding: Vec<usize>,
    grid_start_x: i32,
    grid_finish_x: i32,
    grid_start_y: i32,
    grid_finish_y: i32,
}

impl SimNode {
    fn center_x(&self) -> f64 {
        self.left + self.width / 2.0
    }

    fn center_y(&self) -> f64 {
        self.top + self.height / 2.0
    }

    fn move_by(&mut self, dx: f64, dy: f64) {
        self.left += dx;
        self.top += dy;
    }

    fn half_w(&self) -> f64 {
        self.width / 2.0
    }

    fn half_h(&self) -> f64 {
        self.height / 2.0
    }

    fn right(&self) -> f64 {
        self.left + self.width
    }

    fn bottom(&self) -> f64 {
        self.top + self.height
    }

    fn bound_left(&self) -> f64 {
        self.left - self.bounds_extras.left.max(0.0)
    }

    fn bound_right(&self) -> f64 {
        self.right() + self.bounds_extras.right.max(0.0)
    }

    fn bound_top(&self) -> f64 {
        self.top - self.bounds_extras.top.max(0.0)
    }

    fn bound_bottom(&self) -> f64 {
        self.bottom() + self.bounds_extras.bottom.max(0.0)
    }
}

fn imath_sign(value: f64) -> f64 {
    // layout-base `IMath.sign`: returns 1, -1, or 0 (and yields 0 for NaN).
    if value > 0.0 {
        1.0
    } else if value < 0.0 {
        -1.0
    } else {
        0.0
    }
}

fn lca_owner_idx(nodes: &[SimNode], root_owner_idx: usize, a: usize, b: usize) -> usize {
    let mut seen: Vec<bool> = vec![false; nodes.len() + 1];

    let mut cur = nodes.get(a).map(|n| n.owner_idx).unwrap_or(root_owner_idx);
    loop {
        if cur >= seen.len() {
            break;
        }
        seen[cur] = true;
        if cur == root_owner_idx {
            break;
        }
        cur = nodes
            .get(cur)
            .map(|n| n.owner_idx)
            .unwrap_or(root_owner_idx);
    }

    let mut cur = nodes.get(b).map(|n| n.owner_idx).unwrap_or(root_owner_idx);
    loop {
        if cur < seen.len() && seen[cur] {
            return cur;
        }
        if cur == root_owner_idx {
            break;
        }
        cur = nodes
            .get(cur)
            .map(|n| n.owner_idx)
            .unwrap_or(root_owner_idx);
    }

    root_owner_idx
}

fn node_in_lca_idx(
    nodes: &[SimNode],
    root_owner_idx: usize,
    node_idx: usize,
    lca_owner: usize,
) -> usize {
    let Some(node) = nodes.get(node_idx) else {
        return node_idx;
    };
    if node.owner_idx == lca_owner {
        return node_idx;
    }

    let mut owner = node.owner_idx;
    while owner != root_owner_idx {
        let Some(parent_owner) = nodes.get(owner).map(|n| n.owner_idx) else {
            break;
        };
        if parent_owner == lca_owner {
            return owner;
        }
        owner = parent_owner;
    }

    node_idx
}

#[derive(Debug, Clone, Copy)]
struct SimEdge {
    a: usize,
    b: usize,
    // Cache LCA-lifted endpoints for spring forces (layout-base `LEdge.getSourceInLca/getTargetInLca`).
    //
    // For inter-graph edges, CoSE applies spring forces between the *immediate children* of the
    // edge's lowest common ancestor owner graph (often compound nodes), rather than pulling the
    // original leaf endpoints across compound boundaries.
    a_in_lca: usize,
    b_in_lca: usize,
    #[allow(dead_code)]
    a_anchor: Option<Anchor>,
    #[allow(dead_code)]
    b_anchor: Option<Anchor>,
    base_ideal_length: f64,
    ideal_length: f64,
    elasticity: f64,
    label_width: Option<f64>,
    label_height: Option<f64>,
}

#[derive(Debug, Clone)]
struct Constraints {
    align_horizontal: Vec<Vec<usize>>,
    align_vertical: Vec<Vec<usize>>,
    relative: Vec<RelConstraint>,
}

#[derive(Debug, Clone, Copy)]
struct RelConstraint {
    left: Option<usize>,
    right: Option<usize>,
    top: Option<usize>,
    bottom: Option<usize>,
    gap: f64,
}

impl Constraints {
    fn from_opts(sim: &SimGraph, opts: &FcoseOptions) -> Self {
        let mut align_horizontal: Vec<Vec<usize>> = Vec::new();
        let mut align_vertical: Vec<Vec<usize>> = Vec::new();

        if let Some(a) = opts.alignment_constraint.as_ref() {
            align_horizontal = map_align_lists(sim, &a.horizontal);
            align_vertical = map_align_lists(sim, &a.vertical);
        }

        let mut relative: Vec<RelConstraint> = Vec::new();
        for r in &opts.relative_placement_constraint {
            relative.push(RelConstraint {
                left: r
                    .left
                    .as_deref()
                    .and_then(|id| sim.id_to_idx.get(id).copied()),
                right: r
                    .right
                    .as_deref()
                    .and_then(|id| sim.id_to_idx.get(id).copied()),
                top: r
                    .top
                    .as_deref()
                    .and_then(|id| sim.id_to_idx.get(id).copied()),
                bottom: r
                    .bottom
                    .as_deref()
                    .and_then(|id| sim.id_to_idx.get(id).copied()),
                gap: r.gap.max(0.0),
            });
        }

        Self {
            align_horizontal,
            align_vertical,
            relative,
        }
    }
}

#[derive(Debug, Clone, Copy)]
enum Axis {
    Horizontal,
    Vertical,
}

#[derive(Debug, Clone)]
struct ConstraintRuntime {
    horizontal: AxisConstraintRuntime,
    vertical: AxisConstraintRuntime,
}

#[derive(Debug, Clone)]
struct AxisConstraintRuntime {
    node_count: usize,
    dummy_to_nodes: Vec<Vec<usize>>,
    #[allow(dead_code)]
    node_to_dummy: Vec<Option<usize>>,
    fixed_nodes: IndexSet<usize>,
    nodes_in_relative: Vec<usize>,
    rel_map: Vec<Vec<AxisRelAdj>>,
    temp_pos: Vec<f64>,
}

#[derive(Debug, Clone, Copy)]
enum AxisRelAdj {
    Right { node: usize, gap: f64 },
    Left { node: usize, gap: f64 },
    Bottom { node: usize, gap: f64 },
    Top { node: usize, gap: f64 },
}

impl ConstraintRuntime {
    fn new(nodes: &[SimNode], c: &Constraints) -> Option<Self> {
        if c.relative.is_empty() {
            return None;
        }
        Some(Self {
            horizontal: AxisConstraintRuntime::new_axis(
                nodes,
                &c.align_vertical,
                &c.relative,
                Axis::Horizontal,
            ),
            vertical: AxisConstraintRuntime::new_axis(
                nodes,
                &c.align_horizontal,
                &c.relative,
                Axis::Vertical,
            ),
        })
    }

    fn update_displacements(
        &mut self,
        _nodes: &[SimNode],
        c: &Constraints,
        disps: &mut [(f64, f64)],
        total_iterations: usize,
        _max_d: f64,
        rng: &mut XorShift64Star,
    ) {
        // Fixed nodes (not currently exposed by our public API).
        for &idx in &self.horizontal.fixed_nodes {
            if idx < disps.len() {
                disps[idx].0 = 0.0;
            }
        }
        for &idx in &self.vertical.fixed_nodes {
            if idx < disps.len() {
                disps[idx].1 = 0.0;
            }
        }

        // Alignments (match `cose-base` updateDisplacements): average displacements per group.
        for group in &c.align_vertical {
            if group.len() <= 1 {
                continue;
            }
            let mut sum = 0.0;
            for &idx in group {
                sum += disps[idx].0;
            }
            let avg = sum / (group.len() as f64);
            for &idx in group {
                disps[idx].0 = avg;
            }
        }
        for group in &c.align_horizontal {
            if group.len() <= 1 {
                continue;
            }
            let mut sum = 0.0;
            for &idx in group {
                sum += disps[idx].1;
            }
            let avg = sum / (group.len() as f64);
            for &idx in group {
                disps[idx].1 = avg;
            }
        }

        // Relative placements (match `cose-base` relax-movement mode).
        // Upstream keeps `nodeToTempPositionMap*` as a persistent accumulator across iterations:
        // it starts from node centers and is advanced by the chosen displacements each tick.
        // Do not re-seed from node centers here, or the relaxation order differs.
        if total_iterations.is_multiple_of(10) {
            self.horizontal.shuffle_tail_third(rng);
            self.vertical.shuffle_tail_third(rng);
        }

        self.horizontal
            .apply_relative_relaxation(disps, Axis::Horizontal);
        self.vertical
            .apply_relative_relaxation(disps, Axis::Vertical);
    }
}

impl AxisConstraintRuntime {
    fn new_axis(
        nodes: &[SimNode],
        axis_alignment_groups: &[Vec<usize>],
        rel: &[RelConstraint],
        axis: Axis,
    ) -> Self {
        let n = nodes.len();
        let d = axis_alignment_groups.len();

        let mut node_to_dummy: Vec<Option<usize>> = vec![None; n];
        let mut dummy_to_nodes: Vec<Vec<usize>> = Vec::with_capacity(d);
        for (i, group) in axis_alignment_groups.iter().enumerate() {
            let dummy_key = n + i;
            dummy_to_nodes.push(group.clone());
            for &idx in group {
                if idx < n {
                    node_to_dummy[idx] = Some(dummy_key);
                }
            }
        }

        let key_count = n + d;
        let mut rel_map: Vec<Vec<AxisRelAdj>> = vec![Vec::new(); key_count];
        let mut nodes_in_relative_set: IndexSet<usize> = IndexSet::new();

        for r in rel {
            match axis {
                Axis::Horizontal => {
                    let (Some(left), Some(right)) = (r.left, r.right) else {
                        continue;
                    };
                    let lk = node_to_dummy.get(left).copied().flatten().unwrap_or(left);
                    let rk = node_to_dummy.get(right).copied().flatten().unwrap_or(right);
                    nodes_in_relative_set.insert(lk);
                    nodes_in_relative_set.insert(rk);
                    rel_map[lk].push(AxisRelAdj::Right {
                        node: rk,
                        gap: r.gap,
                    });
                    rel_map[rk].push(AxisRelAdj::Left {
                        node: lk,
                        gap: r.gap,
                    });
                }
                Axis::Vertical => {
                    let (Some(top), Some(bottom)) = (r.top, r.bottom) else {
                        continue;
                    };
                    let tk = node_to_dummy.get(top).copied().flatten().unwrap_or(top);
                    let bk = node_to_dummy
                        .get(bottom)
                        .copied()
                        .flatten()
                        .unwrap_or(bottom);
                    nodes_in_relative_set.insert(tk);
                    nodes_in_relative_set.insert(bk);
                    rel_map[tk].push(AxisRelAdj::Bottom {
                        node: bk,
                        gap: r.gap,
                    });
                    rel_map[bk].push(AxisRelAdj::Top {
                        node: tk,
                        gap: r.gap,
                    });
                }
            }
        }

        let mut rt = Self {
            node_count: n,
            dummy_to_nodes,
            node_to_dummy,
            fixed_nodes: IndexSet::new(),
            nodes_in_relative: nodes_in_relative_set.into_iter().collect(),
            rel_map,
            temp_pos: vec![0.0; key_count],
        };
        rt.refresh_temp_positions(nodes, axis);
        rt
    }

    fn refresh_temp_positions(&mut self, nodes: &[SimNode], axis: Axis) {
        let n = self.node_count;
        for key in 0..self.temp_pos.len() {
            let v = if key < n {
                match axis {
                    Axis::Horizontal => nodes[key].center_x(),
                    Axis::Vertical => nodes[key].center_y(),
                }
            } else {
                let dummy_idx = key - n;
                let first = self.dummy_to_nodes[dummy_idx]
                    .first()
                    .copied()
                    .unwrap_or(0)
                    .min(n.saturating_sub(1));
                match axis {
                    Axis::Horizontal => nodes[first].center_x(),
                    Axis::Vertical => nodes[first].center_y(),
                }
            };
            self.temp_pos[key] = v;
        }
    }

    fn shuffle_tail_third(&mut self, rng: &mut XorShift64Star) {
        let len = self.nodes_in_relative.len();
        if len <= 1 {
            return;
        }
        // Upstream (`cose-base`) uses:
        //
        // `for (i = len - 1; i >= (2 * len / 3); i--)`
        //
        // where `(2 * len / 3)` is a JS Number (not integer division). Therefore the effective
        // lower bound is `ceil(2 * len / 3)`.
        let start = (2 * len + 2) / 3;
        for i in (start..len).rev() {
            let j = rng.next_usize(i + 1);
            self.nodes_in_relative.swap(i, j);
        }
    }

    fn apply_relative_relaxation(&mut self, disps: &mut [(f64, f64)], axis: Axis) {
        let n = self.node_count;

        for &key in &self.nodes_in_relative {
            if self.fixed_nodes.contains(&key) {
                continue;
            }

            let mut displacement = if key < n {
                match axis {
                    Axis::Horizontal => disps[key].0,
                    Axis::Vertical => disps[key].1,
                }
            } else {
                let dummy_idx = key - n;
                let first = self.dummy_to_nodes[dummy_idx]
                    .first()
                    .copied()
                    .unwrap_or(0)
                    .min(n.saturating_sub(1));
                match axis {
                    Axis::Horizontal => disps[first].0,
                    Axis::Vertical => disps[first].1,
                }
            };

            for adj in &self.rel_map[key] {
                match (*adj, axis) {
                    (AxisRelAdj::Right { node, gap }, Axis::Horizontal) => {
                        let diff = (self.temp_pos[node] - self.temp_pos[key]) - displacement;
                        if diff < gap {
                            displacement -= gap - diff;
                        }
                    }
                    (AxisRelAdj::Left { node, gap }, Axis::Horizontal) => {
                        let diff = (self.temp_pos[key] - self.temp_pos[node]) + displacement;
                        if diff < gap {
                            displacement += gap - diff;
                        }
                    }
                    (AxisRelAdj::Bottom { node, gap }, Axis::Vertical) => {
                        let diff = (self.temp_pos[node] - self.temp_pos[key]) - displacement;
                        if diff < gap {
                            displacement -= gap - diff;
                        }
                    }
                    (AxisRelAdj::Top { node, gap }, Axis::Vertical) => {
                        let diff = (self.temp_pos[key] - self.temp_pos[node]) + displacement;
                        if diff < gap {
                            displacement += gap - diff;
                        }
                    }
                    _ => {}
                }
            }

            self.temp_pos[key] += displacement;

            if key < n {
                match axis {
                    Axis::Horizontal => disps[key].0 = displacement,
                    Axis::Vertical => disps[key].1 = displacement,
                }
            } else {
                let dummy_idx = key - n;
                for &idx in &self.dummy_to_nodes[dummy_idx] {
                    if idx >= disps.len() {
                        continue;
                    }
                    match axis {
                        Axis::Horizontal => disps[idx].0 = displacement,
                        Axis::Vertical => disps[idx].1 = displacement,
                    }
                }
            }
        }
    }
}

fn map_align_lists(sim: &SimGraph, groups: &[Vec<String>]) -> Vec<Vec<usize>> {
    // Preserve Mermaid/Cytoscape ordering (and duplicates) for alignment arrays.
    //
    // Upstream `ConstraintHandler` uses the *first* node id in each alignment group as the seed
    // for dummy-node positions in the relative-placement enforcement phase. Sorting/deduping
    // here changes that seed and can shift the entire layout in parity-root mode.
    let mut out: Vec<Vec<usize>> = Vec::new();
    for g in groups {
        let idxs: Vec<usize> = g
            .iter()
            .filter_map(|id| sim.id_to_idx.get(id.as_str()).copied())
            .collect();
        if idxs.len() > 1 {
            out.push(idxs);
        }
    }
    out
}

#[derive(Debug)]
struct SimGraph {
    nodes: Vec<SimNode>,
    edges: Vec<SimEdge>,
    id_to_idx: FxHashMap<String, usize>,
    compound_parent: FxHashMap<String, Option<String>>,
    compound_ids_in_order: Vec<String>,
    leaf_count: usize,
    // Owner graph identity for repulsion/gravity: each node belongs to the child graph of its
    // parent compound, or the root graph.
    root_owner_idx: usize,
    // Immediate children list for each owner graph (owner idx is a compound node idx, or
    // `root_owner_idx` for the root graph).
    children_by_owner: Vec<Vec<usize>>,
    // Compound node indices in descending inclusion depth (deepest first), for updateBounds.
    compounds_deep_first: Vec<usize>,
    // Descendant leaf indices for each node (empty for leaves).
    descendant_leaves: Vec<Vec<usize>>,
    // Estimated size for each owner graph (static; computed from node sizes).
    owner_estimated_size: Vec<f64>,
    // layout-base `LNode.inclusionTreeDepth` (root-level nodes depth=1).
    inclusion_depth: Vec<usize>,
}

#[derive(Debug, Clone)]
struct OwnerBounds {
    left: Vec<f64>,
    right: Vec<f64>,
    top: Vec<f64>,
    bottom: Vec<f64>,
}

impl SimGraph {
    const DEFAULT_EDGE_LENGTH: f64 = 50.0;
    const DEFAULT_SPRING_STRENGTH: f64 = 0.45;
    const DEFAULT_REPULSION_STRENGTH: f64 = 4500.0;
    // cytoscape-fcose default (overrides layout-base default 0.4 via `options.gravity`).
    const DEFAULT_GRAVITY_STRENGTH: f64 = 0.25;
    const DEFAULT_COMPOUND_GRAVITY_STRENGTH: f64 = 1.0; // layout-base `FDLayoutConstants.DEFAULT_COMPOUND_GRAVITY_STRENGTH`
    const DEFAULT_GRAVITY_RANGE_FACTOR: f64 = 3.8; // layout-base `FDLayoutConstants.DEFAULT_GRAVITY_RANGE_FACTOR`
    const DEFAULT_COMPOUND_GRAVITY_RANGE_FACTOR: f64 = 1.5; // layout-base `FDLayoutConstants.DEFAULT_COMPOUND_GRAVITY_RANGE_FACTOR`
    const DEFAULT_GRAPH_MARGIN: f64 = 15.0; // layout-base `LayoutConstants.DEFAULT_GRAPH_MARGIN`
    const EMPTY_COMPOUND_NODE_SIZE: f64 = 40.0; // layout-base `LayoutConstants.EMPTY_COMPOUND_NODE_SIZE`
    const SIMPLE_NODE_SIZE: f64 = 40.0; // layout-base `LayoutConstants.SIMPLE_NODE_SIZE`
    const PER_LEVEL_IDEAL_EDGE_LENGTH_FACTOR: f64 = 0.1; // layout-base `FDLayoutConstants.PER_LEVEL_IDEAL_EDGE_LENGTH_FACTOR`
    const DEFAULT_COOLING_FACTOR_INCREMENTAL: f64 = 0.3; // layout-base `FDLayoutConstants.DEFAULT_COOLING_FACTOR_INCREMENTAL`
    const FINAL_TEMPERATURE: f64 = 0.04; // cose-base `CoSELayout.initSpringEmbedder()`
    const GRID_CALCULATION_CHECK_PERIOD: usize = 10; // layout-base `FDLayoutConstants.GRID_CALCULATION_CHECK_PERIOD`

    const MAX_ITERATIONS: usize = 2500;
    const CONVERGENCE_CHECK_PERIOD: usize = 100;
    #[allow(dead_code)]
    const MAX_NODE_DISPLACEMENT: f64 = 300.0;
    const MAX_NODE_DISPLACEMENT_INCREMENTAL: f64 = 100.0; // layout-base `FDLayoutConstants.MAX_NODE_DISPLACEMENT_INCREMENTAL`
    fn from_graph(graph: &Graph) -> Self {
        let leaf_count = graph.nodes.len();
        let mut nodes: Vec<SimNode> = Vec::with_capacity(graph.nodes.len() + graph.compounds.len());
        let mut id_to_idx: FxHashMap<String, usize> = FxHashMap::default();
        id_to_idx.reserve((graph.nodes.len() + graph.compounds.len()).saturating_mul(2));

        for (idx, n) in graph.nodes.iter().enumerate() {
            let w = n.width.max(1.0);
            let h = n.height.max(1.0);
            nodes.push(SimNode {
                id: n.id.clone(),
                parent: n.parent.clone(),
                owner_idx: usize::MAX,
                is_compound: false,
                width: w,
                height: h,
                bounds_extras: n.bounds_extras,
                estimated_size: 0.0,
                left: n.x - w / 2.0,
                top: n.y - h / 2.0,
                spring_fx: 0.0,
                spring_fy: 0.0,
                repulsion_fx: 0.0,
                repulsion_fy: 0.0,
                gravitation_fx: 0.0,
                gravitation_fy: 0.0,
                no_of_children: 1.0,
                padding: 0.0,
                surrounding: Vec::new(),
                grid_start_x: 0,
                grid_finish_x: 0,
                grid_start_y: 0,
                grid_finish_y: 0,
            });
            id_to_idx.insert(n.id.clone(), idx);
        }

        let mut compound_parent: FxHashMap<String, Option<String>> = FxHashMap::default();
        compound_parent.reserve(graph.compounds.len().saturating_mul(2));
        let mut compound_ids_in_order: Vec<String> = Vec::with_capacity(graph.compounds.len());
        for c in &graph.compounds {
            compound_parent.insert(c.id.clone(), c.parent.clone());
            compound_ids_in_order.push(c.id.clone());
        }

        // Materialize compound nodes as layout nodes (Cytoscape parent nodes).
        for c in &graph.compounds {
            let idx = nodes.len();
            nodes.push(SimNode {
                id: c.id.clone(),
                parent: c.parent.clone(),
                owner_idx: usize::MAX,
                is_compound: true,
                width: Self::EMPTY_COMPOUND_NODE_SIZE,
                height: Self::EMPTY_COMPOUND_NODE_SIZE,
                bounds_extras: BoundsExtras::default(),
                estimated_size: Self::EMPTY_COMPOUND_NODE_SIZE,
                left: 0.0,
                top: 0.0,
                spring_fx: 0.0,
                spring_fy: 0.0,
                repulsion_fx: 0.0,
                repulsion_fy: 0.0,
                gravitation_fx: 0.0,
                gravitation_fy: 0.0,
                no_of_children: 1.0,
                padding: 0.0,
                surrounding: Vec::new(),
                grid_start_x: 0,
                grid_finish_x: 0,
                grid_start_y: 0,
                grid_finish_y: 0,
            });
            id_to_idx.insert(c.id.clone(), idx);
        }

        let mut edges: Vec<SimEdge> = Vec::new();
        for e in &graph.edges {
            let Some(&a) = id_to_idx.get(e.source.as_str()) else {
                continue;
            };
            let Some(&b) = id_to_idx.get(e.target.as_str()) else {
                continue;
            };
            if a == b {
                continue;
            }

            let ideal = if e.ideal_length.is_finite() && e.ideal_length > 0.0 {
                e.ideal_length
            } else {
                Self::DEFAULT_EDGE_LENGTH
            };
            let elasticity = if e.elasticity.is_finite() && e.elasticity > 0.0 {
                e.elasticity
            } else {
                Self::DEFAULT_SPRING_STRENGTH
            };
            edges.push(SimEdge {
                a,
                b,
                a_in_lca: a,
                b_in_lca: b,
                a_anchor: e.source_anchor,
                b_anchor: e.target_anchor,
                base_ideal_length: ideal.max(1.0),
                ideal_length: ideal.max(1.0),
                elasticity,
                label_width: e.label_width.filter(|v| v.is_finite() && *v > 0.0),
                label_height: e.label_height.filter(|v| v.is_finite() && *v > 0.0),
            });
        }

        let root_owner_idx = nodes.len();

        // Resolve owner graph identities (`node.getOwner()` in layout-base): nodes repel only
        // within the same owner graph (i.e. same parent compound).
        for n in &mut nodes {
            let owner_idx = n
                .parent
                .as_deref()
                .and_then(|p| id_to_idx.get(p).copied())
                .unwrap_or(root_owner_idx);
            n.owner_idx = owner_idx;
        }

        let mut children_by_owner: Vec<Vec<usize>> = vec![Vec::new(); nodes.len() + 1];
        // Preserve Cytoscape insertion order within each owner graph:
        // - parent (compound) nodes are created before non-parent nodes
        // - within each category, relative order follows Mermaid's `addGroups(...)` and
        //   `addServices/addJunctions(...)` array iteration order
        //
        // This ordering is observable in `graphManager.getGraphs()/getAllNodes()` iteration and
        // affects deterministic parity for FR-grid repulsion (processed set ordering).
        for c in &graph.compounds {
            if let Some(&idx) = id_to_idx.get(c.id.as_str()) {
                let owner = nodes
                    .get(idx)
                    .map(|n| n.owner_idx)
                    .unwrap_or(root_owner_idx);
                if owner < children_by_owner.len() {
                    children_by_owner[owner].push(idx);
                }
            }
        }
        for n in &graph.nodes {
            if let Some(&idx) = id_to_idx.get(n.id.as_str()) {
                let owner = nodes
                    .get(idx)
                    .map(|n| n.owner_idx)
                    .unwrap_or(root_owner_idx);
                if owner < children_by_owner.len() {
                    children_by_owner[owner].push(idx);
                }
            }
        }

        // Compute compound inclusion depths (root-level nodes depth=1), and build a stable
        // deepest-first compound node order for updateBounds.
        let mut inclusion_depth: Vec<usize> = vec![1; nodes.len()];
        fn depth_of(
            idx: usize,
            nodes: &[SimNode],
            memo: &mut [Option<usize>],
            guard: &mut usize,
        ) -> usize {
            if let Some(v) = memo[idx] {
                return v;
            }
            *guard += 1;
            if *guard > 4096 {
                return 1;
            }
            let d = match nodes[idx].owner_idx {
                owner if owner >= nodes.len() => 1,
                parent_idx => depth_of(parent_idx, nodes, memo, guard).saturating_add(1),
            };
            memo[idx] = Some(d);
            d
        }
        let mut memo: Vec<Option<usize>> = vec![None; nodes.len()];
        for i in 0..nodes.len() {
            let mut guard = 0usize;
            inclusion_depth[i] = depth_of(i, &nodes, &mut memo, &mut guard);
        }
        let mut compounds_deep_first: Vec<usize> = nodes
            .iter()
            .enumerate()
            .filter_map(|(idx, n)| n.is_compound.then_some(idx))
            .collect();
        compounds_deep_first.sort_by_key(|&idx| std::cmp::Reverse(inclusion_depth[idx]));

        // Compute `no_of_children` weights and descendant leaf lists.
        let mut descendant_leaves: Vec<Vec<usize>> = vec![Vec::new(); nodes.len()];
        let mut no_of_children: Vec<f64> = vec![1.0; nodes.len()];

        // Initialize leaf descendants for leaves.
        for idx in 0..nodes.len() {
            if !nodes[idx].is_compound {
                descendant_leaves[idx] = vec![idx];
                no_of_children[idx] = 1.0;
            }
        }

        // For compounds, aggregate descendant leaves from immediate children (postorder).
        for &cidx in &compounds_deep_first {
            let children = &children_by_owner[cidx];
            let mut leaves: Vec<usize> = Vec::new();
            for &child in children {
                leaves.extend(descendant_leaves[child].iter().copied());
            }
            leaves.sort_unstable();
            leaves.dedup();
            if leaves.is_empty() {
                descendant_leaves[cidx] = Vec::new();
                no_of_children[cidx] = 1.0;
            } else {
                no_of_children[cidx] = leaves.len() as f64;
                descendant_leaves[cidx] = leaves;
            }
        }

        // Compute estimated sizes (used for gravity ranges, and to match layout-base defaults).
        let mut est_size: Vec<f64> = vec![0.0; nodes.len()];
        for idx in 0..nodes.len() {
            if !nodes[idx].is_compound {
                est_size[idx] = (nodes[idx].width + nodes[idx].height) / 2.0;
            }
        }
        // Deepest-first postorder (children first).
        for &cidx in &compounds_deep_first {
            let children = &children_by_owner[cidx];
            let sum: f64 = children.iter().map(|&ch| est_size[ch]).sum();
            let size = if children.is_empty() {
                Self::EMPTY_COMPOUND_NODE_SIZE
            } else {
                (sum / (children.len() as f64).sqrt()).max(1.0)
            };
            est_size[cidx] = size;
        }
        // layout-base `LNode.calcEstimatedSize()` also sets compound node `rect.width/height` to
        // the estimated size. This is later overwritten by `updateBounds()`, but it affects
        // early spring-embedder iterations (repulsion ranges, smart ideal edge length, etc.).
        for &cidx in &compounds_deep_first {
            let s = est_size[cidx].max(1.0);
            nodes[cidx].width = s;
            nodes[cidx].height = s;
        }

        for idx in 0..nodes.len() {
            nodes[idx].estimated_size = est_size[idx].max(1.0);
        }

        let mut owner_estimated_size: Vec<f64> =
            vec![Self::EMPTY_COMPOUND_NODE_SIZE; nodes.len() + 1];
        // For compound owners, estimated size is the compound node's estimated size.
        for &cidx in &compounds_deep_first {
            owner_estimated_size[cidx] = est_size[cidx].max(1.0);
        }
        // Root owner estimated size is computed from its immediate children.
        {
            let children = &children_by_owner[root_owner_idx];
            let sum: f64 = children.iter().map(|&ch| est_size[ch]).sum();
            owner_estimated_size[root_owner_idx] = if children.is_empty() {
                Self::EMPTY_COMPOUND_NODE_SIZE
            } else {
                (sum / (children.len() as f64).sqrt()).max(1.0)
            };
        }

        for (idx, n) in nodes.iter_mut().enumerate() {
            n.no_of_children = no_of_children[idx].max(1.0);
        }

        Self {
            nodes,
            edges,
            id_to_idx,
            compound_parent,
            compound_ids_in_order,
            leaf_count,
            root_owner_idx,
            children_by_owner,
            compounds_deep_first,
            descendant_leaves,
            owner_estimated_size,
            inclusion_depth,
        }
    }

    fn translate(&mut self, dx: f64, dy: f64) {
        for n in &mut self.nodes {
            n.left += dx;
            n.top += dy;
        }
    }

    fn bounding_box_center_rects(&self) -> Option<(f64, f64)> {
        if self.nodes.is_empty() {
            return None;
        }
        let mut min_x = f64::INFINITY;
        let mut min_y = f64::INFINITY;
        let mut max_x = f64::NEG_INFINITY;
        let mut max_y = f64::NEG_INFINITY;
        for n in &self.nodes {
            min_x = min_x.min(n.left);
            min_y = min_y.min(n.top);
            max_x = max_x.max(n.right());
            max_y = max_y.max(n.bottom());
        }
        if !(min_x.is_finite() && min_y.is_finite() && max_x.is_finite() && max_y.is_finite()) {
            return None;
        }
        Some(((min_x + max_x) / 2.0, (min_y + max_y) / 2.0))
    }

    fn bounding_box_center_nodes_with_extras(&self) -> Option<(f64, f64)> {
        if self.nodes.is_empty() {
            return None;
        }
        let mut min_x = f64::INFINITY;
        let mut min_y = f64::INFINITY;
        let mut max_x = f64::NEG_INFINITY;
        let mut max_y = f64::NEG_INFINITY;
        for n in &self.nodes {
            min_x = min_x.min(n.bound_left());
            min_y = min_y.min(n.bound_top());
            max_x = max_x.max(n.bound_right());
            max_y = max_y.max(n.bound_bottom());
        }
        if !(min_x.is_finite() && min_y.is_finite() && max_x.is_finite() && max_y.is_finite()) {
            return None;
        }
        Some(((min_x + max_x) / 2.0, (min_y + max_y) / 2.0))
    }

    fn bounding_box_center_eles(&self, run_idx: usize) -> Option<(f64, f64)> {
        if self.nodes.is_empty() {
            return None;
        }

        // Cytoscape edge bboxes are inflated by a small padding (see `edge.boundingBox()`).
        // Mermaid Architecture baselines empirically match ~2.5px here.
        const EDGE_BBOX_PAD: f64 = 2.5;
        // Cytoscape compound nodes (with `padding`) end up slightly larger than the naive
        // "child bbox + padding" model, largely due to renderer bbox padding.
        // Mermaid Architecture baselines match adding an additional ~1.5px on each side.
        const COMPOUND_BBOX_EXTRA: f64 = 1.5;

        let mut min_x = f64::INFINITY;
        let mut min_y = f64::INFINITY;
        let mut max_x = f64::NEG_INFINITY;
        let mut max_y = f64::NEG_INFINITY;

        // Nodes (incl. labels). Cytoscape `eles.boundingBox()` treats compound nodes as wrappers
        // around their child graphs, and `compound-sizing-wrt-labels: include` causes descendant
        // label extents to affect compound bounds.
        //
        // Model this by:
        // - using leaf `bound_*` (includes label extras) as the base primitive
        // - computing compound bboxes bottom-up from immediate children (so compound padding
        //   stacks across deep nesting, as observed in Mermaid/Cytoscape)
        // - inflating each compound by `padding + COMPOUND_BBOX_EXTRA`
        //
        // This keeps layout rects (used by the spring embedder) unchanged while making the
        // relocation origin (`eles.boundingBox()` center) match upstream.
        #[derive(Debug, Clone, Copy)]
        struct Bb {
            x1: f64,
            y1: f64,
            x2: f64,
            y2: f64,
        }

        impl Bb {
            fn union(self, other: Bb) -> Bb {
                Bb {
                    x1: self.x1.min(other.x1),
                    y1: self.y1.min(other.y1),
                    x2: self.x2.max(other.x2),
                    y2: self.y2.max(other.y2),
                }
            }

            fn inflate(self, pad: f64) -> Bb {
                Bb {
                    x1: self.x1 - pad,
                    y1: self.y1 - pad,
                    x2: self.x2 + pad,
                    y2: self.y2 + pad,
                }
            }
        }

        fn leaf_bbox(n: &SimNode) -> Bb {
            let x1 = n.bound_left();
            let y1 = n.bound_top();
            let x2 = n.bound_right();
            let y2 = n.bound_bottom();
            Bb { x1, y1, x2, y2 }
        }

        let mut bbs: Vec<Option<Bb>> = vec![None; self.nodes.len()];
        for (idx, n) in self.nodes.iter().enumerate() {
            if !n.is_compound {
                bbs[idx] = Some(leaf_bbox(n));
            }
        }

        for &cidx in &self.compounds_deep_first {
            let Some(n) = self.nodes.get(cidx) else {
                continue;
            };
            if !n.is_compound {
                continue;
            }
            let children = self
                .children_by_owner
                .get(cidx)
                .map(|v| v.as_slice())
                .unwrap_or(&[]);
            if children.is_empty() {
                // Empty compound: fall back to its rect (no label extras tracked for compounds).
                bbs[cidx] = Some(Bb {
                    x1: n.left,
                    y1: n.top,
                    x2: n.right(),
                    y2: n.bottom(),
                });
                continue;
            }
            let mut bb: Option<Bb> = None;
            for &ch in children {
                let ch_bb = bbs.get(ch).and_then(|v| *v).unwrap_or_else(|| {
                    let Some(cn) = self.nodes.get(ch) else {
                        return Bb {
                            x1: 0.0,
                            y1: 0.0,
                            x2: 0.0,
                            y2: 0.0,
                        };
                    };
                    Bb {
                        x1: cn.left,
                        y1: cn.top,
                        x2: cn.right(),
                        y2: cn.bottom(),
                    }
                });
                bb = Some(bb.map(|b| b.union(ch_bb)).unwrap_or(ch_bb));
            }
            let pad = n.padding.max(0.0) + COMPOUND_BBOX_EXTRA;
            bbs[cidx] = bb.map(|b| b.inflate(pad));
        }

        let top_level = self
            .children_by_owner
            .get(self.root_owner_idx)
            .map(|v| v.as_slice())
            .unwrap_or(&[]);
        for &idx in top_level {
            let Some(bb) = bbs.get(idx).and_then(|v| *v).or_else(|| {
                self.nodes.get(idx).map(|n| Bb {
                    x1: n.left,
                    y1: n.top,
                    x2: n.right(),
                    y2: n.bottom(),
                })
            }) else {
                continue;
            };
            min_x = min_x.min(bb.x1);
            min_y = min_y.min(bb.y1);
            max_x = max_x.max(bb.x2);
            max_y = max_y.max(bb.y2);
        }

        // Edges: Cytoscape `eles.boundingBox()` includes edge geometry. For Mermaid Architecture,
        // edge endpoints are manually specified as `{ 0/50%/100% }` offsets (see
        // `source-endpoint`/`target-endpoint` in Mermaid's Cytoscape stylesheet).
        //
        // Cytoscape resolves those endpoints by adding the pixel offset to `node.position()`
        // (even though `position()` represents the node center for normal geometry). Mermaid then
        // reuses that same `position()` value as the SVG top-left `translate(x,y)` for the icon.
        //
        // In our port we mirror upstream by treating `SimNode.center_{x,y}` as that top-left
        // anchor, and compute endpoint points as offsets from it (not as shape intersections).
        fn endpoint(n: &SimNode, anchor: Option<Anchor>) -> (f64, f64) {
            let ox = n.center_x();
            let oy = n.center_y();
            let w = n.width;
            let h = n.height;
            match anchor {
                Some(Anchor::Left) => (ox, oy + (h / 2.0)),
                Some(Anchor::Right) => (ox + w, oy + (h / 2.0)),
                Some(Anchor::Top) => (ox + (w / 2.0), oy),
                Some(Anchor::Bottom) => (ox + (w / 2.0), oy + h),
                None => (ox + (w / 2.0), oy + (h / 2.0)),
            }
        }

        fn polyline_midpoint(points: &[(f64, f64)]) -> Option<(f64, f64)> {
            if points.len() < 2 {
                return None;
            }
            let mut total = 0.0f64;
            let mut seg_lens: Vec<f64> = Vec::with_capacity(points.len().saturating_sub(1));
            for w in points.windows(2) {
                let dx = w[1].0 - w[0].0;
                let dy = w[1].1 - w[0].1;
                let len = (dx * dx + dy * dy).sqrt();
                seg_lens.push(len);
                total += len;
            }
            if !total.is_finite() || total <= 0.0 {
                return Some(points[0]);
            }
            let target = total / 2.0;
            let mut acc = 0.0f64;
            for (i, &len) in seg_lens.iter().enumerate() {
                if !len.is_finite() || len <= 0.0 {
                    continue;
                }
                if acc + len >= target {
                    let t = ((target - acc) / len).clamp(0.0, 1.0);
                    let (x0, y0) = points[i];
                    let (x1, y1) = points[i + 1];
                    return Some((x0 + (x1 - x0) * t, y0 + (y1 - y0) * t));
                }
                acc += len;
            }
            points.last().copied()
        }

        for e in &self.edges {
            let Some(a) = self.nodes.get(e.a) else {
                continue;
            };
            let Some(b) = self.nodes.get(e.b) else {
                continue;
            };
            let (sx, sy) = endpoint(a, e.a_anchor);
            let (tx, ty) = endpoint(b, e.b_anchor);

            let mut path_points: Vec<(f64, f64)> = vec![(sx, sy), (tx, ty)];
            let mut label_point_override: Option<(f64, f64)> = None;

            fn include_point(
                min_x: &mut f64,
                min_y: &mut f64,
                max_x: &mut f64,
                max_y: &mut f64,
                x: f64,
                y: f64,
            ) {
                *min_x = (*min_x).min(x);
                *min_y = (*min_y).min(y);
                *max_x = (*max_x).max(x);
                *max_y = (*max_y).max(y);
            }

            include_point(
                &mut min_x,
                &mut min_y,
                &mut max_x,
                &mut max_y,
                sx - EDGE_BBOX_PAD,
                sy - EDGE_BBOX_PAD,
            );
            include_point(
                &mut min_x,
                &mut min_y,
                &mut max_x,
                &mut max_y,
                sx + EDGE_BBOX_PAD,
                sy + EDGE_BBOX_PAD,
            );
            include_point(
                &mut min_x,
                &mut min_y,
                &mut max_x,
                &mut max_y,
                tx - EDGE_BBOX_PAD,
                ty - EDGE_BBOX_PAD,
            );
            include_point(
                &mut min_x,
                &mut min_y,
                &mut max_x,
                &mut max_y,
                tx + EDGE_BBOX_PAD,
                ty + EDGE_BBOX_PAD,
            );

            // Mermaid styles diagonal (XY) edges as Cytoscape `curve-style: segments` with
            // `segment-weights: 0` and `segment-distances: 0.5px` in the pre-layout state. This
            // creates a subtle asymmetric bbox (e.g. ±0.353553... for 45deg edges).
            if run_idx == 0 && sx != tx && sy != ty {
                const SEG_DIST: f64 = 0.5;
                let dx = tx - sx;
                let dy = ty - sy;
                let len = (dx * dx + dy * dy).sqrt();
                if len.is_finite() && len > 0.0 {
                    // Left-hand perpendicular, normalized.
                    let off_x = (-dy / len) * SEG_DIST;
                    let off_y = (dx / len) * SEG_DIST;
                    // `segment-weights: 0` => base point at source endpoint.
                    let px = sx + off_x;
                    let py = sy + off_y;
                    path_points.insert(1, (px, py));
                    // Cytoscape's pre-layout `segments` curve places the edge label near the
                    // segment control point. Using that point for bbox purposes matches the
                    // upstream `edge.boundingBox()` extents for diagonal Architecture edges.
                    label_point_override = Some((px, py));
                    include_point(
                        &mut min_x,
                        &mut min_y,
                        &mut max_x,
                        &mut max_y,
                        px - EDGE_BBOX_PAD,
                        py - EDGE_BBOX_PAD,
                    );
                    include_point(
                        &mut min_x,
                        &mut min_y,
                        &mut max_x,
                        &mut max_y,
                        px + EDGE_BBOX_PAD,
                        py + EDGE_BBOX_PAD,
                    );
                }
            }

            // After the first run, Mermaid updates segment weights/distances so diagonal edges
            // become orthogonal with a single bend at either `(sx, ty)` or `(tx, sy)` depending on
            // the source direction.
            if run_idx > 0 && sx != tx && sy != ty {
                let (bx, by) = match e.a_anchor {
                    Some(Anchor::Top) | Some(Anchor::Bottom) => (sx, ty),
                    _ => (tx, sy),
                };
                path_points.insert(1, (bx, by));
                // After Mermaid updates segment weights/distances, the SVG renderer treats this
                // bend point as the "midpoint" of the orthogonal polyline. Cytoscape's edge label
                // placement for the same style is closest to this bend as well, so use it for the
                // bbox approximation.
                label_point_override = Some((bx, by));
                include_point(
                    &mut min_x,
                    &mut min_y,
                    &mut max_x,
                    &mut max_y,
                    bx - EDGE_BBOX_PAD,
                    by - EDGE_BBOX_PAD,
                );
                include_point(
                    &mut min_x,
                    &mut min_y,
                    &mut max_x,
                    &mut max_y,
                    bx + EDGE_BBOX_PAD,
                    by + EDGE_BBOX_PAD,
                );
            }

            // Edge labels: Cytoscape includes label geometry inside `edge.boundingBox()`, and
            // `eles.boundingBox()` unions it into the overall component bbox.
            if let (Some(lw), Some(lh)) = (e.label_width, e.label_height) {
                let lw = lw.max(0.0);
                let lh = lh.max(0.0);
                if lw.is_finite() && lw > 0.0 && lh.is_finite() && lh > 0.0 {
                    let mp = label_point_override.or_else(|| polyline_midpoint(&path_points));
                    if let Some((mx, my)) = mp {
                        let hw = lw / 2.0;
                        let hh = lh / 2.0;
                        include_point(
                            &mut min_x,
                            &mut min_y,
                            &mut max_x,
                            &mut max_y,
                            mx - hw - EDGE_BBOX_PAD,
                            my - hh - EDGE_BBOX_PAD,
                        );
                        include_point(
                            &mut min_x,
                            &mut min_y,
                            &mut max_x,
                            &mut max_y,
                            mx + hw + EDGE_BBOX_PAD,
                            my + hh + EDGE_BBOX_PAD,
                        );
                    }
                }
            }
        }

        if !(min_x.is_finite() && min_y.is_finite() && max_x.is_finite() && max_y.is_finite()) {
            return None;
        }
        Some(((min_x + max_x) / 2.0, (min_y + max_y) / 2.0))
    }

    fn bounding_box_center_leaves(&self) -> Option<(f64, f64)> {
        if self.leaf_count == 0 || self.nodes.is_empty() {
            return None;
        }
        let end = self.leaf_count.min(self.nodes.len());
        let mut min_x = f64::INFINITY;
        let mut min_y = f64::INFINITY;
        let mut max_x = f64::NEG_INFINITY;
        let mut max_y = f64::NEG_INFINITY;
        for n in &self.nodes[..end] {
            min_x = min_x.min(n.left);
            min_y = min_y.min(n.top);
            max_x = max_x.max(n.right());
            max_y = max_y.max(n.bottom());
        }
        if !(min_x.is_finite() && min_y.is_finite() && max_x.is_finite() && max_y.is_finite()) {
            return None;
        }
        Some(((min_x + max_x) / 2.0, (min_y + max_y) / 2.0))
    }

    fn update_bounds(&mut self) -> OwnerBounds {
        debug_assert_eq!(self.root_owner_idx, self.nodes.len());

        let owner_count = self.nodes.len() + 1;
        let mut left: Vec<f64> = vec![f64::INFINITY; owner_count];
        let mut right: Vec<f64> = vec![f64::NEG_INFINITY; owner_count];
        let mut top: Vec<f64> = vec![f64::INFINITY; owner_count];
        let mut bottom: Vec<f64> = vec![f64::NEG_INFINITY; owner_count];

        // Mirror layout-base `graphManager.updateBounds()`:
        // - update child compound bounds first
        // - then compute each graph bounds with a margin derived from parent compound padding
        for &cidx in &self.compounds_deep_first {
            let children = &self.children_by_owner[cidx];
            if children.is_empty() {
                // Empty compound: keep its current rect as-is.
                left[cidx] = self.nodes[cidx].left - self.nodes[cidx].padding;
                right[cidx] = self.nodes[cidx].right() + self.nodes[cidx].padding;
                top[cidx] = self.nodes[cidx].top - self.nodes[cidx].padding;
                bottom[cidx] = self.nodes[cidx].bottom() + self.nodes[cidx].padding;
                continue;
            }

            let mut min_x = f64::INFINITY;
            let mut min_y = f64::INFINITY;
            let mut max_x = f64::NEG_INFINITY;
            let mut max_y = f64::NEG_INFINITY;
            for &ch in children {
                min_x = min_x.min(self.nodes[ch].left);
                min_y = min_y.min(self.nodes[ch].top);
                max_x = max_x.max(self.nodes[ch].right());
                max_y = max_y.max(self.nodes[ch].bottom());
            }
            let margin = self.nodes[cidx].padding.max(0.0);
            min_x -= margin;
            min_y -= margin;
            max_x += margin;
            max_y += margin;

            left[cidx] = min_x;
            right[cidx] = max_x;
            top[cidx] = min_y;
            bottom[cidx] = max_y;

            // Update compound node rect to wrap its child graph.
            self.nodes[cidx].left = min_x;
            self.nodes[cidx].top = min_y;
            self.nodes[cidx].width = (max_x - min_x).max(1.0);
            self.nodes[cidx].height = (max_y - min_y).max(1.0);
        }

        // Root graph bounds (margin defaults to LayoutConstants.DEFAULT_GRAPH_MARGIN).
        {
            let children = &self.children_by_owner[self.root_owner_idx];
            if children.is_empty() {
                left[self.root_owner_idx] = 0.0;
                right[self.root_owner_idx] = 0.0;
                top[self.root_owner_idx] = 0.0;
                bottom[self.root_owner_idx] = 0.0;
            } else {
                let mut min_x = f64::INFINITY;
                let mut min_y = f64::INFINITY;
                let mut max_x = f64::NEG_INFINITY;
                let mut max_y = f64::NEG_INFINITY;
                for &ch in children {
                    min_x = min_x.min(self.nodes[ch].left);
                    min_y = min_y.min(self.nodes[ch].top);
                    max_x = max_x.max(self.nodes[ch].right());
                    max_y = max_y.max(self.nodes[ch].bottom());
                }
                let margin = Self::DEFAULT_GRAPH_MARGIN;
                left[self.root_owner_idx] = min_x - margin;
                right[self.root_owner_idx] = max_x + margin;
                top[self.root_owner_idx] = min_y - margin;
                bottom[self.root_owner_idx] = max_y + margin;
            }
        }

        OwnerBounds {
            left,
            right,
            top,
            bottom,
        }
    }

    fn max_compound_nesting_depth(&self) -> usize {
        if self.compound_parent.is_empty() {
            return 0;
        }

        let mut max_depth: usize = 0;
        for n in &self.nodes {
            let mut depth: usize = 0;
            let mut cur = n.parent.as_deref();
            let mut guard: usize = 0;
            while let Some(id) = cur {
                depth += 1;
                max_depth = max_depth.max(depth);
                guard += 1;
                if guard > 1024 {
                    break;
                }
                cur = self.compound_parent.get(id).and_then(|p| p.as_deref());
            }
        }
        max_depth
    }

    fn all_nodes_layout_order(&self) -> Vec<usize> {
        // layout-base `graphManager.getAllNodes()` returns a flat list created by concatenating
        // `graph.getNodes()` over `graphManager.getGraphs()` in graph creation order. Graphs are
        // created recursively: root graph first, then each compound's child graph when that
        // compound node is encountered.
        //
        // Reconstruct that order by visiting owner graphs in pre-order, following the
        // `children_by_owner` inclusion tree.
        let mut out: Vec<usize> = Vec::with_capacity(self.nodes.len());
        let mut visited_graph: Vec<bool> = vec![false; self.nodes.len() + 1];

        fn visit(owner: usize, sim: &SimGraph, out: &mut Vec<usize>, visited: &mut [bool]) {
            if owner >= visited.len() {
                return;
            }
            if visited[owner] {
                return;
            }
            visited[owner] = true;

            let nodes = sim
                .children_by_owner
                .get(owner)
                .map(|v| v.as_slice())
                .unwrap_or(&[]);
            for &idx in nodes {
                out.push(idx);
            }
            for &idx in nodes {
                let is_compound = sim.nodes.get(idx).is_some_and(|n| n.is_compound);
                let has_children = sim
                    .children_by_owner
                    .get(idx)
                    .is_some_and(|v| !v.is_empty());
                if is_compound && has_children {
                    visit(idx, sim, out, visited);
                }
            }
        }

        visit(self.root_owner_idx, self, &mut out, &mut visited_graph);
        out
    }

    fn run_spring_embedder(
        &mut self,
        constraints: &Constraints,
        opts: &FcoseOptions,
        rng: &mut XorShift64Star,
        mut timings: Option<&mut FcoseSpringTimings>,
    ) -> SpringStats {
        if self.nodes.is_empty() {
            return SpringStats::default();
        }

        let timing_enabled = timings.is_some();
        let debug_positions = std::env::var("MANATEE_FCOSE_DEBUG_POSITIONS")
            .ok()
            .as_deref()
            == Some("1");
        let debug_positions_all = std::env::var("MANATEE_FCOSE_DEBUG_POSITIONS_ALL")
            .ok()
            .as_deref()
            == Some("1");
        let leaf_count = self.leaf_count;
        let dump_positions = |tag: &str, nodes: &[SimNode]| {
            if !debug_positions {
                return;
            }
            // Default to leaf nodes for readability; optionally include compounds.
            let n = if debug_positions_all {
                nodes.len()
            } else {
                leaf_count.min(nodes.len())
            };
            eprintln!("[manatee-fcose-pos] {tag} leaf_count={n}");
            for i in 0..n {
                eprintln!(
                    "[manatee-fcose-pos] {tag} id={} compound={} center=({:.6},{:.6}) left_top=({:.6},{:.6}) size=({:.3},{:.3}) owner={}",
                    nodes[i].id,
                    nodes[i].is_compound,
                    nodes[i].center_x(),
                    nodes[i].center_y(),
                    nodes[i].left,
                    nodes[i].top,
                    nodes[i].width,
                    nodes[i].height,
                    nodes[i].owner_idx
                );
            }
        };

        let opts_prep_start = timing_enabled.then(std::time::Instant::now);
        // Recompute per-edge ideal lengths (layout-base `FDLayout.calcIdealEdgeLengths`).
        // This must be re-applied on each run because Mermaid runs FCoSE twice.
        //
        // Important: the "global" default edge length constant (used for several heuristics) is
        // derived from the *base* `idealEdgeLength` option before the smart inter-graph
        // adjustments are applied. Reset first so the second run starts from the same baseline.
        self.reset_edge_ideal_lengths();

        // layout-base/CoSE uses a *global* `DEFAULT_EDGE_LENGTH` for multiple heuristics (minimum
        // repulsion distance, overlap separation buffer, repulsion grid range, convergence
        // thresholds, etc.). In upstream Cytoscape FCoSE this value is derived from the
        // `idealEdgeLength` option (before per-edge nesting/smart adjustments).
        let default_edge_length = opts
            .default_edge_length
            .filter(|v| v.is_finite() && *v > 0.0)
            .unwrap_or_else(|| {
                if self.edges.is_empty() {
                    Self::DEFAULT_EDGE_LENGTH
                } else {
                    let sum: f64 = self.edges.iter().map(|e| e.ideal_length).sum();
                    (sum / (self.edges.len() as f64)).max(1.0)
                }
            });
        let half_default_edge_length = default_edge_length / 2.0;

        self.adjust_intergraph_ideal_edge_lengths();
        if std::env::var("MANATEE_FCOSE_DEBUG_EDGE_LENGTHS")
            .ok()
            .as_deref()
            == Some("1")
        {
            let mut inter = 0usize;
            for e in &self.edges {
                let intergraph = self.nodes[e.a].owner_idx != self.nodes[e.b].owner_idx;
                if intergraph {
                    inter += 1;
                }
                eprintln!(
                    "[manatee-fcose-edge] a={} b={} inter={} base={:.6} ideal={:.6} elasticity={:.6}",
                    e.a, e.b, intergraph, e.base_ideal_length, e.ideal_length, e.elasticity
                );
            }
            eprintln!(
                "[manatee-fcose-edge] edges={} intergraph={}",
                self.edges.len(),
                inter
            );
        }
        // CoSE updates `MIN_REPULSION_DIST` based on the effective `DEFAULT_EDGE_LENGTH` when
        // `idealEdgeLength` is set. For Mermaid Architecture this is always set (as a function),
        // so we scale the minimum repulsion distance with the average ideal length.
        let min_repulsion_dist = (default_edge_length / 10.0).max(0.0005);
        if let (Some(t), Some(s)) = (timings.as_deref_mut(), opts_prep_start) {
            t.opts_prep = s.elapsed();
        }

        // Apply uniform compound padding (Cytoscape style `padding`).
        let compound_padding = opts.compound_padding.unwrap_or(0.0).max(0.0);
        for n in &mut self.nodes {
            if n.is_compound {
                n.padding = compound_padding;
            }
        }

        // FCoSE performs a spectral initialization when `randomize=true` (Mermaid defaults to
        // `randomize: true`). The upstream JS implementation relies on `Math.random`; in Rust we
        // make this explicit and deterministic via `random_seed`.
        let spectral_start = timing_enabled.then(std::time::Instant::now);
        let spectral_edges: Vec<SimEdge> = self
            .edges
            .iter()
            .copied()
            .filter(|e| e.a < self.leaf_count && e.b < self.leaf_count)
            .collect();
        let spectral_applied = spectral::apply_spectral_start_positions(
            &mut self.nodes[..self.leaf_count],
            &spectral_edges,
            &self.compound_parent,
            &self.compound_ids_in_order,
            rng,
        );
        dump_positions("spectral_raw", &self.nodes);
        if let (Some(t), Some(s)) = (timings.as_deref_mut(), spectral_start) {
            t.spectral = s.elapsed();
        }

        dump_positions("spectral", &self.nodes);

        let gravity_constant = Self::DEFAULT_GRAVITY_STRENGTH;
        let debug_forces = std::env::var("MANATEE_FCOSE_DEBUG_FORCES").ok().as_deref() == Some("1");
        let debug_edge_forces = std::env::var("MANATEE_FCOSE_DEBUG_EDGE_FORCES")
            .ok()
            .as_deref()
            == Some("1");
        let debug_disps = std::env::var("MANATEE_FCOSE_DEBUG_DISPS").ok().as_deref() == Some("1");

        // Upstream CoSE applies gravitational forces only to nodes that belong to a *disconnected*
        // owner graph (see `CoSELayout.calculateNodesToApplyGravitationTo()`), where "connected"
        // is evaluated in the compound graph using `LEdge.getOtherEndInGraph(...)` (i.e. edges
        // incident to descendants can connect immediate children via ancestor lifting).
        //
        // Applying gravity to all nodes (a common simplification) makes sparse cross-compound
        // Architecture graphs significantly more compact than Mermaid/Cytoscape, which in turn
        // changes the root `viewBox/max-width` in parity-root comparisons.
        let disable_gravity = std::env::var("MANATEE_FCOSE_DISABLE_GRAVITY")
            .ok()
            .as_deref()
            == Some("1");
        let apply_gravity: Vec<bool> = if disable_gravity {
            vec![false; self.nodes.len()]
        } else {
            self.nodes_apply_gravity_mask()
        };

        // Match `cose-base` repulsion cutoff (`CoSELayout.calcRepulsionRange()`):
        //
        // `repulsionRange = 2 * (level + 1) * idealEdgeLength`
        //
        // Cytoscape FCoSE runs the compound graph in a single CoSE pass with `level=0`, so this
        // reduces to `2 * idealEdgeLength`.
        let repulsion_range = (2.0 * default_edge_length).max(1.0);

        // layout-base uses the FR-grid repulsion variant by default, which caches each node's
        // surrounding set and refreshes it every `GRID_CALCULATION_CHECK_PERIOD` iterations.
        let mut repulsion_grid: Option<RepulsionGrid> = None;

        // Fallback for degenerate cases where spectral is skipped (e.g. very small graphs).
        if self.edges.is_empty() && !spectral_applied {
            let collapse_start = timing_enabled.then(std::time::Instant::now);
            self.collapse_start_positions(default_edge_length, rng);
            if let (Some(t), Some(s)) = (timings.as_deref_mut(), collapse_start) {
                t.collapse_start_positions = s.elapsed();
            }
        }

        // Upstream `cose-base` runs a dedicated constraint handler before the spring embedder.
        // This can rotate/reflect the draft layout and enforce alignment/relative-placement
        // constraints in position space, which strongly affects overall orientation and the
        // parity-root root viewport.
        let disable_pre = std::env::var("MANATEE_FCOSE_DISABLE_PRE_CONSTRAINTS")
            .ok()
            .as_deref()
            == Some("1");
        if !disable_pre
            && !(constraints.align_horizontal.is_empty()
                && constraints.align_vertical.is_empty()
                && constraints.relative.is_empty())
        {
            let pre_constraints_start = timing_enabled.then(std::time::Instant::now);
            handle_constraints_pre_layout(&mut self.nodes[..self.leaf_count], constraints);
            dump_positions("pre_constraints", &self.nodes);
            if let (Some(t), Some(s)) = (timings.as_deref_mut(), pre_constraints_start) {
                t.pre_constraints = s.elapsed();
            }
        }

        let constraint_rt_start = timing_enabled.then(std::time::Instant::now);
        let mut constraint_rt = ConstraintRuntime::new(&self.nodes, constraints);
        if let (Some(t), Some(s)) = (timings.as_deref_mut(), constraint_rt_start) {
            t.constraint_rt = s.elapsed();
        }

        let n = self.nodes.len() as f64;
        let displacement_threshold_per_node = (3.0 * default_edge_length) / 100.0;
        let total_displacement_threshold = displacement_threshold_per_node * n;

        // cytoscape-fcose postprocessing (`cose.js`) forces CoSE incremental mode on by setting
        // `LayoutConstants.DEFAULT_INCREMENTAL = true`. This means we start with the incremental
        // cooling factor and max displacement values, even when `randomize=true`.
        //
        // This is a major contributor to parity-root `viewBox/max-width` stability for sparse
        // graphs (notably the Architecture fixtures).
        let initial_cooling_factor = Self::DEFAULT_COOLING_FACTOR_INCREMENTAL;
        let mut cooling_factor = initial_cooling_factor;
        let max_node_displacement = Self::MAX_NODE_DISPLACEMENT_INCREMENTAL;
        let max_iterations = Self::MAX_ITERATIONS.max(self.nodes.len() * 5);
        let max_cooling_cycle = (max_iterations as f64) / (Self::CONVERGENCE_CHECK_PERIOD as f64);
        let final_temperature = Self::FINAL_TEMPERATURE;
        let mut cooling_cycle = 0.0f64;

        let mut total_iterations = 0usize;
        let mut old_total_displacement = 0.0f64;
        let mut last_total_displacement = 0.0f64;

        let disable_convergence = std::env::var("MANATEE_FCOSE_DISABLE_CONVERGENCE")
            .ok()
            .as_deref()
            == Some("1");

        let iterations_start = timing_enabled.then(std::time::Instant::now);
        let mut processed: Vec<bool> = vec![false; self.nodes.len()];
        let mut disps: Vec<(f64, f64)> = vec![(0.0, 0.0); self.nodes.len()];
        let all_nodes_in_layout_order = self.all_nodes_layout_order();
        loop {
            total_iterations += 1;
            if total_iterations == max_iterations {
                break;
            }

            if total_iterations.is_multiple_of(Self::CONVERGENCE_CHECK_PERIOD) {
                let oscilating = total_iterations > (max_iterations / 3)
                    && (last_total_displacement - old_total_displacement).abs() < 2.0;
                let converged = last_total_displacement < total_displacement_threshold;

                old_total_displacement = last_total_displacement;

                if !disable_convergence && (converged || oscilating) {
                    break;
                }

                cooling_cycle += 1.0;

                let numerator = (100.0 * (initial_cooling_factor - final_temperature)).ln();
                let denominator = max_cooling_cycle.ln().max(1e-9);
                let power = numerator / denominator;
                let schedule = cooling_cycle.powf(power) / 100.0;
                cooling_factor = (initial_cooling_factor - schedule).max(final_temperature);
            }

            let mut total_displacement = 0.0f64;

            // Match `cose-base` tick order: update compound bounds (with padding) before forces.
            let bounds = self.update_bounds();

            // Spring forces (per-edge ideal lengths).
            for e in &self.edges {
                // layout-base spring forces act between the edge's actual endpoints
                // (`edge.getSource()/getTarget()`), even for inter-graph edges. LCA-lifted
                // endpoints are used for *ideal edge length* adjustments, not for force
                // application.
                let (a, b) = (e.a, e.b);
                if a == b {
                    continue;
                }
                if a >= self.nodes.len() || b >= self.nodes.len() {
                    continue;
                }
                if rects_intersect(&self.nodes[a], &self.nodes[b]) {
                    continue;
                }
                let (ax, ay, bx, by) = rect_clip_points(&self.nodes[a], &self.nodes[b]);
                let mut lx = bx - ax;
                let mut ly = by - ay;
                let raw_lx = lx;
                let raw_ly = ly;
                // layout-base `LEdge.updateLength()` clamps very small deltas to {-1, 0, 1}
                // using `IMath.sign`, to avoid divide-by-zero instability.
                if lx.abs() < 1.0 {
                    lx = imath_sign(lx);
                }
                if ly.abs() < 1.0 {
                    ly = imath_sign(ly);
                }
                let len = (lx * lx + ly * ly).sqrt();
                if len == 0.0 {
                    continue;
                }

                // In Cytoscape CoSE/FCoSE, the spring force is scaled by the effective
                // `edgeElasticity` option. Mermaid Architecture sets this to `0.45` for
                // same-parent edges and `0.001` for edges that cross a group boundary.
                let spring_force = e.elasticity * (len - e.ideal_length.max(1.0));
                let sfx = spring_force * (lx / len);
                let sfy = spring_force * (ly / len);
                if debug_edge_forces && total_iterations == 1 {
                    let ida = self.nodes.get(a).map(|n| n.id.as_str()).unwrap_or("<oob>");
                    let idb = self.nodes.get(b).map(|n| n.id.as_str()).unwrap_or("<oob>");
                    eprintln!(
                        "[manatee-fcose-edge-force] iter=1 a={} b={} a_ctr=({:.15},{:.15}) b_ctr=({:.15},{:.15}) clip_a=({:.15},{:.15}) clip_b=({:.15},{:.15}) raw_lx={:.20} raw_ly={:.20} raw_ly_bits={:#018x} lx={:.12} ly={:.12} len={:.12} ideal={:.12} elast={:.6} spring={:.12} sfx={:.12} sfy={:.12}",
                        ida,
                        idb,
                        self.nodes[a].center_x(),
                        self.nodes[a].center_y(),
                        self.nodes[b].center_x(),
                        self.nodes[b].center_y(),
                        ax,
                        ay,
                        bx,
                        by,
                        raw_lx,
                        raw_ly,
                        raw_ly.to_bits(),
                        lx,
                        ly,
                        len,
                        e.ideal_length.max(1.0),
                        e.elasticity,
                        spring_force,
                        sfx,
                        sfy
                    );
                }
                self.nodes[a].spring_fx += sfx;
                self.nodes[a].spring_fy += sfy;
                self.nodes[b].spring_fx -= sfx;
                self.nodes[b].spring_fy -= sfy;
            }

            // Repulsion forces (layout-base FR grid variant, with cached surrounding lists).
            //
            // Upstream refreshes the grid + surrounding lists when `totalIterations % 10 == 1`,
            // then reuses those "stale" surrounding lists for the next 9 iterations.
            let refresh_surrounding = (total_iterations % Self::GRID_CALCULATION_CHECK_PERIOD) == 1;
            if refresh_surrounding {
                let (l, r, t, b) = (
                    bounds.left[self.root_owner_idx],
                    bounds.right[self.root_owner_idx],
                    bounds.top[self.root_owner_idx],
                    bounds.bottom[self.root_owner_idx],
                );
                repulsion_grid = RepulsionGrid::build(
                    l,
                    t,
                    r,
                    b,
                    &mut self.nodes,
                    repulsion_range,
                    &all_nodes_in_layout_order,
                );
            }

            if repulsion_range.is_finite() && repulsion_range > 0.0 {
                processed.fill(false);
                for &i in &all_nodes_in_layout_order {
                    if i >= self.nodes.len() {
                        continue;
                    }
                    if refresh_surrounding {
                        if let Some(g) = &repulsion_grid {
                            g.refresh_node_surrounding(
                                i,
                                &mut self.nodes,
                                &processed,
                                repulsion_range,
                            );
                        } else {
                            self.nodes[i].surrounding.clear();
                        }
                    }

                    let surrounding = std::mem::take(&mut self.nodes[i].surrounding);
                    for &j in &surrounding {
                        if i == j {
                            continue;
                        }
                        if j >= self.nodes.len() {
                            continue;
                        }
                        if self.nodes[i].owner_idx != self.nodes[j].owner_idx {
                            continue;
                        }
                        let (rfx, rfy) = calc_repulsion_force(
                            &self.nodes[i],
                            &self.nodes[j],
                            min_repulsion_dist,
                            half_default_edge_length,
                        );
                        // Apply a symmetric pairwise force.
                        //
                        // Unlike `i < j` index-based deduping, upstream CoSE/FCoSE dedupes via a
                        // processed set in `getAllNodes()` order. Here `surrounding` is already
                        // filtered by `processed`, so we must not skip pairs where `j < i`.
                        self.nodes[i].repulsion_fx += rfx;
                        self.nodes[i].repulsion_fy += rfy;
                        self.nodes[j].repulsion_fx -= rfx;
                        self.nodes[j].repulsion_fy -= rfy;
                    }
                    self.nodes[i].surrounding = surrounding;
                    processed[i] = true;
                }
            } else {
                // Fallback: unbounded repulsion (all pairs).
                for i in 0..self.nodes.len() {
                    for j in (i + 1)..self.nodes.len() {
                        if self.nodes[i].owner_idx != self.nodes[j].owner_idx {
                            continue;
                        }
                        let (rfx, rfy) = calc_repulsion_force(
                            &self.nodes[i],
                            &self.nodes[j],
                            min_repulsion_dist,
                            half_default_edge_length,
                        );
                        self.nodes[i].repulsion_fx += rfx;
                        self.nodes[i].repulsion_fy += rfy;
                        self.nodes[j].repulsion_fx -= rfx;
                        self.nodes[j].repulsion_fy -= rfy;
                    }
                }
            }

            // Gravity forces (layout-base `FDLayout.calcGravitationalForce`), per owner graph.
            for (idx, n) in self.nodes.iter_mut().enumerate() {
                n.gravitation_fx = 0.0;
                n.gravitation_fy = 0.0;
                if !apply_gravity.get(idx).copied().unwrap_or(false) {
                    continue;
                }

                let owner = n.owner_idx;
                let (l, r, t, b) = (
                    bounds.left.get(owner).copied().unwrap_or(0.0),
                    bounds.right.get(owner).copied().unwrap_or(0.0),
                    bounds.top.get(owner).copied().unwrap_or(0.0),
                    bounds.bottom.get(owner).copied().unwrap_or(0.0),
                );
                if !(l.is_finite() && r.is_finite() && t.is_finite() && b.is_finite()) {
                    continue;
                }
                let cx = (l + r) / 2.0;
                let cy = (t + b) / 2.0;

                let dx = n.center_x() - cx;
                let dy = n.center_y() - cy;
                let abs_dx = dx.abs() + n.half_w();
                let abs_dy = dy.abs() + n.half_h();

                let (range_factor, compound_mul) = if owner == self.root_owner_idx {
                    (Self::DEFAULT_GRAVITY_RANGE_FACTOR, 1.0)
                } else {
                    (
                        Self::DEFAULT_COMPOUND_GRAVITY_RANGE_FACTOR,
                        Self::DEFAULT_COMPOUND_GRAVITY_STRENGTH,
                    )
                };
                let estimated =
                    self.owner_estimated_size.get(owner).copied().unwrap_or(0.0) * range_factor;
                if estimated.is_finite() && estimated > 0.0 {
                    if abs_dx > estimated || abs_dy > estimated {
                        n.gravitation_fx = -gravity_constant * dx * compound_mul;
                        n.gravitation_fy = -gravity_constant * dy * compound_mul;
                    }
                }
            }

            if debug_forces && total_iterations == 1 {
                eprintln!("[manatee-fcose-force] iter=1 nodes={}", self.nodes.len());
                for (idx, n) in self.nodes.iter().enumerate() {
                    eprintln!(
                        "[manatee-fcose-force] idx={} id={} owner={} compound={} spring=({:.6},{:.6}) rep=({:.6},{:.6}) grav=({:.6},{:.6})",
                        idx,
                        n.id,
                        n.owner_idx,
                        n.is_compound,
                        n.spring_fx,
                        n.spring_fy,
                        n.repulsion_fx,
                        n.repulsion_fy,
                        n.gravitation_fx,
                        n.gravitation_fy,
                    );
                }
            }

            // Move nodes (with constraints applied to displacements).
            //
            // Upstream `cose-base` computes displacements from forces, then applies constraint
            // handling that *updates those displacements* (rather than hard-projecting node
            // positions after the move). Hard projection tends to over-separate constrained nodes
            // and can noticeably inflate root viewBox/max-width in parity-root mode.
            let max_d = cooling_factor * max_node_displacement;
            disps.fill((0.0, 0.0));
            // Port of `CoSELayout.moveNodes()`:
            // - displacements are calculated in `getAllNodes()` order
            // - compound displacements are propagated to (leaf) descendants before those leaves
            //   compute/clamp their own displacement
            for &idx in &all_nodes_in_layout_order {
                let Some(n) = self.nodes.get(idx) else {
                    continue;
                };

                let denom = n.no_of_children.max(1.0);
                let dx = cooling_factor * (n.spring_fx + n.repulsion_fx + n.gravitation_fx) / denom;
                let dy = cooling_factor * (n.spring_fy + n.repulsion_fy + n.gravitation_fy) / denom;

                if let Some(slot) = disps.get_mut(idx) {
                    slot.0 += dx;
                    slot.1 += dy;

                    if slot.0.abs() > max_d {
                        slot.0 = max_d * imath_sign(slot.0);
                    }
                    if slot.1.abs() > max_d {
                        slot.1 = max_d * imath_sign(slot.1);
                    }
                }

                let is_non_empty_compound = n.is_compound
                    && self
                        .children_by_owner
                        .get(idx)
                        .is_some_and(|v| !v.is_empty());
                if !is_non_empty_compound {
                    continue;
                }

                let (pdx, pdy) = disps.get(idx).copied().unwrap_or((0.0, 0.0));
                if pdx == 0.0 && pdy == 0.0 {
                    continue;
                }
                for &leaf in self.descendant_leaves.get(idx).into_iter().flatten() {
                    if leaf >= disps.len() {
                        continue;
                    }
                    disps[leaf].0 += pdx;
                    disps[leaf].1 += pdy;
                }
            }

            let disps_before_constraints =
                (debug_disps && total_iterations == 1).then(|| disps.clone());

            if let Some(rt) = constraint_rt.as_mut() {
                rt.update_displacements(
                    &self.nodes,
                    constraints,
                    &mut disps,
                    total_iterations,
                    max_d,
                    rng,
                );
            } else {
                apply_constraints_to_displacements(&self.nodes, constraints, &mut disps, max_d);
            }

            if debug_disps && total_iterations == 1 {
                if let Some(before) = disps_before_constraints {
                    eprintln!("[manatee-fcose-disp] iter=1 phase=before_constraints");
                    for (idx, (dx, dy)) in before.iter().copied().enumerate() {
                        let id = self
                            .nodes
                            .get(idx)
                            .map(|n| n.id.as_str())
                            .unwrap_or("<oob>");
                        eprintln!(
                            "[manatee-fcose-disp] idx={} id={} disp=({:.6},{:.6})",
                            idx, id, dx, dy
                        );
                    }
                }
                eprintln!("[manatee-fcose-disp] iter=1 phase=after_constraints");
                for (idx, (dx, dy)) in disps.iter().copied().enumerate() {
                    let id = self
                        .nodes
                        .get(idx)
                        .map(|n| n.id.as_str())
                        .unwrap_or("<oob>");
                    eprintln!(
                        "[manatee-fcose-disp] idx={} id={} disp=({:.6},{:.6})",
                        idx, id, dx, dy
                    );
                }
            }

            for (idx, n) in self.nodes.iter_mut().enumerate() {
                let (mdx, mdy) = disps.get(idx).copied().unwrap_or((0.0, 0.0));
                let is_non_empty_compound = n.is_compound
                    && !self
                        .children_by_owner
                        .get(idx)
                        .is_some_and(|v| v.is_empty());
                if !is_non_empty_compound {
                    n.move_by(mdx, mdy);
                    total_displacement += mdx.abs() + mdy.abs();
                }

                n.spring_fx = 0.0;
                n.spring_fy = 0.0;
                n.repulsion_fx = 0.0;
                n.repulsion_fy = 0.0;
                n.gravitation_fx = 0.0;
                n.gravitation_fy = 0.0;
            }

            last_total_displacement = total_displacement;

            if debug_positions && total_iterations == 1 {
                dump_positions("iter1", &self.nodes);
            }
        }
        if let (Some(t), Some(s)) = (timings.as_deref_mut(), iterations_start) {
            t.iterations = s.elapsed();
        }

        // Ensure compound rectangles reflect the final leaf positions before callers compute
        // component bbox/centering (e.g. `aux.relocateComponent(...)` parity).
        let _ = self.update_bounds();
        dump_positions("final", &self.nodes);

        SpringStats {
            iterations: total_iterations,
            spectral_applied,
        }
    }

    fn nodes_apply_gravity_mask(&self) -> Vec<bool> {
        let owner_connected = self.owner_graph_connected_mask();
        self.nodes
            .iter()
            .map(|n| !owner_connected.get(n.owner_idx).copied().unwrap_or(true))
            .collect()
    }

    fn owner_graph_connected_mask(&self) -> Vec<bool> {
        let owner_count = self.nodes.len() + 1;
        let mut edges_by_node: Vec<Vec<usize>> = vec![Vec::new(); self.nodes.len()];
        for (eidx, e) in self.edges.iter().enumerate() {
            if e.a < self.nodes.len() {
                edges_by_node[e.a].push(eidx);
            }
            if e.b < self.nodes.len() {
                edges_by_node[e.b].push(eidx);
            }
        }

        let mut connected: Vec<bool> = vec![true; owner_count];
        for owner in 0..owner_count {
            let nodes_in_graph = self
                .children_by_owner
                .get(owner)
                .map(|v| v.as_slice())
                .unwrap_or(&[]);
            if nodes_in_graph.is_empty() {
                connected[owner] = true;
                continue;
            }

            let mut visited: Vec<bool> = vec![false; self.nodes.len()];
            let mut queue: std::collections::VecDeque<usize> = std::collections::VecDeque::new();

            let push_with_children =
                |start: usize,
                 visited: &mut [bool],
                 queue: &mut std::collections::VecDeque<usize>| {
                    let mut stack: Vec<usize> = vec![start];
                    while let Some(cur) = stack.pop() {
                        if cur >= visited.len() {
                            continue;
                        }
                        if visited[cur] {
                            continue;
                        }
                        visited[cur] = true;
                        queue.push_back(cur);
                        if let Some(ch) = self.children_by_owner.get(cur) {
                            for &kid in ch {
                                stack.push(kid);
                            }
                        }
                    }
                };

            push_with_children(nodes_in_graph[0], &mut visited, &mut queue);

            while let Some(cur) = queue.pop_front() {
                if cur >= self.nodes.len() {
                    continue;
                }
                for &eidx in edges_by_node.get(cur).map(|v| v.as_slice()).unwrap_or(&[]) {
                    let e = &self.edges[eidx];
                    let other = if e.a == cur {
                        e.b
                    } else if e.b == cur {
                        e.a
                    } else {
                        continue;
                    };
                    let Some(mapped) = self.map_node_to_owner_graph(other, owner) else {
                        continue;
                    };
                    if !visited.get(mapped).copied().unwrap_or(false) {
                        push_with_children(mapped, &mut visited, &mut queue);
                    }
                }
            }

            connected[owner] = nodes_in_graph
                .iter()
                .all(|&nidx| visited.get(nidx).copied().unwrap_or(false));
        }

        connected
    }

    fn map_node_to_owner_graph(
        &self,
        mut node_idx: usize,
        owner_graph_idx: usize,
    ) -> Option<usize> {
        let root_owner_idx = self.root_owner_idx;
        loop {
            if node_idx >= self.nodes.len() {
                return None;
            }
            if self.nodes[node_idx].owner_idx == owner_graph_idx {
                return Some(node_idx);
            }
            let owner = self.nodes[node_idx].owner_idx;
            if owner == root_owner_idx {
                break;
            }
            node_idx = owner;
        }
        None
    }

    fn reset_edge_ideal_lengths(&mut self) {
        for e in &mut self.edges {
            e.ideal_length = e.base_ideal_length;
        }
    }

    fn adjust_intergraph_ideal_edge_lengths(&mut self) {
        if self.edges.is_empty() || self.nodes.is_empty() {
            return;
        }

        let nodes: &[SimNode] = &self.nodes;
        let inclusion_depth: &[usize] = &self.inclusion_depth;
        let root_owner_idx = self.root_owner_idx;

        for e in &mut self.edges {
            // Cache LCA-lifted endpoints for spring forces.
            let lca_owner = lca_owner_idx(nodes, root_owner_idx, e.a, e.b);
            let src_in_lca = node_in_lca_idx(nodes, root_owner_idx, e.a, lca_owner);
            let tgt_in_lca = node_in_lca_idx(nodes, root_owner_idx, e.b, lca_owner);
            e.a_in_lca = src_in_lca;
            e.b_in_lca = tgt_in_lca;

            if nodes[e.a].owner_idx == nodes[e.b].owner_idx {
                continue;
            }

            let original = e.base_ideal_length.max(1.0);

            let lca_depth = if lca_owner == root_owner_idx {
                1usize
            } else {
                inclusion_depth.get(lca_owner).copied().unwrap_or(1).max(1)
            };

            // layout-base `DEFAULT_USE_SMART_IDEAL_EDGE_LENGTH_CALCULATION = true`.
            let size_src = nodes
                .get(src_in_lca)
                .map(|n| n.estimated_size)
                .unwrap_or(Self::SIMPLE_NODE_SIZE);
            let size_tgt = nodes
                .get(tgt_in_lca)
                .map(|n| n.estimated_size)
                .unwrap_or(Self::SIMPLE_NODE_SIZE);
            e.ideal_length += size_src + size_tgt - 2.0 * Self::SIMPLE_NODE_SIZE;

            let src_depth = inclusion_depth.get(e.a).copied().unwrap_or(1).max(1);
            let tgt_depth = inclusion_depth.get(e.b).copied().unwrap_or(1).max(1);
            let hops = (src_depth + tgt_depth).saturating_sub(2 * lca_depth);
            e.ideal_length += original * Self::PER_LEVEL_IDEAL_EDGE_LENGTH_FACTOR * (hops as f64);

            if !e.ideal_length.is_finite() || e.ideal_length <= 0.0 {
                e.ideal_length = 1.0;
            }
        }
    }

    fn estimated_size(&self) -> f64 {
        // layout-base `LGraph.calcEstimatedSize()` for a flat graph:
        // - each node estimated size is (w + h) / 2
        // - graph estimated size is sum / sqrt(n)
        let n = self.nodes.len() as f64;
        if n <= 0.0 {
            return 0.0;
        }
        let sum: f64 = self.nodes.iter().map(|n| (n.width + n.height) / 2.0).sum();
        (sum / n.sqrt()).max(1.0)
    }

    fn collapse_start_positions(&mut self, scale: f64, rng: &mut XorShift64Star) {
        if self.nodes.len() <= 2 {
            return;
        }
        // Keep starts close to the origin (we relocate later).
        let jitter = (0.01 * scale).max(0.01);
        for n in self.nodes.iter_mut() {
            let jx = rng.next_f64_signed() * jitter;
            let jy = rng.next_f64_signed() * jitter;
            n.left = jx;
            n.top = jy;
        }
    }
}

fn handle_constraints_pre_layout(nodes: &mut [SimNode], c: &Constraints) {
    if nodes.is_empty() {
        return;
    }

    let mut x: Vec<f64> = nodes.iter().map(|n| n.center_x()).collect();
    let mut y: Vec<f64> = nodes.iter().map(|n| n.center_y()).collect();

    // Match `cose-base` ConstraintHandler: rotate/reflect the draft layout using an orthogonal
    // Procrustes transform derived from alignment constraints, then vote-based reflection for
    // relative placement directionality.
    if !c.align_vertical.is_empty() || !c.align_horizontal.is_empty() {
        if let Some(t) = procrustes_transform_for_alignments(&x, &y, c) {
            let tt = t.transpose();
            for i in 0..x.len() {
                let v = na::Vector2::new(x[i], y[i]);
                let r = tt * v;
                x[i] = r.x;
                y[i] = r.y;
            }
            if !c.relative.is_empty() {
                apply_reflection_for_relative_placement(&mut x, &mut y, &c.relative);
            }
        }
    } else if !c.relative.is_empty() {
        // `ConstraintHandler` also applies a relative-only transform when there are no alignment
        // constraints: it finds the largest weakly-connected component in the relative-placement
        // DAG and uses it to derive a Procrustes rotation (plus a reflection vote).
        //
        // This has an outsized effect on overall orientation and thus the parity-root viewport.
        handle_relative_only_transform(&mut x, &mut y, &c.relative);
    }

    // Enforce alignment constraints in position space.
    for group in &c.align_vertical {
        if group.len() <= 1 {
            continue;
        }
        let mut sum = 0.0;
        for &idx in group {
            sum += x[idx];
        }
        let target = sum / (group.len() as f64);
        for &idx in group {
            x[idx] = target;
        }
    }
    for group in &c.align_horizontal {
        if group.len() <= 1 {
            continue;
        }
        let mut sum = 0.0;
        for &idx in group {
            sum += y[idx];
        }
        let target = sum / (group.len() as f64);
        for &idx in group {
            y[idx] = target;
        }
    }

    // Enforce relative placement constraints in position space.
    if !c.relative.is_empty() {
        enforce_relative_placement(&mut x, &mut y, c);
    }

    for (i, n) in nodes.iter_mut().enumerate() {
        n.left = x[i] - n.width / 2.0;
        n.top = y[i] - n.height / 2.0;
    }
}

fn handle_relative_only_transform(x: &mut [f64], y: &mut [f64], rel: &[RelConstraint]) {
    use std::collections::VecDeque;

    #[derive(Debug, Clone, Copy)]
    struct Edge {
        id: usize,
        gap: f64,
    }

    let n_total = x.len().min(y.len());
    if n_total == 0 {
        return;
    }

    let mut undirected: Vec<Vec<usize>> = vec![Vec::new(); n_total];
    let mut present: Vec<bool> = vec![false; n_total];
    for r in rel {
        let (a, b) = if let (Some(left), Some(right)) = (r.left, r.right) {
            (left, right)
        } else if let (Some(top), Some(bottom)) = (r.top, r.bottom) {
            (top, bottom)
        } else {
            continue;
        };
        if a >= n_total || b >= n_total {
            continue;
        }
        undirected[a].push(b);
        undirected[b].push(a);
        present[a] = true;
        present[b] = true;
    }

    let present_count = present.iter().filter(|&&v| v).count();
    if present_count == 0 {
        return;
    }

    fn find_components(g: &[Vec<usize>], present: &[bool], node_count: usize) -> Vec<Vec<usize>> {
        let mut visited: Vec<bool> = vec![false; node_count];
        let mut out: Vec<Vec<usize>> = Vec::new();
        for start in 0..node_count {
            if !present[start] || visited[start] {
                continue;
            }

            let mut q: VecDeque<usize> = VecDeque::new();
            let mut comp: Vec<usize> = Vec::new();
            visited[start] = true;
            q.push_back(start);
            while let Some(cur) = q.pop_front() {
                comp.push(cur);
                for &n in &g[cur] {
                    if n >= node_count {
                        continue;
                    }
                    if !visited[n] {
                        visited[n] = true;
                        q.push_back(n);
                    }
                }
            }
            out.push(comp);
        }
        out
    }

    fn find_appropriate_positions(
        nodes_sorted: &[usize],
        in_comp: &[bool],
        graph: &[Vec<Edge>],
        axis: Axis,
        x: &[f64],
        y: &[f64],
    ) -> Vec<f64> {
        let node_count = x.len().min(y.len());
        let mut indeg: Vec<usize> = vec![0; node_count];
        for &src in nodes_sorted {
            if src >= node_count {
                continue;
            }
            for e in &graph[src] {
                if e.id >= node_count || !in_comp[e.id] {
                    continue;
                }
                indeg[e.id] = indeg[e.id].saturating_add(1);
            }
        }

        let mut pos: Vec<f64> = vec![f64::NEG_INFINITY; node_count];
        let mut q: VecDeque<usize> = VecDeque::new();
        for &node in nodes_sorted {
            if node >= node_count {
                continue;
            }
            if indeg[node] == 0 {
                q.push_back(node);
                pos[node] = match axis {
                    Axis::Horizontal => x[node],
                    Axis::Vertical => y[node],
                };
            }
        }

        while let Some(cur) = q.pop_front() {
            let cur_pos = pos.get(cur).copied().unwrap_or(f64::NEG_INFINITY);
            for e in graph.get(cur).into_iter().flatten() {
                if e.id >= node_count || !in_comp[e.id] {
                    continue;
                }
                let next_pos = cur_pos + e.gap;
                if pos[e.id] < next_pos {
                    pos[e.id] = next_pos;
                }
                if let Some(v) = indeg.get_mut(e.id) {
                    *v = v.saturating_sub(1);
                    if *v == 0 {
                        q.push_back(e.id);
                    }
                }
            }
        }

        pos
    }

    let components = find_components(&undirected, &present, n_total);
    if components.is_empty() {
        return;
    }

    let mut largest_idx = 0usize;
    let mut largest_sz = 0usize;
    for (i, c) in components.iter().enumerate() {
        if c.len() > largest_sz {
            largest_sz = c.len();
            largest_idx = i;
        }
    }

    if largest_sz * 2 < present_count {
        apply_reflection_for_relative_placement(x, y, rel);
        return;
    }

    let largest = &components[largest_idx];
    let mut in_comp: Vec<bool> = vec![false; n_total];
    for &idx in largest {
        if idx < n_total {
            in_comp[idx] = true;
        }
    }

    let mut nodes_sorted: Vec<usize> = largest.clone();
    nodes_sorted.sort_unstable();

    // Apply reflection votes based only on edges inside the dominant component (upstream behavior).
    let mut in_comp_constraints: Vec<RelConstraint> = Vec::new();
    let mut dag_h: Vec<Vec<Edge>> = vec![Vec::new(); n_total];
    let mut dag_v: Vec<Vec<Edge>> = vec![Vec::new(); n_total];
    for r in rel {
        if let (Some(left), Some(right)) = (r.left, r.right) {
            if left < n_total && right < n_total && in_comp[left] && in_comp[right] {
                dag_h[left].push(Edge {
                    id: right,
                    gap: r.gap,
                });
                in_comp_constraints.push(*r);
            }
        } else if let (Some(top), Some(bottom)) = (r.top, r.bottom) {
            if top < n_total && bottom < n_total && in_comp[top] && in_comp[bottom] {
                dag_v[top].push(Edge {
                    id: bottom,
                    gap: r.gap,
                });
                in_comp_constraints.push(*r);
            }
        }
    }
    apply_reflection_for_relative_placement(x, y, &in_comp_constraints);

    // Build axis DAGs and compute an "appropriate" coordinate per node using a topological
    // relaxation similar to `findAppropriatePositionForRelativePlacement`.
    let pos_h = find_appropriate_positions(&nodes_sorted, &in_comp, &dag_h, Axis::Horizontal, x, y);
    let pos_v = find_appropriate_positions(&nodes_sorted, &in_comp, &dag_v, Axis::Vertical, x, y);

    let mut source: Vec<na::Vector2<f64>> = Vec::with_capacity(largest.len());
    let mut target: Vec<na::Vector2<f64>> = Vec::with_capacity(largest.len());
    for &idx in largest {
        if idx >= n_total {
            continue;
        }
        source.push(na::Vector2::new(x[idx], y[idx]));
        let tx = pos_h.get(idx).copied().unwrap_or(x[idx]);
        let ty = pos_v.get(idx).copied().unwrap_or(y[idx]);
        target.push(na::Vector2::new(tx, ty));
    }

    if let Some(t) = procrustes_transform_from_pairs(&source, &target) {
        let tt = t.transpose();
        for i in 0..x.len().min(y.len()) {
            let v = na::Vector2::new(x[i], y[i]);
            let r = tt * v;
            x[i] = r.x;
            y[i] = r.y;
        }
    }
}

fn procrustes_transform_for_alignments(
    x: &[f64],
    y: &[f64],
    c: &Constraints,
) -> Option<na::Matrix2<f64>> {
    let mut source: Vec<na::Vector2<f64>> = Vec::new();
    let mut target: Vec<na::Vector2<f64>> = Vec::new();

    for group in &c.align_vertical {
        if group.is_empty() {
            continue;
        }
        let mut sum_x = 0.0;
        for &idx in group {
            sum_x += x[idx];
        }
        let x_pos = sum_x / (group.len() as f64);
        for &idx in group {
            source.push(na::Vector2::new(x[idx], y[idx]));
            target.push(na::Vector2::new(x_pos, y[idx]));
        }
    }

    for group in &c.align_horizontal {
        if group.is_empty() {
            continue;
        }
        let mut sum_y = 0.0;
        for &idx in group {
            sum_y += y[idx];
        }
        let y_pos = sum_y / (group.len() as f64);
        for &idx in group {
            source.push(na::Vector2::new(x[idx], y[idx]));
            target.push(na::Vector2::new(x[idx], y_pos));
        }
    }

    if source.len() <= 1 || target.len() != source.len() {
        return None;
    }

    procrustes_transform_from_pairs(&source, &target)
}

fn procrustes_transform_from_pairs(
    source: &[na::Vector2<f64>],
    target: &[na::Vector2<f64>],
) -> Option<na::Matrix2<f64>> {
    if source.len() <= 1 || target.len() != source.len() {
        return None;
    }

    let mut mean_s = na::Vector2::new(0.0, 0.0);
    let mut mean_t = na::Vector2::new(0.0, 0.0);
    for (s, t) in source.iter().zip(target.iter()) {
        mean_s += s;
        mean_t += t;
    }
    let inv_n = 1.0 / (source.len() as f64);
    mean_s *= inv_n;
    mean_t *= inv_n;

    // `ConstraintHandler` forms `tempMatrix = A'B` where A is target, B is source (mean-centered).
    let mut m = na::Matrix2::zeros();
    for (s, t) in source.iter().zip(target.iter()) {
        let sc = s - mean_s;
        let tc = t - mean_t;
        m += tc * sc.transpose();
    }

    if !(m[(0, 0)].is_finite()
        && m[(0, 1)].is_finite()
        && m[(1, 0)].is_finite()
        && m[(1, 1)].is_finite())
    {
        return None;
    }

    // Mirror layout-base `ConstraintHandler`:
    //
    // - `tempMatrix = A'B` where A is target, B is source (mean-centered)
    // - `SVD(tempMatrix) = U S V'` (JamaJS-derived routine in layout-base)
    // - `transformationMatrix = V U'`
    //
    // Use the same JamaJS-derived SVD port we already depend on for spectral layout, to avoid
    // subtle numeric drift that can break parity on symmetric constraint sets.
    let m_in = vec![vec![m[(0, 0)], m[(0, 1)]], vec![m[(1, 0)], m[(1, 1)]]];
    let svd = spectral::svd_jama(&m_in)?;
    if svd.u.len() < 2 || svd.v.len() < 2 {
        return None;
    }
    let u = &svd.u;
    let v = &svd.v;

    // T = V * U^T
    let t00 = v[0][0] * u[0][0] + v[0][1] * u[0][1];
    let t01 = v[0][0] * u[1][0] + v[0][1] * u[1][1];
    let t10 = v[1][0] * u[0][0] + v[1][1] * u[0][1];
    let t11 = v[1][0] * u[1][0] + v[1][1] * u[1][1];
    Some(na::Matrix2::new(t00, t01, t10, t11))
}

fn apply_reflection_for_relative_placement(x: &mut [f64], y: &mut [f64], rel: &[RelConstraint]) {
    let mut reflect_on_y = 0;
    let mut not_reflect_on_y = 0;
    let mut reflect_on_x = 0;
    let mut not_reflect_on_x = 0;

    for r in rel {
        if let (Some(left), Some(right)) = (r.left, r.right) {
            if x[left] - x[right] >= 0.0 {
                reflect_on_y += 1;
            } else {
                not_reflect_on_y += 1;
            }
        } else if let (Some(top), Some(bottom)) = (r.top, r.bottom) {
            if y[top] - y[bottom] >= 0.0 {
                reflect_on_x += 1;
            } else {
                not_reflect_on_x += 1;
            }
        }
    }

    if reflect_on_y > not_reflect_on_y && reflect_on_x > not_reflect_on_x {
        for i in 0..x.len() {
            x[i] = -x[i];
            y[i] = -y[i];
        }
    } else if reflect_on_y > not_reflect_on_y {
        for v in x.iter_mut() {
            *v = -*v;
        }
    } else if reflect_on_x > not_reflect_on_x {
        for v in y.iter_mut() {
            *v = -*v;
        }
    }
}

fn enforce_relative_placement(x: &mut [f64], y: &mut [f64], c: &Constraints) {
    #[derive(Debug, Clone, Copy)]
    struct Neighbor {
        id: usize,
        gap: f64,
    }

    let n = x.len().min(y.len());
    if n == 0 {
        return;
    }

    fn enforce_relative_placement_no_align_small(
        x: &mut [f64],
        y: &mut [f64],
        rel: &[RelConstraint],
        n: usize,
    ) {
        use std::collections::VecDeque;

        fn build_axis_dag_keys(
            axis: Axis,
            rel: &[RelConstraint],
            n: usize,
        ) -> (Vec<usize>, Vec<Vec<Neighbor>>) {
            let mut keys: Vec<usize> = Vec::new();
            let mut seen: Vec<bool> = vec![false; n];
            let mut dag: Vec<Vec<Neighbor>> = vec![Vec::new(); n];

            for r in rel {
                match axis {
                    Axis::Horizontal => {
                        let (Some(left), Some(right)) = (r.left, r.right) else {
                            continue;
                        };
                        if left >= n || right >= n {
                            continue;
                        }
                        if !seen[left] {
                            seen[left] = true;
                            keys.push(left);
                        }
                        if !seen[right] {
                            seen[right] = true;
                            keys.push(right);
                        }
                        dag[left].push(Neighbor {
                            id: right,
                            gap: r.gap,
                        });
                    }
                    Axis::Vertical => {
                        let (Some(top), Some(bottom)) = (r.top, r.bottom) else {
                            continue;
                        };
                        if top >= n || bottom >= n {
                            continue;
                        }
                        if !seen[top] {
                            seen[top] = true;
                            keys.push(top);
                        }
                        if !seen[bottom] {
                            seen[bottom] = true;
                            keys.push(bottom);
                        }
                        dag[top].push(Neighbor {
                            id: bottom,
                            gap: r.gap,
                        });
                    }
                }
            }

            (keys, dag)
        }

        fn build_rev(keys: &[usize], dag: &[Vec<Neighbor>], n: usize) -> Vec<Vec<Neighbor>> {
            let mut rev: Vec<Vec<Neighbor>> = vec![Vec::new(); n];
            for &src in keys {
                if src >= n {
                    continue;
                }
                for e in &dag[src] {
                    if e.id >= n {
                        continue;
                    }
                    rev[e.id].push(Neighbor {
                        id: src,
                        gap: e.gap,
                    });
                }
            }
            rev
        }

        fn pos_before(key: usize, axis: Axis, x: &[f64], y: &[f64]) -> f64 {
            match axis {
                Axis::Horizontal => x[key],
                Axis::Vertical => y[key],
            }
        }

        fn component_sources(
            keys: &[usize],
            dag: &[Vec<Neighbor>],
            rev: &[Vec<Neighbor>],
            n: usize,
        ) -> Vec<Vec<usize>> {
            let mut undirected: Vec<Vec<usize>> = vec![Vec::new(); n];
            for &src in keys {
                if src >= n {
                    continue;
                }
                for e in &dag[src] {
                    if e.id >= n {
                        continue;
                    }
                    undirected[src].push(e.id);
                    undirected[e.id].push(src);
                }
            }

            let mut visited: Vec<bool> = vec![false; n];
            let mut out: Vec<Vec<usize>> = Vec::new();
            for &start in keys {
                if start >= n || visited[start] {
                    continue;
                }
                let mut q: VecDeque<usize> = VecDeque::new();
                let mut comp: Vec<usize> = Vec::new();
                visited[start] = true;
                q.push_back(start);
                while let Some(cur) = q.pop_front() {
                    comp.push(cur);
                    for &next in &undirected[cur] {
                        if next < n && !visited[next] {
                            visited[next] = true;
                            q.push_back(next);
                        }
                    }
                }

                let mut sources: Vec<usize> = Vec::new();
                for &node in &comp {
                    if node < n && rev[node].is_empty() {
                        sources.push(node);
                    }
                }
                out.push(sources);
            }
            out
        }

        fn find_appropriate_positions(
            keys: &[usize],
            dag: &[Vec<Neighbor>],
            axis: Axis,
            n: usize,
            x: &[f64],
            y: &[f64],
            sources: &[Vec<usize>],
        ) -> Vec<f64> {
            let mut in_deg: Vec<usize> = vec![0; n];
            for &src in keys {
                for e in &dag[src] {
                    in_deg[e.id] = in_deg[e.id].saturating_add(1);
                }
            }

            let mut position: Vec<f64> = vec![0.0; n];
            let mut past_bits: Vec<u64> = vec![0; n];
            let mut past_order: Vec<Vec<usize>> = vec![Vec::new(); n];
            let mut q: VecDeque<usize> = VecDeque::new();

            for &k in keys {
                position[k] = f64::NEG_INFINITY;
                if in_deg[k] == 0 {
                    q.push_back(k);
                }
                past_bits[k] = 1u64 << (k as u64);
                past_order[k] = vec![k];
            }

            for component in sources {
                if component.is_empty() {
                    continue;
                }
                let mut sum = 0.0;
                for &node in component {
                    sum += pos_before(node, axis, x, y);
                }
                let avg = sum / (component.len() as f64);
                for &node in component {
                    position[node] = avg;
                }
            }

            while let Some(cur) = q.pop_front() {
                let cur_pos = position[cur];
                for neigh in &dag[cur] {
                    let want = cur_pos + neigh.gap;
                    if position[neigh.id] < want {
                        position[neigh.id] = want;
                    }
                    in_deg[neigh.id] = in_deg[neigh.id].saturating_sub(1);
                    if in_deg[neigh.id] == 0 {
                        q.push_back(neigh.id);
                    }

                    let mut merged_bits = past_bits[cur];
                    let mut merged_order: Vec<usize> = past_order[cur].clone();
                    for &v in &past_order[neigh.id] {
                        let bit = 1u64 << (v as u64);
                        if (merged_bits & bit) == 0 {
                            merged_bits |= bit;
                            merged_order.push(v);
                        }
                    }
                    past_bits[neigh.id] = merged_bits;
                    past_order[neigh.id] = merged_order;
                }
            }

            let mut sink_nodes: Vec<usize> = Vec::new();
            for &k in keys {
                if dag[k].is_empty() {
                    sink_nodes.push(k);
                }
            }

            let mut comp_bits: Vec<u64> = Vec::new();
            let mut comp_order: Vec<Vec<usize>> = Vec::new();
            for &k in keys {
                if !sink_nodes.contains(&k) || past_order[k].is_empty() {
                    continue;
                }
                let first = past_order[k][0];
                let first_bit = 1u64 << (first as u64);
                if let Some(idx) = comp_bits.iter().position(|b| (*b & first_bit) != 0) {
                    let mut bits = comp_bits[idx];
                    let mut order = comp_order[idx].clone();
                    for &v in &past_order[k] {
                        let bit = 1u64 << (v as u64);
                        if (bits & bit) == 0 {
                            bits |= bit;
                            order.push(v);
                        }
                    }
                    comp_bits[idx] = bits;
                    comp_order[idx] = order;
                } else {
                    comp_bits.push(past_bits[k]);
                    comp_order.push(past_order[k].clone());
                }
            }

            for comp in comp_order {
                let mut min_before = f64::INFINITY;
                let mut max_before = f64::NEG_INFINITY;
                let mut min_after = f64::INFINITY;
                let mut max_after = f64::NEG_INFINITY;
                for &node in &comp {
                    let before = pos_before(node, axis, x, y);
                    let after = position[node];
                    min_before = min_before.min(before);
                    max_before = max_before.max(before);
                    min_after = min_after.min(after);
                    max_after = max_after.max(after);
                }
                let diff = ((min_before + max_before) / 2.0) - ((min_after + max_after) / 2.0);
                for &node in &comp {
                    position[node] += diff;
                }
            }

            position
        }

        let (keys_h, dag_h) = build_axis_dag_keys(Axis::Horizontal, rel, n);
        if !keys_h.is_empty() {
            let rev_h = build_rev(&keys_h, &dag_h, n);
            let sources = component_sources(&keys_h, &dag_h, &rev_h, n);
            let pos =
                find_appropriate_positions(&keys_h, &dag_h, Axis::Horizontal, n, x, y, &sources);
            for &k in &keys_h {
                x[k] = pos[k];
            }
        }

        let (keys_v, dag_v) = build_axis_dag_keys(Axis::Vertical, rel, n);
        if !keys_v.is_empty() {
            let rev_v = build_rev(&keys_v, &dag_v, n);
            let sources = component_sources(&keys_v, &dag_v, &rev_v, n);
            let pos =
                find_appropriate_positions(&keys_v, &dag_v, Axis::Vertical, n, x, y, &sources);
            for &k in &keys_v {
                y[k] = pos[k];
            }
        }
    }

    if c.align_vertical.is_empty() && c.align_horizontal.is_empty() && n <= 64 {
        enforce_relative_placement_no_align_small(x, y, &c.relative, n);
        return;
    }

    // Dummy mappings for alignment constraints (per-axis, matching `ConstraintHandler`).
    let mut dummy_to_nodes_for_vertical_alignment: Vec<Vec<usize>> = Vec::new();
    let mut node_to_dummy_for_vertical_alignment: Vec<Option<usize>> = vec![None; n];
    for (i, group) in c.align_vertical.iter().enumerate() {
        let dummy = n + i;
        dummy_to_nodes_for_vertical_alignment.push(group.clone());
        for &idx in group {
            if idx < n {
                node_to_dummy_for_vertical_alignment[idx] = Some(dummy);
            }
        }
    }
    let mut dummy_pos_for_vertical_alignment: Vec<f64> = dummy_to_nodes_for_vertical_alignment
        .iter()
        .map(|g| x[*g.first().unwrap_or(&0)])
        .collect();

    let mut dummy_to_nodes_for_horizontal_alignment: Vec<Vec<usize>> = Vec::new();
    let mut node_to_dummy_for_horizontal_alignment: Vec<Option<usize>> = vec![None; n];
    for (i, group) in c.align_horizontal.iter().enumerate() {
        let dummy = n + i;
        dummy_to_nodes_for_horizontal_alignment.push(group.clone());
        for &idx in group {
            if idx < n {
                node_to_dummy_for_horizontal_alignment[idx] = Some(dummy);
            }
        }
    }
    let mut dummy_pos_for_horizontal_alignment: Vec<f64> = dummy_to_nodes_for_horizontal_alignment
        .iter()
        .map(|g| y[*g.first().unwrap_or(&0)])
        .collect();

    fn ensure_key(map: &mut IndexMap<usize, Vec<Neighbor>>, key: usize) {
        if !map.contains_key(&key) {
            map.insert(key, Vec::new());
        }
    }

    let mut dag_h: IndexMap<usize, Vec<Neighbor>> = IndexMap::new();
    let mut dag_v: IndexMap<usize, Vec<Neighbor>> = IndexMap::new();
    for r in &c.relative {
        if let (Some(left), Some(right)) = (r.left, r.right) {
            let src = node_to_dummy_for_vertical_alignment[left].unwrap_or(left);
            let dst = node_to_dummy_for_vertical_alignment[right].unwrap_or(right);
            ensure_key(&mut dag_h, src);
            ensure_key(&mut dag_h, dst);
            dag_h.get_mut(&src).unwrap().push(Neighbor {
                id: dst,
                gap: r.gap,
            });
        } else if let (Some(top), Some(bottom)) = (r.top, r.bottom) {
            let src = node_to_dummy_for_horizontal_alignment[top].unwrap_or(top);
            let dst = node_to_dummy_for_horizontal_alignment[bottom].unwrap_or(bottom);
            ensure_key(&mut dag_v, src);
            ensure_key(&mut dag_v, dst);
            dag_v.get_mut(&src).unwrap().push(Neighbor {
                id: dst,
                gap: r.gap,
            });
        }
    }

    fn dag_to_undirected(dag: &IndexMap<usize, Vec<Neighbor>>) -> IndexMap<usize, Vec<Neighbor>> {
        let mut u: IndexMap<usize, Vec<Neighbor>> = IndexMap::new();
        for (&k, _) in dag.iter() {
            u.insert(k, Vec::new());
        }
        for (&k, neigh) in dag.iter() {
            for n in neigh {
                u.get_mut(&k).unwrap().push(*n);
                u.get_mut(&n.id)
                    .unwrap()
                    .push(Neighbor { id: k, gap: n.gap });
            }
        }
        u
    }

    fn dag_to_reversed(dag: &IndexMap<usize, Vec<Neighbor>>) -> IndexMap<usize, Vec<Neighbor>> {
        let mut r: IndexMap<usize, Vec<Neighbor>> = IndexMap::new();
        for (&k, _) in dag.iter() {
            r.insert(k, Vec::new());
        }
        for (&k, neigh) in dag.iter() {
            for n in neigh {
                r.get_mut(&n.id)
                    .unwrap()
                    .push(Neighbor { id: k, gap: n.gap });
            }
        }
        r
    }

    fn find_components(undirected: &IndexMap<usize, Vec<Neighbor>>) -> Vec<Vec<usize>> {
        use std::collections::{HashSet, VecDeque};
        let mut visited: HashSet<usize> = HashSet::new();
        let mut out: Vec<Vec<usize>> = Vec::new();
        for (&k, _) in undirected.iter() {
            if visited.contains(&k) {
                continue;
            }
            let mut q: VecDeque<usize> = VecDeque::new();
            let mut comp: Vec<usize> = Vec::new();
            q.push_back(k);
            visited.insert(k);
            while let Some(cur) = q.pop_front() {
                comp.push(cur);
                for n in &undirected[&cur] {
                    if visited.insert(n.id) {
                        q.push_back(n.id);
                    }
                }
            }
            out.push(comp);
        }
        out
    }

    fn component_sources(
        dag: &IndexMap<usize, Vec<Neighbor>>,
        rev: &IndexMap<usize, Vec<Neighbor>>,
    ) -> Vec<Vec<usize>> {
        let undirected = dag_to_undirected(dag);
        let comps = find_components(&undirected);
        let mut out: Vec<Vec<usize>> = Vec::new();
        for comp in comps {
            let mut sources: Vec<usize> = Vec::new();
            for node in comp {
                if rev.get(&node).is_none_or(|v| v.is_empty()) {
                    sources.push(node);
                }
            }
            out.push(sources);
        }
        out
    }

    fn pos_before(key: usize, axis: Axis, n: usize, x: &[f64], y: &[f64], dummy: &[f64]) -> f64 {
        if key < n {
            match axis {
                Axis::Horizontal => x[key],
                Axis::Vertical => y[key],
            }
        } else {
            dummy[key - n]
        }
    }

    fn find_appropriate_positions(
        dag: &IndexMap<usize, Vec<Neighbor>>,
        axis: Axis,
        n: usize,
        x: &[f64],
        y: &[f64],
        dummy_pos: &[f64],
        component_sources: &[Vec<usize>],
    ) -> IndexMap<usize, f64> {
        use std::collections::VecDeque;

        let mut in_deg: IndexMap<usize, usize> = IndexMap::new();
        for (&k, _) in dag.iter() {
            in_deg.insert(k, 0);
        }
        for (&_k, neigh) in dag.iter() {
            for n2 in neigh {
                *in_deg.get_mut(&n2.id).unwrap() += 1;
            }
        }

        let mut position: IndexMap<usize, f64> = IndexMap::new();
        let mut past: IndexMap<usize, IndexSet<usize>> = IndexMap::new();
        let mut q: VecDeque<usize> = VecDeque::new();

        for (&k, &deg) in in_deg.iter() {
            position.insert(k, f64::NEG_INFINITY);
            if deg == 0 {
                q.push_back(k);
            }
            past.insert(k, IndexSet::from([k]));
        }

        // Align sources of each component (enforcement path, empty fixed-node set).
        for component in component_sources {
            if component.is_empty() {
                continue;
            }
            let mut sum = 0.0;
            for &node in component {
                sum += pos_before(node, axis, n, x, y, dummy_pos);
            }
            let avg = sum / (component.len() as f64);
            for &node in component {
                position.insert(node, avg);
            }
        }

        while let Some(cur) = q.pop_front() {
            let cur_pos = position[&cur];
            for neigh in &dag[&cur] {
                let want = cur_pos + neigh.gap;
                if position[&neigh.id] < want {
                    position.insert(neigh.id, want);
                }
                let deg = in_deg.get_mut(&neigh.id).unwrap();
                *deg = deg.saturating_sub(1);
                if *deg == 0 {
                    q.push_back(neigh.id);
                }
                let mut merged: IndexSet<usize> = past[&cur].clone();
                for v in past[&neigh.id].iter().copied() {
                    merged.insert(v);
                }
                past.insert(neigh.id, merged);
            }
        }

        // Readjust position after enforcement.
        let mut sink_nodes: IndexSet<usize> = IndexSet::new();
        for (&k, neigh) in dag.iter() {
            if neigh.is_empty() {
                sink_nodes.insert(k);
            }
        }

        let mut components: Vec<IndexSet<usize>> = Vec::new();
        for (&k, set) in past.iter() {
            if !sink_nodes.contains(&k) || set.is_empty() {
                continue;
            }
            let first = *set.iter().next().unwrap();
            if let Some(idx) = components.iter().position(|c| c.contains(&first)) {
                let mut merged = components[idx].clone();
                for v in set.iter().copied() {
                    merged.insert(v);
                }
                components[idx] = merged;
            } else {
                components.push(set.clone());
            }
        }

        for comp in components {
            let mut min_before = f64::INFINITY;
            let mut max_before = f64::NEG_INFINITY;
            let mut min_after = f64::INFINITY;
            let mut max_after = f64::NEG_INFINITY;
            for &node in comp.iter() {
                let before = pos_before(node, axis, n, x, y, dummy_pos);
                let after = position[&node];
                min_before = min_before.min(before);
                max_before = max_before.max(before);
                min_after = min_after.min(after);
                max_after = max_after.max(after);
            }
            let diff = ((min_before + max_before) / 2.0) - ((min_after + max_after) / 2.0);
            for &node in comp.iter() {
                position.insert(node, position[&node] + diff);
            }
        }

        position
    }

    if !dag_h.is_empty() {
        let rev = dag_to_reversed(&dag_h);
        let sources = component_sources(&dag_h, &rev);
        let pos = find_appropriate_positions(
            &dag_h,
            Axis::Horizontal,
            n,
            x,
            y,
            &dummy_pos_for_vertical_alignment,
            &sources,
        );
        for (&key, &v) in pos.iter() {
            if key < n {
                x[key] = v;
            } else {
                let di = key - n;
                for &idx in &dummy_to_nodes_for_vertical_alignment[di] {
                    x[idx] = v;
                }
                dummy_pos_for_vertical_alignment[di] = v;
            }
        }
    }

    if !dag_v.is_empty() {
        let rev = dag_to_reversed(&dag_v);
        let sources = component_sources(&dag_v, &rev);
        let pos = find_appropriate_positions(
            &dag_v,
            Axis::Vertical,
            n,
            x,
            y,
            &dummy_pos_for_horizontal_alignment,
            &sources,
        );
        for (&key, &v) in pos.iter() {
            if key < n {
                y[key] = v;
            } else {
                let di = key - n;
                for &idx in &dummy_to_nodes_for_horizontal_alignment[di] {
                    y[idx] = v;
                }
                dummy_pos_for_horizontal_alignment[di] = v;
            }
        }
    }
}

fn apply_constraints_to_displacements(
    nodes: &[SimNode],
    c: &Constraints,
    disps: &mut [(f64, f64)],
    max_d: f64,
) {
    // Alignments: enforce exact alignment by adjusting displacements to a shared target line.
    for group in &c.align_horizontal {
        if group.len() <= 1 {
            continue;
        }
        let mut sum = 0.0;
        let mut cnt = 0.0;
        for &idx in group {
            sum += nodes[idx].center_y() + disps[idx].1;
            cnt += 1.0;
        }
        if cnt > 0.0 {
            let target = sum / cnt;
            for &idx in group {
                disps[idx].1 += target - (nodes[idx].center_y() + disps[idx].1);
            }
        }
    }
    for group in &c.align_vertical {
        if group.len() <= 1 {
            continue;
        }
        let mut sum = 0.0;
        let mut cnt = 0.0;
        for &idx in group {
            sum += nodes[idx].center_x() + disps[idx].0;
            cnt += 1.0;
        }
        if cnt > 0.0 {
            let target = sum / cnt;
            for &idx in group {
                disps[idx].0 += target - (nodes[idx].center_x() + disps[idx].0);
            }
        }
    }

    // Relative placements: iteratively relax displacements to satisfy minimum center gaps.
    // This is a small, deterministic approximation of `cose-base` constraint handling.
    for _ in 0..4 {
        let mut changed = false;
        for r in &c.relative {
            if let (Some(left), Some(right)) = (r.left, r.right) {
                let new_gap = (nodes[right].center_x() + disps[right].0)
                    - (nodes[left].center_x() + disps[left].0);
                if new_gap < r.gap {
                    let delta = r.gap - new_gap;
                    disps[left].0 -= delta / 2.0;
                    disps[right].0 += delta / 2.0;
                    changed = true;
                }
            }
            if let (Some(top), Some(bottom)) = (r.top, r.bottom) {
                let new_gap = (nodes[bottom].center_y() + disps[bottom].1)
                    - (nodes[top].center_y() + disps[top].1);
                if new_gap < r.gap {
                    let delta = r.gap - new_gap;
                    disps[top].1 -= delta / 2.0;
                    disps[bottom].1 += delta / 2.0;
                    changed = true;
                }
            }
        }
        if !changed {
            break;
        }
    }

    // Re-apply per-axis displacement caps (matching the upstream `calculateDisplacement` clamp).
    if max_d.is_finite() && max_d > 0.0 {
        for (dx, dy) in disps {
            if dx.abs() > max_d {
                *dx = max_d * imath_sign(*dx);
            }
            if dy.abs() > max_d {
                *dy = max_d * imath_sign(*dy);
            }
        }
    }
}

#[derive(Debug, Clone)]
struct XorShift64Star {
    state: u64,
    calls: u64,
}

impl XorShift64Star {
    fn new(seed: u64) -> Self {
        Self {
            state: seed.max(1),
            calls: 0,
        }
    }

    #[allow(dead_code)]
    fn mix_u64(&mut self, v: u64) {
        // One-way mix to decorrelate node indices.
        self.state ^= v.wrapping_mul(0x9E3779B97F4A7C15_u64);
        let _ = self.next_u64();
    }

    fn next_u64(&mut self) -> u64 {
        self.calls = self.calls.wrapping_add(1);
        let mut x = self.state;
        x ^= x >> 12;
        x ^= x << 25;
        x ^= x >> 27;
        self.state = x;
        x.wrapping_mul(0x2545F4914F6CDD1D_u64)
    }

    fn next_f64_signed(&mut self) -> f64 {
        // Map to [-1, 1] (exclusive).
        let u = self.next_u64() >> 11;
        let v = (u as f64) / ((1u64 << 53) as f64);
        (v * 2.0) - 1.0
    }

    fn next_f64_unit(&mut self) -> f64 {
        // Map to [0, 1) with 53 bits of precision.
        let u = self.next_u64() >> 11;
        (u as f64) / ((1u64 << 53) as f64)
    }

    fn calls(&self) -> u64 {
        self.calls
    }

    fn next_usize(&mut self, upper: usize) -> usize {
        if upper <= 1 {
            return 0;
        }
        // Match the seeded upstream baselines which override `Math.random()` with a 53-bit float
        // derived from `nextU64() >> 11`, then select indices via
        // `Math.floor(Math.random() * upper)`.
        //
        // Using `% upper` introduces modulo bias and (more importantly for parity) can yield a
        // different first sample pivot for small graphs (e.g. upper=3), which cascades into a
        // different spectral embedding orientation.
        let v = self.next_f64_unit();
        let idx = (v * (upper as f64)).floor() as usize;
        idx.min(upper - 1)
    }
}

#[cfg(test)]
mod tests {
    use super::{
        BoundsExtras, Constraints, RelConstraint, RepulsionGrid, SimNode, XorShift64Star,
        apply_reflection_for_relative_placement, procrustes_transform_for_alignments,
    };
    use nalgebra as na;

    fn node_at(left: f64, top: f64, w: f64, h: f64) -> SimNode {
        SimNode {
            id: "n".to_string(),
            parent: None,
            owner_idx: 0,
            is_compound: false,
            width: w,
            height: h,
            bounds_extras: BoundsExtras::default(),
            estimated_size: (w + h) / 2.0,
            left,
            top,
            spring_fx: 0.0,
            spring_fy: 0.0,
            repulsion_fx: 0.0,
            repulsion_fy: 0.0,
            gravitation_fx: 0.0,
            gravitation_fy: 0.0,
            no_of_children: 1.0,
            padding: 0.0,
            surrounding: Vec::new(),
            grid_start_x: 0,
            grid_finish_x: 0,
            grid_start_y: 0,
            grid_finish_y: 0,
        }
    }

    #[test]
    fn xorshift64star_next_f64_unit_matches_seeded_upstream_baseline() {
        // Mirrors the JS prelude in `xtask` used to generate deterministic upstream SVGs:
        //
        // - xorshift64* (same shift/multiply constants)
        // - `Math.random = () => Number(nextU64() >> 11n) / 2^53`
        let mut rng = XorShift64Star::new(1);
        let expected = [
            0.28083505005035947,
            0.6711372530266764,
            0.7258461452833668,
            0.303529299965799,
            0.056176763098259475,
        ];
        for (i, &e) in expected.iter().enumerate() {
            let v = rng.next_f64_unit();
            assert!(
                (v - e).abs() < 1e-15,
                "unexpected rng value at {i}: got {v}, expected {e}"
            );
        }
    }

    #[test]
    fn xorshift64star_next_usize_matches_js_floor_random_times_upper() {
        // For seed=1, the first `Math.random()` value is ~0.2808 so `floor(r * 3) == 0`.
        // Using `% 3` on the underlying u64 yields `1`, which would diverge from the upstream
        // spectral sampling path for small graphs.
        let mut rng = XorShift64Star::new(1);
        assert_eq!(rng.next_usize(3), 0);
    }

    #[test]
    fn repulsion_grid_surrounding_excludes_processed_nodes() {
        // Build a tiny 1D-ish layout:
        //
        // - node0 and node1 are exactly within range
        // - node2 is far outside range
        let repulsion_range = 10.0;
        let mut nodes = vec![
            node_at(0.0, 0.0, 10.0, 10.0),
            node_at(20.0, 0.0, 10.0, 10.0),
            node_at(200.0, 0.0, 10.0, 10.0),
        ];
        let mut left = f64::INFINITY;
        let mut top = f64::INFINITY;
        let mut right = f64::NEG_INFINITY;
        let mut bottom = f64::NEG_INFINITY;
        for n in &nodes {
            left = left.min(n.left);
            top = top.min(n.top);
            right = right.max(n.left + n.width);
            bottom = bottom.max(n.top + n.height);
        }
        let node_order = [0usize, 1, 2];
        let grid = RepulsionGrid::build(
            left,
            top,
            right,
            bottom,
            &mut nodes,
            repulsion_range,
            &node_order,
        )
        .expect("grid");

        let mut processed = vec![false; nodes.len()];
        grid.refresh_node_surrounding(0, &mut nodes, &processed, repulsion_range);
        assert_eq!(nodes[0].surrounding, vec![1]);

        processed[0] = true;
        grid.refresh_node_surrounding(1, &mut nodes, &processed, repulsion_range);
        assert!(
            !nodes[1].surrounding.contains(&0),
            "node1 should not include already-processed node0"
        );
    }

    #[test]
    fn relative_placement_gap_is_center_to_center() {
        use super::{Constraints, RelConstraint, apply_constraints_to_displacements};

        let nodes = vec![
            node_at(0.0, 0.0, 10.0, 10.0),  // center_x = 5
            node_at(20.0, 0.0, 10.0, 10.0), // center_x = 25
        ];
        let mut disps = vec![(0.0, 0.0); nodes.len()];

        let c = Constraints {
            align_horizontal: Vec::new(),
            align_vertical: Vec::new(),
            relative: vec![RelConstraint {
                left: Some(0),
                right: Some(1),
                top: None,
                bottom: None,
                gap: 50.0,
            }],
        };

        apply_constraints_to_displacements(&nodes, &c, &mut disps, 1e9);
        let gap = (nodes[1].center_x() + disps[1].0) - (nodes[0].center_x() + disps[0].0);
        assert!((gap - 50.0).abs() < 1e-9, "gap: got {gap}");
    }

    #[test]
    fn rect_clip_points_matches_layout_base_igeometry_getintersection2() {
        // Expected values computed via layout-base@2.0.1:
        //
        // `IGeometry.getIntersection(rectA, rectB, out)` where:
        // - rectA = (-274.090946,-129.901919,80,80)
        // - rectB = (512.630977,-782.722296,80,80)
        let a = node_at(-274.090_946, -129.901_919, 80.0, 80.0);
        let b = node_at(512.630_977, -782.722_296, 80.0, 80.0);
        let (ax, ay, bx, by) = super::rect_clip_points(&a, &b);

        let eps = 1e-6;
        assert!((ax - -194.090_946).abs() < eps, "ax: got {ax}");
        assert!((ay - -123.093_844_020_246_31).abs() < eps, "ay: got {ay}");
        assert!((bx - 512.630_977).abs() < eps, "bx: got {bx}");
        assert!((by - -709.530_370_979_753_7).abs() < eps, "by: got {by}");
    }

    #[test]
    fn constraint_procrustes_transform_matches_upstream_fixture_025_checkpoint() {
        // Ground truth extracted via `tools/debug/arch_probe_fcose_vs_upstream_025.js`:
        //
        // - `draft.debug.recomputed`: raw spectral coordinates (pre-relocation)
        //
        // Upstream `ConstraintHandler` applies a Procrustes + reflection transform directly to the
        // raw coordinates; component relocation (`aux.relocateComponent(componentCenter, ...)`) is
        // performed later by cytoscape-fcose and shows up in `draft.pos` / `fromSpectral.*`.
        //
        // This test intentionally isolates the transform-only step on `draft.debug.recomputed`.
        //
        // This test guards against subtle transpose/sign mistakes in our Procrustes port.
        let ids = ["a", "b", "c", "d", "e", "f"];
        let draft = [
            (-69.77618192016361, 79.87553327881355),
            (34.28258770643722, 100.36650015929253),
            (104.06591551872783, 20.494458759991097),
            (69.78035458033064, -79.87753496744543),
            (-34.28895079233283, -100.37063982916823),
            (-104.06372509299923, -20.48831740148356),
        ];
        let expected = [
            (-63.516289670902054, 84.938197098671),
            (41.796999300167016, 97.47687430142939),
            (105.32067914788601, 12.54161697884377),
            (63.52029847371475, -84.94050953250945),
            (-41.80365806342419, -97.48051938039694),
            (-105.31802918744155, -12.535659466037792),
        ];

        let mut nodes: Vec<SimNode> = Vec::new();
        for (i, (x, y)) in draft.iter().copied().enumerate() {
            nodes.push(SimNode {
                id: ids[i].to_string(),
                parent: None,
                owner_idx: i,
                is_compound: false,
                width: 80.0,
                height: 80.0,
                bounds_extras: BoundsExtras::default(),
                estimated_size: 80.0,
                left: x - 40.0,
                top: y - 40.0,
                spring_fx: 0.0,
                spring_fy: 0.0,
                repulsion_fx: 0.0,
                repulsion_fy: 0.0,
                gravitation_fx: 0.0,
                gravitation_fy: 0.0,
                no_of_children: 1.0,
                padding: 0.0,
                surrounding: Vec::new(),
                grid_start_x: 0,
                grid_finish_x: 0,
                grid_start_y: 0,
                grid_finish_y: 0,
            });
        }

        let c = Constraints {
            align_horizontal: vec![vec![0, 5], vec![2, 3]],
            align_vertical: vec![vec![1, 2], vec![3, 4]],
            relative: vec![
                RelConstraint {
                    left: Some(0),
                    right: Some(5),
                    top: None,
                    bottom: None,
                    gap: 120.0,
                },
                RelConstraint {
                    left: Some(4),
                    right: Some(1),
                    top: None,
                    bottom: None,
                    gap: 120.0,
                },
                RelConstraint {
                    left: None,
                    right: None,
                    top: Some(1),
                    bottom: Some(2),
                    gap: 120.0,
                },
                RelConstraint {
                    left: None,
                    right: None,
                    top: Some(4),
                    bottom: Some(3),
                    gap: 120.0,
                },
                RelConstraint {
                    left: Some(3),
                    right: Some(2),
                    top: None,
                    bottom: None,
                    gap: 120.0,
                },
            ],
        };

        let mut x: Vec<f64> = nodes.iter().map(|n| n.center_x()).collect();
        let mut y: Vec<f64> = nodes.iter().map(|n| n.center_y()).collect();

        let t = procrustes_transform_for_alignments(&x, &y, &c).expect("transform");
        let tt = t.transpose();
        for i in 0..x.len() {
            let v = na::Vector2::new(x[i], y[i]);
            let r = tt * v;
            x[i] = r.x;
            y[i] = r.y;
        }
        apply_reflection_for_relative_placement(&mut x, &mut y, &c.relative);

        for i in 0..ids.len() {
            let (ex, ey) = expected[i];
            let dx = (x[i] - ex).abs();
            let dy = (y[i] - ey).abs();
            assert!(
                dx < 1e-9 && dy < 1e-9,
                "mismatch for {}: got=({:.12},{:.12}) expected=({:.12},{:.12}) d=({:.3e},{:.3e})",
                ids[i],
                x[i],
                y[i],
                ex,
                ey,
                dx,
                dy
            );
        }
    }
}

fn rects_intersect(a: &SimNode, b: &SimNode) -> bool {
    // Mirror layout-base `RectangleD.intersects`: touching edges count as intersection.
    !(a.right() < b.left || a.bottom() < b.top || b.right() < a.left || b.bottom() < a.top)
}

fn get_cardinal_direction(slope: f64, slope_prime: f64, line: i32) -> i32 {
    if slope > slope_prime {
        line
    } else {
        1 + (line % 4)
    }
}

fn rect_clip_points(a: &SimNode, b: &SimNode) -> (f64, f64, f64, f64) {
    // Port of layout-base `IGeometry.getIntersection2(rectA, rectB, result)`.
    //
    // result[0-1] contains clip point on rectA; result[2-3] contains clip point on rectB.
    let p1x = a.center_x();
    let p1y = a.center_y();
    let p2x = b.center_x();
    let p2y = b.center_y();

    if rects_intersect(a, b) {
        return (p1x, p1y, p2x, p2y);
    }

    let top_left_ax = a.left;
    let top_left_ay = a.top;
    let top_right_ax = a.right();
    let bottom_left_ax = a.left;
    let bottom_left_ay = a.bottom();
    let bottom_right_ax = a.right();
    let half_width_a = a.half_w();
    let half_height_a = a.half_h();

    let top_left_bx = b.left;
    let top_left_by = b.top;
    let top_right_bx = b.right();
    let bottom_left_bx = b.left;
    let bottom_left_by = b.bottom();
    let bottom_right_bx = b.right();
    let half_width_b = b.half_w();
    let half_height_b = b.half_h();

    let mut clip_ax = p1x;
    let mut clip_ay = p1y;
    let mut clip_bx = p2x;
    let mut clip_by = p2y;

    if p1x == p2x {
        if p1y > p2y {
            return (p1x, top_left_ay, p2x, bottom_left_by);
        } else if p1y < p2y {
            return (p1x, bottom_left_ay, p2x, top_left_by);
        }
    } else if p1y == p2y {
        if p1x > p2x {
            return (top_left_ax, p1y, top_right_bx, p2y);
        } else if p1x < p2x {
            return (top_right_ax, p1y, top_left_bx, p2y);
        }
    } else {
        let slope_a = a.height / a.width;
        let slope_b = b.height / b.width;
        let slope_prime = (p2y - p1y) / (p2x - p1x);

        let mut clip_a_found = false;
        let mut clip_b_found = false;

        if -slope_a == slope_prime {
            if p1x > p2x {
                clip_ax = bottom_left_ax;
                clip_ay = bottom_left_ay;
                clip_a_found = true;
            } else {
                clip_ax = top_right_ax;
                clip_ay = top_left_ay;
                clip_a_found = true;
            }
        } else if slope_a == slope_prime {
            if p1x > p2x {
                clip_ax = top_left_ax;
                clip_ay = top_left_ay;
                clip_a_found = true;
            } else {
                clip_ax = bottom_right_ax;
                clip_ay = bottom_left_ay;
                clip_a_found = true;
            }
        }

        if -slope_b == slope_prime {
            if p2x > p1x {
                clip_bx = bottom_left_bx;
                clip_by = bottom_left_by;
                clip_b_found = true;
            } else {
                clip_bx = top_right_bx;
                clip_by = top_left_by;
                clip_b_found = true;
            }
        } else if slope_b == slope_prime {
            if p2x > p1x {
                clip_bx = top_left_bx;
                clip_by = top_left_by;
                clip_b_found = true;
            } else {
                clip_bx = bottom_right_bx;
                clip_by = bottom_left_by;
                clip_b_found = true;
            }
        }

        if !clip_a_found || !clip_b_found {
            let (card_a, card_b) = if p1x > p2x {
                if p1y > p2y {
                    (
                        get_cardinal_direction(slope_a, slope_prime, 4),
                        get_cardinal_direction(slope_b, slope_prime, 2),
                    )
                } else {
                    (
                        get_cardinal_direction(-slope_a, slope_prime, 3),
                        get_cardinal_direction(-slope_b, slope_prime, 1),
                    )
                }
            } else if p1y > p2y {
                (
                    get_cardinal_direction(-slope_a, slope_prime, 1),
                    get_cardinal_direction(-slope_b, slope_prime, 3),
                )
            } else {
                (
                    get_cardinal_direction(slope_a, slope_prime, 2),
                    get_cardinal_direction(slope_b, slope_prime, 4),
                )
            };

            if !clip_a_found {
                match card_a {
                    1 => {
                        clip_ay = top_left_ay;
                        clip_ax = p1x + -half_height_a / slope_prime;
                    }
                    2 => {
                        clip_ax = bottom_right_ax;
                        clip_ay = p1y + half_width_a * slope_prime;
                    }
                    3 => {
                        clip_ay = bottom_left_ay;
                        clip_ax = p1x + half_height_a / slope_prime;
                    }
                    4 => {
                        clip_ax = bottom_left_ax;
                        clip_ay = p1y + -half_width_a * slope_prime;
                    }
                    _ => {}
                }
            }

            if !clip_b_found {
                match card_b {
                    1 => {
                        clip_by = top_left_by;
                        clip_bx = p2x + -half_height_b / slope_prime;
                    }
                    2 => {
                        clip_bx = bottom_right_bx;
                        clip_by = p2y + half_width_b * slope_prime;
                    }
                    3 => {
                        clip_by = bottom_left_by;
                        clip_bx = p2x + half_height_b / slope_prime;
                    }
                    4 => {
                        clip_bx = bottom_left_bx;
                        clip_by = p2y + -half_width_b * slope_prime;
                    }
                    _ => {}
                }
            }
        }
    }

    (clip_ax, clip_ay, clip_bx, clip_by)
}

fn calc_repulsion_force(
    a: &SimNode,
    b: &SimNode,
    min_repulsion_dist: f64,
    separation_buffer: f64,
) -> (f64, f64) {
    if rects_intersect(a, b) {
        let (ox, oy) = calc_separation_amount(a, b, separation_buffer);
        let repulsion_fx = 2.0 * ox;
        let repulsion_fy = 2.0 * oy;

        // layout-base: scale overlap separation by a children constant so large compounds move
        // more slowly than leaves (and to reduce oscillation).
        let denom = (a.no_of_children + b.no_of_children).max(1.0);
        let children_constant = (a.no_of_children * b.no_of_children) / denom;

        // Return a force delta to be applied as:
        // - nodeA += rfx/rfy
        // - nodeB -= rfx/rfy
        (
            -children_constant * repulsion_fx,
            -children_constant * repulsion_fy,
        )
    } else {
        let (ax, ay, bx, by) = rect_clip_points(a, b);
        let mut dx = bx - ax;
        let mut dy = by - ay;

        if dx.abs() < min_repulsion_dist {
            dx = imath_sign(dx) * min_repulsion_dist;
        }
        if dy.abs() < min_repulsion_dist {
            dy = imath_sign(dy) * min_repulsion_dist;
        }

        let dist_sq = dx * dx + dy * dy;
        let dist = dist_sq.sqrt();
        if dist_sq == 0.0 || dist == 0.0 {
            return (0.0, 0.0);
        }

        // layout-base:
        // `(nodeA.nodeRepulsion/2 + nodeB.nodeRepulsion/2) * noOfChildrenA * noOfChildrenB / dist^2`.
        // Default node repulsion is 4500 for both nodes.
        let repulsion_force =
            SimGraph::DEFAULT_REPULSION_STRENGTH * a.no_of_children * b.no_of_children / dist_sq;
        let repulsion_fx = repulsion_force * dx / dist;
        let repulsion_fy = repulsion_force * dy / dist;
        (-repulsion_fx, -repulsion_fy)
    }
}

#[derive(Debug, Clone)]
struct RepulsionGrid {
    left: f64,
    top: f64,
    size_x: i32,
    size_y: i32,
    // Flat grid: cells[x * size_y + y] contains node indices.
    cells: Vec<Vec<usize>>,
}

impl RepulsionGrid {
    fn idx(&self, x: i32, y: i32) -> usize {
        (x as usize) * (self.size_y as usize) + (y as usize)
    }

    fn cell(&self, x: i32, y: i32) -> &[usize] {
        &self.cells[self.idx(x, y)]
    }

    fn build(
        left: f64,
        top: f64,
        right: f64,
        bottom: f64,
        nodes: &mut [SimNode],
        repulsion_range: f64,
        node_order: &[usize],
    ) -> Option<Self> {
        if nodes.is_empty() {
            return None;
        }
        if !repulsion_range.is_finite() || repulsion_range <= 0.0 {
            return None;
        }
        if !(left.is_finite() && top.is_finite() && right.is_finite() && bottom.is_finite()) {
            return None;
        }

        let w = (right - left).max(1.0);
        let h = (bottom - top).max(1.0);
        if !(w.is_finite() && h.is_finite()) {
            return None;
        }

        // layout-base `FDLayout.calcGrid`: size = ceil((graph.right - graph.left) / repulsionRange).
        let size_x = ((w / repulsion_range).ceil() as i32).max(1);
        let size_y = ((h / repulsion_range).ceil() as i32).max(1);
        let mut cells: Vec<Vec<usize>> = vec![Vec::new(); (size_x as usize) * (size_y as usize)];

        // Mirror layout-base `addNodeToGrid`: push the node into every cell that intersects the
        // node's rect, using top-left anchored coordinates.
        //
        // Important: layout-base inserts nodes into the grid in `getAllNodes()` order (see
        // `FDLayout.updateGrid()`), which is observable because the surrounding list is built
        // by iterating over the grid cells and preserving insertion order. Matching this order
        // reduces floating-point accumulation drift in parity tests.
        for &idx in node_order {
            let Some(n) = nodes.get_mut(idx) else {
                continue;
            };
            let start_x = ((n.left - left) / repulsion_range).floor() as i32;
            let finish_x = ((n.right() - left) / repulsion_range).floor() as i32;
            let start_y = ((n.top - top) / repulsion_range).floor() as i32;
            let finish_y = ((n.bottom() - top) / repulsion_range).floor() as i32;

            n.grid_start_x = start_x;
            n.grid_finish_x = finish_x;
            n.grid_start_y = start_y;
            n.grid_finish_y = finish_y;

            for gx in start_x..=finish_x {
                if gx < 0 || gx >= size_x {
                    continue;
                }
                for gy in start_y..=finish_y {
                    if gy < 0 || gy >= size_y {
                        continue;
                    }
                    let cell_idx = (gx as usize) * (size_y as usize) + (gy as usize);
                    cells[cell_idx].push(idx);
                }
            }
        }

        Some(Self {
            left,
            top,
            size_x,
            size_y,
            cells,
        })
    }

    fn refresh_node_surrounding(
        &self,
        node_idx: usize,
        nodes: &mut [SimNode],
        processed: &[bool],
        repulsion_range: f64,
    ) {
        let start_x = nodes[node_idx].grid_start_x;
        let finish_x = nodes[node_idx].grid_finish_x;
        let start_y = nodes[node_idx].grid_start_y;
        let finish_y = nodes[node_idx].grid_finish_y;

        let mut seen: Vec<bool> = vec![false; nodes.len()];
        let mut surrounding: Vec<usize> = Vec::new();

        for gx in (start_x - 1)..=(finish_x + 1) {
            if gx < 0 || gx >= self.size_x {
                continue;
            }
            for gy in (start_y - 1)..=(finish_y + 1) {
                if gy < 0 || gy >= self.size_y {
                    continue;
                }
                for &other in self.cell(gx, gy) {
                    if other == node_idx {
                        continue;
                    }
                    if processed.get(other).copied().unwrap_or(false) {
                        continue;
                    }
                    if seen[other] {
                        continue;
                    }
                    if nodes[node_idx].owner_idx != nodes[other].owner_idx {
                        continue;
                    }

                    let dx = (nodes[node_idx].center_x() - nodes[other].center_x()).abs()
                        - (nodes[node_idx].half_w() + nodes[other].half_w());
                    let dy = (nodes[node_idx].center_y() - nodes[other].center_y()).abs()
                        - (nodes[node_idx].half_h() + nodes[other].half_h());
                    if dx <= repulsion_range && dy <= repulsion_range {
                        seen[other] = true;
                        surrounding.push(other);
                    }
                }
            }
        }

        nodes[node_idx].surrounding = surrounding;
    }
}

fn calc_separation_amount(a: &SimNode, b: &SimNode, separation_buffer: f64) -> (f64, f64) {
    debug_assert!(rects_intersect(a, b));

    let (dir_x, dir_y) = decide_directions_for_overlapping_nodes(a, b);

    // Port of layout-base `IGeometry.calcSeparationAmount` overlap logic used by FDLayout.
    let mut overlap_x = a.right().min(b.right()) - a.left.max(b.left);
    let mut overlap_y = a.bottom().min(b.bottom()) - a.top.max(b.top);

    if (a.left <= b.left) && (a.right() >= b.right()) {
        overlap_x += (b.left - a.left).min(a.right() - b.right());
    } else if (b.left <= a.left) && (b.right() >= a.right()) {
        overlap_x += (a.left - b.left).min(b.right() - a.right());
    }
    if (a.top <= b.top) && (a.bottom() >= b.bottom()) {
        overlap_y += (b.top - a.top).min(a.bottom() - b.bottom());
    } else if (b.top <= a.top) && (b.bottom() >= a.bottom()) {
        overlap_y += (a.top - b.top).min(b.bottom() - a.bottom());
    }

    let mut slope = ((b.center_y() - a.center_y()) / (b.center_x() - a.center_x())).abs();
    if (b.center_y() == a.center_y()) && (b.center_x() == a.center_x()) {
        slope = 1.0;
    }

    let mut move_by_y = slope * overlap_x;
    let mut move_by_x = overlap_y / slope;
    if overlap_x < move_by_x {
        move_by_x = overlap_x;
    } else {
        move_by_y = overlap_y;
    }

    let dx = -dir_x * ((move_by_x / 2.0) + separation_buffer);
    let dy = -dir_y * ((move_by_y / 2.0) + separation_buffer);
    (dx, dy)
}

fn decide_directions_for_overlapping_nodes(a: &SimNode, b: &SimNode) -> (f64, f64) {
    let dir_x = if a.center_x() < b.center_x() {
        -1.0
    } else {
        1.0
    };
    let dir_y = if a.center_y() < b.center_y() {
        -1.0
    } else {
        1.0
    };
    (dir_x, dir_y)
}