From 591a5d9fe76fa2539c33948734b4672b9d722480 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Thu, 7 May 2026 03:15:29 +0100 Subject: [PATCH] Add Barnes-Hut 3D repulsion --- README.md | 15 +-- RESEARCH.md | 8 +- src/lib.rs | 303 +++++++++++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 314 insertions(+), 12 deletions(-) diff --git a/README.md b/README.md index 07c5b34..64b9b34 100644 --- a/README.md +++ b/README.md @@ -98,13 +98,15 @@ bureaucratic head injury. `Layout3dConfig::repulsion_mode` controls the expensive part: -- `RepulsionMode::SpatialGrid` is the default realtime path. +- `RepulsionMode::BarnesHut` is the default realtime path. +- `RepulsionMode::SpatialGrid` keeps a simpler local approximation available. - `RepulsionMode::Exact` keeps all-pairs repulsion for small graphs or quality comparisons. -The grid path is approximate local repulsion, not Barnes-Hut. It is the first -practical step toward interactive scale; large graphs still want multilevel -coarsening and a better far-field approximation. +Barnes-Hut uses an octree and `barnes_hut_theta` as the opening-angle knob. +Lower theta values are more accurate; higher values are faster. Large graphs +still want multilevel coarsening, because even good approximation does not make +every visual decision worth rendering individually. ## Graph Analysis And Folding @@ -164,6 +166,5 @@ The first release favors a small deterministic core over a sprawling layout cathedral. The public API exposes generated constraints so callers can inspect, reuse, tune, or replace the force relaxation phase later. -The current repulsion pass is exact pairwise repulsion. That is simple and -stable for modest graphs; larger graphs should move toward a Barnes-Hut or grid -approximation before pretending heroism is a scalability plan. +The realtime solver defaults to Barnes-Hut repulsion. Exact pairwise repulsion +is still available for small graphs and quality comparisons. diff --git a/RESEARCH.md b/RESEARCH.md index 3f44a03..c9d57c1 100644 --- a/RESEARCH.md +++ b/RESEARCH.md @@ -16,7 +16,8 @@ These features are implemented in the crate now: - Persistent `Layout3dSolver` for realtime frame-by-frame iteration. - Warm-start positions for graph edits and interaction. -- Spatial-grid approximate repulsion for the 3D solver. +- Barnes-Hut octree repulsion for the 3D solver. +- Spatial-grid approximate repulsion as a simpler fallback. - SCC detection for cyclic bodies. - weak components for broad spatial islands. - articulation points and bridges for hinge-like connectors. @@ -72,13 +73,14 @@ crate is adaptive step control and Barnes-Hut for large interactive layouts. Implementation direction: - Replace fixed cooling with adaptive speed. -- Upgrade spatial-grid repulsion to Barnes-Hut for far-field approximation. +- Tune Barnes-Hut and add benchmarks against exact/grid modes. - Expose live-tick APIs for Bevy interaction instead of only batch layout. Current state: - `Layout3dSolver` exists. -- `RepulsionMode::SpatialGrid` is the default realtime path. +- `RepulsionMode::BarnesHut` is the default realtime path. +- `RepulsionMode::SpatialGrid` remains available as a simpler local approximation. - `RepulsionMode::Exact` remains available for small graphs and comparisons. Source: diff --git a/src/lib.rs b/src/lib.rs index 8094b58..45c3f3c 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -153,6 +153,10 @@ pub struct Layout3dConfig { pub grid_cell_size: f32, /// Number of neighboring cells searched around each node. pub grid_radius: i32, + /// Barnes-Hut opening angle. Lower is more accurate; higher is faster. + pub barnes_hut_theta: f32, + /// Maximum Barnes-Hut octree depth. + pub barnes_hut_max_depth: usize, } impl Default for Layout3dConfig { @@ -172,9 +176,11 @@ impl Default for Layout3dConfig { community_cohesion_strength: 0.025, core_anchor_strength: 0.018, clustering_cohesion_strength: 0.02, - repulsion_mode: RepulsionMode::SpatialGrid, + repulsion_mode: RepulsionMode::BarnesHut, grid_cell_size: 180.0, grid_radius: 1, + barnes_hut_theta: 0.65, + barnes_hut_max_depth: 24, } } } @@ -183,8 +189,10 @@ impl Default for Layout3dConfig { pub enum RepulsionMode { /// Exact all-pairs repulsion. Best for quality on small graphs. Exact, - /// Approximate local repulsion through a uniform spatial grid. + /// Barnes-Hut octree approximation. Best default for realtime 3D layouts. #[default] + BarnesHut, + /// Approximate local repulsion through a uniform spatial grid. SpatialGrid, } @@ -340,6 +348,9 @@ impl Layout3dSolver { RepulsionMode::Exact => { apply_repulsion_3d_exact(&self.positions, &mut self.forces, &self.config) } + RepulsionMode::BarnesHut => { + apply_repulsion_3d_barnes_hut(&self.positions, &mut self.forces, &self.config) + } RepulsionMode::SpatialGrid => { apply_repulsion_3d_grid(&self.positions, &mut self.forces, &self.config) } @@ -1659,6 +1670,242 @@ fn apply_repulsion_3d_exact( } } +fn apply_repulsion_3d_barnes_hut( + positions: &[(f32, f32, f32)], + forces: &mut [(f32, f32, f32)], + config: &Layout3dConfig, +) { + if positions.len() < 2 { + return; + } + + let tree = BarnesHutTree::new(positions, config.barnes_hut_max_depth); + let theta = config.barnes_hut_theta.max(0.01); + for target in 0..positions.len() { + tree.apply_force(target, theta, positions, forces, config); + } +} + +#[derive(Clone, Debug)] +struct BarnesHutTree { + nodes: Vec, +} + +impl BarnesHutTree { + fn new(positions: &[(f32, f32, f32)], max_depth: usize) -> Self { + let (center, half_size) = barnes_hut_bounds(positions); + let mut tree = Self { + nodes: vec![BarnesHutNode::new(center, half_size)], + }; + + for index in 0..positions.len() { + tree.insert(0, index, positions, 0, max_depth); + } + + tree + } + + fn insert( + &mut self, + node_index: usize, + body: usize, + positions: &[(f32, f32, f32)], + depth: usize, + max_depth: usize, + ) { + self.nodes[node_index].add_mass(positions[body]); + + if self.nodes[node_index].is_leaf() + && (self.nodes[node_index].bodies.is_empty() || depth >= max_depth) + { + self.nodes[node_index].bodies.push(body); + return; + } + + if self.nodes[node_index].is_leaf() { + let existing = std::mem::take(&mut self.nodes[node_index].bodies); + self.subdivide(node_index); + for old_body in existing { + let child = self.child_for_body(node_index, positions[old_body]); + self.insert(child, old_body, positions, depth + 1, max_depth); + } + } + + let child = self.child_for_body(node_index, positions[body]); + self.insert(child, body, positions, depth + 1, max_depth); + } + + fn subdivide(&mut self, node_index: usize) { + let center = self.nodes[node_index].center; + let child_half = self.nodes[node_index].half_size * 0.5; + let mut children = [None; 8]; + + for (slot, child) in children.iter_mut().enumerate() { + let offset = ( + if slot & 1 == 0 { + -child_half + } else { + child_half + }, + if slot & 2 == 0 { + -child_half + } else { + child_half + }, + if slot & 4 == 0 { + -child_half + } else { + child_half + }, + ); + let child_center = ( + center.0 + offset.0, + center.1 + offset.1, + center.2 + offset.2, + ); + let child_index = self.nodes.len(); + self.nodes + .push(BarnesHutNode::new(child_center, child_half)); + *child = Some(child_index); + } + + self.nodes[node_index].children = children; + } + + fn child_for_body(&self, node_index: usize, position: (f32, f32, f32)) -> usize { + let node = &self.nodes[node_index]; + let mut slot = 0usize; + if position.0 >= node.center.0 { + slot |= 1; + } + if position.1 >= node.center.1 { + slot |= 2; + } + if position.2 >= node.center.2 { + slot |= 4; + } + node.children[slot].unwrap_or(node_index) + } + + fn apply_force( + &self, + target: usize, + theta: f32, + positions: &[(f32, f32, f32)], + forces: &mut [(f32, f32, f32)], + config: &Layout3dConfig, + ) { + self.apply_node_force(0, target, theta, positions, forces, config); + } + + fn apply_node_force( + &self, + node_index: usize, + target: usize, + theta: f32, + positions: &[(f32, f32, f32)], + forces: &mut [(f32, f32, f32)], + config: &Layout3dConfig, + ) { + let node = &self.nodes[node_index]; + if node.mass == 0 { + return; + } + + if node.is_leaf() { + for &body in &node.bodies { + if body != target { + apply_pair_repulsion_3d(target, body, positions, forces, config); + } + } + return; + } + + let target_position = positions[target]; + let dx = target_position.0 - node.center_of_mass.0; + let dy = target_position.1 - node.center_of_mass.1; + let dz = target_position.2 - node.center_of_mass.2; + let distance = (dx * dx + dy * dy + dz * dz).sqrt().max(0.01); + let width = node.half_size * 2.0; + + if !node.contains(target_position) && width / distance < theta { + apply_aggregate_repulsion_3d(target, node, positions, forces, config); + return; + } + + for child in node.children.iter().flatten() { + self.apply_node_force(*child, target, theta, positions, forces, config); + } + } +} + +#[derive(Clone, Debug)] +struct BarnesHutNode { + center: (f32, f32, f32), + half_size: f32, + mass: usize, + center_of_mass: (f32, f32, f32), + children: [Option; 8], + bodies: Vec, +} + +impl BarnesHutNode { + fn new(center: (f32, f32, f32), half_size: f32) -> Self { + Self { + center, + half_size, + mass: 0, + center_of_mass: (0.0, 0.0, 0.0), + children: [None; 8], + bodies: Vec::new(), + } + } + + fn add_mass(&mut self, position: (f32, f32, f32)) { + let old_mass = self.mass as f32; + let new_mass = old_mass + 1.0; + self.center_of_mass.0 = (self.center_of_mass.0 * old_mass + position.0) / new_mass; + self.center_of_mass.1 = (self.center_of_mass.1 * old_mass + position.1) / new_mass; + self.center_of_mass.2 = (self.center_of_mass.2 * old_mass + position.2) / new_mass; + self.mass += 1; + } + + fn is_leaf(&self) -> bool { + self.children.iter().all(Option::is_none) + } + + fn contains(&self, position: (f32, f32, f32)) -> bool { + position.0 >= self.center.0 - self.half_size + && position.0 <= self.center.0 + self.half_size + && position.1 >= self.center.1 - self.half_size + && position.1 <= self.center.1 + self.half_size + && position.2 >= self.center.2 - self.half_size + && position.2 <= self.center.2 + self.half_size + } +} + +fn barnes_hut_bounds(positions: &[(f32, f32, f32)]) -> ((f32, f32, f32), f32) { + let mut min = positions[0]; + let mut max = positions[0]; + + for &position in positions.iter().skip(1) { + min.0 = min.0.min(position.0); + min.1 = min.1.min(position.1); + min.2 = min.2.min(position.2); + max.0 = max.0.max(position.0); + max.1 = max.1.max(position.1); + max.2 = max.2.max(position.2); + } + + let center = ( + (min.0 + max.0) * 0.5, + (min.1 + max.1) * 0.5, + (min.2 + max.2) * 0.5, + ); + let span = (max.0 - min.0).max(max.1 - min.1).max(max.2 - min.2); + (center, (span * 0.5).max(1.0) + 0.001) +} + fn apply_repulsion_3d_grid( positions: &[(f32, f32, f32)], forces: &mut [(f32, f32, f32)], @@ -1705,6 +1952,24 @@ fn grid_cell(position: (f32, f32, f32), cell_size: f32) -> (i32, i32, i32) { ) } +fn apply_aggregate_repulsion_3d( + target: usize, + node: &BarnesHutNode, + positions: &[(f32, f32, f32)], + forces: &mut [(f32, f32, f32)], + config: &Layout3dConfig, +) { + let dx = positions[target].0 - node.center_of_mass.0; + let dy = positions[target].1 - node.center_of_mass.1; + let dz = positions[target].2 - node.center_of_mass.2; + let dist2 = (dx * dx + dy * dy + dz * dz).max(0.01); + let dist = dist2.sqrt(); + let magnitude = config.base.repulsion_strength * node.mass as f32 / dist2; + forces[target].0 += dx / dist * magnitude * config.horizontal_repulsion; + forces[target].1 += dy / dist * magnitude * config.vertical_repulsion; + forces[target].2 += dz / dist * magnitude * config.horizontal_repulsion; +} + fn apply_pair_repulsion_3d( a: usize, b: usize, @@ -2099,4 +2364,38 @@ mod tests { assert_eq!(solver.snapshot().nodes.len(), 2); } + + #[test] + fn barnes_hut_repulsion_mode_runs_as_default() { + let mut graph = Graph::new(); + let a = graph.add_node(1.0); + let b = graph.add_node(1.0); + let c = graph.add_node(1.0); + graph.add_edge(a, b); + graph.add_edge(b, c); + + let mut solver = Layout3dSolver::new(graph, Layout3dConfig::default()); + solver.tick(3); + let snapshot = solver.snapshot(); + + assert_eq!(snapshot.nodes.len(), 3); + assert_eq!(snapshot.constraints.ranks.len(), 3); + } + + #[test] + fn spatial_grid_repulsion_mode_still_runs() { + let mut graph = Graph::new(); + let a = graph.add_node(1.0); + let b = graph.add_node(1.0); + graph.add_edge(a, b); + let config = Layout3dConfig { + repulsion_mode: RepulsionMode::SpatialGrid, + ..Layout3dConfig::default() + }; + + let mut solver = Layout3dSolver::new(graph, config); + solver.tick(2); + + assert_eq!(solver.snapshot().nodes.len(), 2); + } }