1use std::collections::BTreeSet;
37
38#[derive(Clone, Debug, PartialEq)]
44pub struct GraphCsr {
45 pub n: usize,
46 pub offsets: Vec<u32>,
48 pub targets: Vec<u32>,
50 pub weights: Vec<f32>,
52}
53
54impl GraphCsr {
55 #[must_use]
65 pub fn from_edges(n: usize, edges: &[(u32, u32)], weights: Option<&[f32]>) -> Self {
66 let mut adj: Vec<BTreeSet<u32>> = vec![BTreeSet::new(); n];
67 for &(a, b) in edges {
68 let (a, b) = (a as usize, b as usize);
69 if a == b || a >= n || b >= n {
70 continue;
71 }
72 adj[a].insert(b as u32);
73 adj[b].insert(a as u32);
74 }
75 let mut offsets = Vec::with_capacity(n + 1);
76 let mut targets = Vec::new();
77 offsets.push(0u32);
78 for nbrs in &adj {
79 targets.extend(nbrs.iter().copied());
80 offsets.push(targets.len() as u32);
81 }
82 let weights = (0..n)
83 .map(|i| weights.and_then(|w| w.get(i)).copied().filter(|w| *w > 0.0).unwrap_or(1.0))
84 .collect();
85 Self { n, offsets, targets, weights }
86 }
87
88 #[must_use]
90 pub fn degree_sum(&self) -> usize {
91 self.targets.len()
92 }
93
94 #[must_use]
96 pub fn neighbours(&self, i: usize) -> &[u32] {
97 let (s, e) = (self.offsets[i] as usize, self.offsets[i + 1] as usize);
98 &self.targets[s..e]
99 }
100}
101
102#[cfg(feature = "arrow")]
103impl GraphCsr {
104 pub fn from_arrow(
110 n: usize,
111 src: &arrow_array::Int64Array,
112 dst: &arrow_array::Int64Array,
113 weight: Option<&arrow_array::Float64Array>,
114 ) -> Self {
115 use arrow_array::Array;
116 let mut edges = Vec::with_capacity(src.len());
117 for i in 0..src.len().min(dst.len()) {
118 if src.is_null(i) || dst.is_null(i) {
119 continue;
120 }
121 let (a, b) = (src.value(i), dst.value(i));
122 if a < 0 || b < 0 {
123 continue;
124 }
125 edges.push((a as u32, b as u32));
126 }
127 let w: Option<Vec<f32>> = weight.map(|w| {
128 (0..n).map(|i| if i < w.len() && !w.is_null(i) { w.value(i) as f32 } else { 1.0 }).collect()
129 });
130 Self::from_edges(n, &edges, w.as_deref())
131 }
132}
133
134#[derive(Clone, Copy, Debug, PartialEq)]
136pub struct LayoutParams {
137 pub repulsion: f32,
139 pub attraction: f32,
141 pub ideal_len: f32,
143 pub gravity: f32,
145 pub damping: f32,
147 pub dt: f32,
149 pub min_dist: f32,
151}
152
153impl Default for LayoutParams {
154 fn default() -> Self {
155 Self {
156 repulsion: 1.0,
157 attraction: 2.0,
158 ideal_len: 1.0,
159 gravity: 0.05,
160 damping: 0.9,
161 dt: 0.05,
162 min_dist: 0.05,
163 }
164 }
165}
166
167#[derive(Clone, Debug, PartialEq)]
170pub struct LayoutState {
171 pub pos: Vec<[f32; 2]>,
172 pub vel: Vec<[f32; 2]>,
173}
174
175impl LayoutState {
176 #[must_use]
179 pub fn seed(n: usize, scale: f32) -> Self {
180 const GOLDEN: f32 = 2.399_963_2; let pos = (0..n)
182 .map(|i| {
183 let r = scale * ((i as f32 + 0.5) / n.max(1) as f32).sqrt();
184 let a = i as f32 * GOLDEN;
185 [r * a.cos(), r * a.sin()]
186 })
187 .collect();
188 Self { pos, vel: vec![[0.0, 0.0]; n] }
189 }
190
191 #[must_use]
194 pub fn kinetic_energy(&self) -> f32 {
195 self.vel.iter().map(|v| v[0] * v[0] + v[1] * v[1]).sum()
196 }
197
198 #[must_use]
201 pub fn total_edge_len(&self, g: &GraphCsr) -> f32 {
202 let mut sum = 0.0;
203 for i in 0..g.n {
204 for &t in g.neighbours(i) {
205 if (t as usize) > i {
206 let d = [self.pos[t as usize][0] - self.pos[i][0], self.pos[t as usize][1] - self.pos[i][1]];
207 sum += (d[0] * d[0] + d[1] * d[1]).sqrt();
208 }
209 }
210 }
211 sum
212 }
213}
214
215#[must_use]
218pub fn step_cpu(s: &LayoutState, g: &GraphCsr, p: &LayoutParams) -> LayoutState {
219 let n = g.n;
220 let mut out = LayoutState { pos: Vec::with_capacity(n), vel: Vec::with_capacity(n) };
221 let min2 = p.min_dist * p.min_dist;
222 for i in 0..n {
223 let pi = s.pos[i];
224 let wi = g.weights[i];
225 let mut force = [0.0f32, 0.0f32];
226 for (j, pj) in s.pos.iter().enumerate() {
228 if j == i {
229 continue;
230 }
231 let d = [pi[0] - pj[0], pi[1] - pj[1]];
232 let dist2 = (d[0] * d[0] + d[1] * d[1]).max(min2);
233 let f = p.repulsion * wi * g.weights[j] / dist2;
234 force[0] += d[0] * f;
235 force[1] += d[1] * f;
236 }
237 for &t in g.neighbours(i) {
239 let pt = s.pos[t as usize];
240 let d = [pt[0] - pi[0], pt[1] - pi[1]];
241 let dist = (d[0] * d[0] + d[1] * d[1]).sqrt().max(1e-4);
242 let f = p.attraction * (dist - p.ideal_len) / dist;
243 force[0] += d[0] * f;
244 force[1] += d[1] * f;
245 }
246 force[0] -= pi[0] * p.gravity;
248 force[1] -= pi[1] * p.gravity;
249 let inv_m = 1.0 / wi.max(1e-4);
251 let vx = (s.vel[i][0] + force[0] * inv_m * p.dt) * p.damping;
252 let vy = (s.vel[i][1] + force[1] * inv_m * p.dt) * p.damping;
253 out.vel.push([vx, vy]);
254 out.pos.push([pi[0] + vx * p.dt, pi[1] + vy * p.dt]);
255 }
256 out
257}
258
259#[must_use]
262pub fn relax_cpu(g: &GraphCsr, p: &LayoutParams, iters: usize) -> LayoutState {
263 let mut s = LayoutState::seed(g.n, p.ideal_len * (g.n as f32).sqrt());
264 for _ in 0..iters {
265 s = step_cpu(&s, g, p);
266 }
267 s
268}
269
270#[cfg(feature = "gpu")]
272mod gpu {
273 use super::{GraphCsr, LayoutParams, LayoutState};
274 use wgpu::util::DeviceExt;
275
276 #[repr(C)]
280 #[derive(Clone, Copy, Debug, Default, PartialEq, bytemuck::Pod, bytemuck::Zeroable)]
281 pub struct GpuNode {
282 pub a: [f32; 4],
283 pub b: [f32; 4],
284 }
285
286 #[repr(C)]
287 #[derive(Clone, Copy, bytemuck::Pod, bytemuck::Zeroable)]
288 struct Uniforms {
289 a: [f32; 4], b: [f32; 4], }
292
293 pub const LAYOUT_WGSL: &str = include_str!("layout.wgsl");
295
296 pub fn pack(s: &LayoutState, weights: &[f32]) -> Vec<GpuNode> {
298 (0..s.pos.len())
299 .map(|i| GpuNode { a: [s.pos[i][0], s.pos[i][1], s.vel[i][0], s.vel[i][1]], b: [weights[i], 0.0, 0.0, 0.0] })
300 .collect()
301 }
302
303 pub struct GpuLayout {
307 pipeline: wgpu::ComputePipeline,
308 bgl: wgpu::BindGroupLayout,
309 uniforms: wgpu::Buffer,
310 offsets: wgpu::Buffer,
312 targets: wgpu::Buffer,
313 src: wgpu::Buffer,
315 dst: wgpu::Buffer,
316 count: u32,
317 }
318
319 impl GpuLayout {
320 pub fn new(device: &wgpu::Device, g: &GraphCsr, seed: &LayoutState) -> Self {
323 let shader = device.create_shader_module(wgpu::ShaderModuleDescriptor {
324 label: Some("graphview_layout"),
325 source: wgpu::ShaderSource::Wgsl(LAYOUT_WGSL.into()),
326 });
327 let storage = |ro: bool| wgpu::BindGroupLayoutEntry {
328 binding: 0,
329 visibility: wgpu::ShaderStages::COMPUTE,
330 ty: wgpu::BindingType::Buffer { ty: wgpu::BufferBindingType::Storage { read_only: ro }, has_dynamic_offset: false, min_binding_size: None },
331 count: None,
332 };
333 let bgl = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor {
334 label: Some("graphview_layout_bgl"),
335 entries: &[
336 wgpu::BindGroupLayoutEntry {
337 binding: 0,
338 visibility: wgpu::ShaderStages::COMPUTE,
339 ty: wgpu::BindingType::Buffer { ty: wgpu::BufferBindingType::Uniform, has_dynamic_offset: false, min_binding_size: None },
340 count: None,
341 },
342 wgpu::BindGroupLayoutEntry { binding: 1, ..storage(true) },
343 wgpu::BindGroupLayoutEntry { binding: 2, ..storage(false) },
344 wgpu::BindGroupLayoutEntry { binding: 3, ..storage(true) },
345 wgpu::BindGroupLayoutEntry { binding: 4, ..storage(true) },
346 ],
347 });
348 let pipeline = device.create_compute_pipeline(&wgpu::ComputePipelineDescriptor {
349 label: Some("graphview_layout_step"),
350 layout: Some(&device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor {
351 label: Some("graphview_layout_pll"),
352 bind_group_layouts: &[Some(&bgl)],
353 immediate_size: 0,
354 })),
355 module: &shader,
356 entry_point: Some("cs_layout"),
357 compilation_options: Default::default(),
358 cache: None,
359 });
360
361 let nodes = pack(seed, &g.weights);
362 let node_usage = wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC | wgpu::BufferUsages::COPY_DST;
363 let src = device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
364 label: Some("graphview_layout_src"),
365 contents: bytemuck::cast_slice(&nodes),
366 usage: node_usage,
367 });
368 let dst = device.create_buffer(&wgpu::BufferDescriptor {
369 label: Some("graphview_layout_dst"),
370 size: (std::mem::size_of::<GpuNode>() * nodes.len().max(1)) as u64,
371 usage: node_usage,
372 mapped_at_creation: false,
373 });
374 let off_data: &[u32] = if g.offsets.is_empty() { &[0] } else { &g.offsets };
376 let tgt_data: Vec<u32> = if g.targets.is_empty() { vec![0] } else { g.targets.clone() };
377 let offsets = device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
378 label: Some("graphview_layout_offsets"),
379 contents: bytemuck::cast_slice(off_data),
380 usage: wgpu::BufferUsages::STORAGE,
381 });
382 let targets = device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
383 label: Some("graphview_layout_targets"),
384 contents: bytemuck::cast_slice(&tgt_data),
385 usage: wgpu::BufferUsages::STORAGE,
386 });
387 let uniforms = device.create_buffer(&wgpu::BufferDescriptor {
388 label: Some("graphview_layout_uniforms"),
389 size: std::mem::size_of::<Uniforms>() as u64,
390 usage: wgpu::BufferUsages::UNIFORM | wgpu::BufferUsages::COPY_DST,
391 mapped_at_creation: false,
392 });
393 Self { pipeline, bgl, uniforms, offsets, targets, src, dst, count: nodes.len() as u32 }
394 }
395
396 pub fn iterate(&mut self, device: &wgpu::Device, queue: &wgpu::Queue, g: &GraphCsr, p: &LayoutParams, iters: usize) {
399 let u = Uniforms {
400 a: [self.count as f32, p.repulsion, p.attraction, p.ideal_len],
401 b: [p.gravity, p.damping, p.dt, p.min_dist],
402 };
403 queue.write_buffer(&self.uniforms, 0, bytemuck::bytes_of(&u));
404 let _ = g; let mut enc = device.create_command_encoder(&Default::default());
406 for _ in 0..iters {
407 let bind = device.create_bind_group(&wgpu::BindGroupDescriptor {
408 label: Some("graphview_layout_bind"),
409 layout: &self.bgl,
410 entries: &[
411 wgpu::BindGroupEntry { binding: 0, resource: self.uniforms.as_entire_binding() },
412 wgpu::BindGroupEntry { binding: 1, resource: self.src.as_entire_binding() },
413 wgpu::BindGroupEntry { binding: 2, resource: self.dst.as_entire_binding() },
414 wgpu::BindGroupEntry { binding: 3, resource: self.offsets.as_entire_binding() },
415 wgpu::BindGroupEntry { binding: 4, resource: self.targets.as_entire_binding() },
416 ],
417 });
418 {
419 let mut cp = enc.begin_compute_pass(&wgpu::ComputePassDescriptor { label: Some("graphview_layout_pass"), timestamp_writes: None });
420 cp.set_pipeline(&self.pipeline);
421 cp.set_bind_group(0, &bind, &[]);
422 cp.dispatch_workgroups(self.count.div_ceil(64), 1, 1);
423 }
424 std::mem::swap(&mut self.src, &mut self.dst);
425 }
426 queue.submit(Some(enc.finish()));
427 device.poll(wgpu::PollType::wait_indefinitely()).ok();
428 }
429
430 pub fn read_positions(&self, device: &wgpu::Device, queue: &wgpu::Queue) -> Vec<[f32; 2]> {
432 let size = (std::mem::size_of::<GpuNode>() as u32 * self.count) as u64;
433 let readback = device.create_buffer(&wgpu::BufferDescriptor {
434 label: Some("graphview_layout_readback"),
435 size,
436 usage: wgpu::BufferUsages::COPY_DST | wgpu::BufferUsages::MAP_READ,
437 mapped_at_creation: false,
438 });
439 let mut enc = device.create_command_encoder(&Default::default());
440 enc.copy_buffer_to_buffer(&self.src, 0, &readback, 0, size);
441 queue.submit(Some(enc.finish()));
442 let slice = readback.slice(..);
443 let (tx, rx) = std::sync::mpsc::channel();
444 slice.map_async(wgpu::MapMode::Read, move |r| { let _ = tx.send(r); });
445 device.poll(wgpu::PollType::wait_indefinitely()).ok();
446 rx.recv().unwrap().unwrap();
447 let data = slice.get_mapped_range();
448 let nodes: Vec<GpuNode> = bytemuck::cast_slice(&data).to_vec();
449 drop(data);
450 readback.unmap();
451 nodes.iter().map(|nd| [nd.a[0], nd.a[1]]).collect()
452 }
453 }
454}
455
456#[cfg(feature = "gpu")]
457pub use gpu::{GpuLayout, GpuNode, LAYOUT_WGSL};
458
459#[cfg(test)]
460mod tests {
461 use super::*;
462
463 fn two_triangles() -> GraphCsr {
465 let edges = [(0, 1), (1, 2), (2, 0), (3, 4), (4, 5), (5, 3), (2, 3)];
466 GraphCsr::from_edges(6, &edges, None)
467 }
468
469 #[test]
470 fn csr_is_symmetric_and_dedups() {
471 let g = GraphCsr::from_edges(3, &[(0, 1), (1, 0), (0, 1), (1, 2)], None);
473 assert_eq!(g.neighbours(0), &[1]);
474 assert_eq!(g.neighbours(1), &[0, 2]);
475 assert_eq!(g.neighbours(2), &[1]);
476 assert_eq!(g.degree_sum(), 4, "two undirected edges → 4 directed entries");
477 let g2 = GraphCsr::from_edges(2, &[(0, 0), (0, 5), (0, 1)], None);
479 assert_eq!(g2.neighbours(0), &[1]);
480 }
481
482 #[test]
483 fn seed_is_deterministic_and_bounded() {
484 let a = LayoutState::seed(50, 3.0);
485 let b = LayoutState::seed(50, 3.0);
486 assert_eq!(a, b, "no RNG: the seed is reproducible");
487 for q in &a.pos {
488 assert!(q[0].is_finite() && q[1].is_finite());
489 assert!((q[0] * q[0] + q[1] * q[1]).sqrt() <= 3.0 + 1e-3, "inside the seed disc");
490 }
491 }
492
493 #[test]
494 fn step_cpu_is_deterministic() {
495 let g = two_triangles();
496 let p = LayoutParams::default();
497 let s = LayoutState::seed(g.n, 2.0);
498 let a = step_cpu(&s, &g, &p);
499 let b = step_cpu(&s, &g, &p);
500 for (x, y) in a.pos.iter().zip(&b.pos) {
501 assert_eq!(x[0].to_bits(), y[0].to_bits(), "bit-identical pos.x");
502 assert_eq!(x[1].to_bits(), y[1].to_bits(), "bit-identical pos.y");
503 }
504 }
505
506 #[test]
507 fn layout_converges_energy_decays() {
508 let g = two_triangles();
509 let p = LayoutParams::default();
510 let mut s = LayoutState::seed(g.n, p.ideal_len * (g.n as f32).sqrt());
511 for _ in 0..15 {
514 s = step_cpu(&s, &g, &p);
515 }
516 let e_early = s.kinetic_energy();
517 for _ in 0..200 {
518 s = step_cpu(&s, &g, &p);
519 }
520 let e_late = s.kinetic_energy();
521 assert!(e_late < e_early * 0.1, "energy decayed an order of magnitude (converged): {e_early} → {e_late}");
522 assert!(e_late < 0.05, "settled to near-rest: {e_late}");
523 for q in &s.pos {
524 assert!(q[0].is_finite() && q[1].is_finite(), "no NaN/inf blow-up");
525 }
526 }
527
528 #[test]
529 fn springs_settle_an_edge_at_the_analytic_equilibrium() {
530 let g = GraphCsr::from_edges(2, &[(0, 1)], None);
534 let p = LayoutParams { gravity: 0.0, ..LayoutParams::default() };
535 let s = relax_cpu(&g, &p, 600);
536 let d = ((s.pos[1][0] - s.pos[0][0]).powi(2) + (s.pos[1][1] - s.pos[0][1]).powi(2)).sqrt();
537 let expected = (1.0 + 3.0_f32.sqrt()) / 2.0;
538 assert!((d - expected).abs() < 0.02, "edge rests at the force equilibrium {expected} (got {d})");
539 }
540
541 #[test]
542 fn relax_is_deterministic_positions_are_data() {
543 let g = two_triangles();
544 let p = LayoutParams::default();
545 let a = relax_cpu(&g, &p, 120);
546 let b = relax_cpu(&g, &p, 120);
547 assert_eq!(a.pos, b.pos, "same params → identical positions (reproducible DATA)");
548 let len = |i: usize, j: usize| ((a.pos[j][0] - a.pos[i][0]).powi(2) + (a.pos[j][1] - a.pos[i][1]).powi(2)).sqrt();
551 let within = len(0, 1).max(len(1, 2)).max(len(3, 4)).max(len(4, 5));
552 assert!(len(2, 3) > within * 0.9, "the two clusters spread apart along the bridge");
553 }
554
555 #[test]
556 fn heavier_nodes_move_less_per_step() {
557 let edges = [(0, 1), (0, 2)];
560 let light = GraphCsr::from_edges(3, &edges, None);
561 let heavy = GraphCsr::from_edges(3, &edges, Some(&[50.0, 1.0, 1.0]));
562 let p = LayoutParams::default();
563 let s = LayoutState::seed(3, 2.0);
564 let dl = step_cpu(&s, &light, &p);
565 let dh = step_cpu(&s, &heavy, &p);
566 let disp = |st: &LayoutState| ((st.pos[0][0] - s.pos[0][0]).powi(2) + (st.pos[0][1] - s.pos[0][1]).powi(2)).sqrt();
567 assert!(disp(&dh) < disp(&dl) * 0.5, "the heavy node barely moves ({} vs {})", disp(&dh), disp(&dl));
568 }
569
570 #[cfg(feature = "gpu")]
573 fn compute_device() -> Option<(wgpu::Device, wgpu::Queue)> {
574 let instance = wgpu::Instance::default();
575 let adapter = pollster::block_on(instance.request_adapter(&wgpu::RequestAdapterOptions {
576 power_preference: wgpu::PowerPreference::default(),
577 force_fallback_adapter: false,
578 compatible_surface: None,
579 }))
580 .ok()?;
581 if !adapter.get_downlevel_capabilities().flags.contains(wgpu::DownlevelFlags::COMPUTE_SHADERS) {
582 return None;
583 }
584 if adapter.limits().max_storage_buffers_per_shader_stage < 4 {
585 return None;
586 }
587 pollster::block_on(adapter.request_device(&wgpu::DeviceDescriptor {
588 label: Some("graphview-layout-proof"),
589 required_features: wgpu::Features::empty(),
590 required_limits: adapter.limits(),
591 memory_hints: wgpu::MemoryHints::default(),
592 experimental_features: wgpu::ExperimentalFeatures::disabled(),
593 trace: wgpu::Trace::Off,
594 }))
595 .ok()
596 }
597
598 #[cfg(feature = "gpu")]
603 #[test]
604 fn gpu_layout_matches_cpu_reference() {
605 let Some((device, queue)) = compute_device() else {
606 eprintln!("[gpu_layout] no compute device — skipping GPU parity proof");
607 return;
608 };
609 let g = two_triangles();
610 let p = LayoutParams::default();
611 let iters = 50;
612 let seed = LayoutState::seed(g.n, p.ideal_len * (g.n as f32).sqrt());
613
614 let mut gl = GpuLayout::new(&device, &g, &seed);
615 gl.iterate(&device, &queue, &g, &p, iters);
616 let gpu = gl.read_positions(&device, &queue);
617
618 let mut cpu = seed.clone();
619 for _ in 0..iters {
620 cpu = step_cpu(&cpu, &g, &p);
621 }
622 assert_eq!(gpu.len(), cpu.pos.len());
623 let mut moved = false;
624 for (gp, cp) in gpu.iter().zip(&cpu.pos) {
625 assert!((gp[0] - cp[0]).abs() < 1e-2, "pos.x parity: gpu {} vs cpu {}", gp[0], cp[0]);
626 assert!((gp[1] - cp[1]).abs() < 1e-2, "pos.y parity: gpu {} vs cpu {}", gp[1], cp[1]);
627 assert!(gp[0].is_finite() && gp[1].is_finite(), "readback is finite DATA");
628 if (gp[0] - seed.pos[0][0]).abs() > 1e-3 {
629 moved = true;
630 }
631 }
632 assert!(moved, "the VRAM iteration actually advanced the layout");
633 }
634
635 #[cfg(feature = "gpu")]
638 #[test]
639 fn layout_shader_entry_and_node_size() {
640 assert!(LAYOUT_WGSL.contains("fn cs_layout"));
641 assert!(LAYOUT_WGSL.contains("@workgroup_size(64)"));
642 assert_eq!(std::mem::size_of::<GpuNode>(), 32);
643 }
644
645 #[cfg(feature = "arrow")]
646 #[test]
647 fn csr_streams_from_an_arrow_edge_batch() {
648 use arrow_array::{Float64Array, Int64Array};
649 let src = Int64Array::from(vec![0, 1, 2, 0]);
650 let dst = Int64Array::from(vec![1, 2, 0, 2]);
651 let w = Float64Array::from(vec![3.0, 1.0, 1.0]);
652 let g = GraphCsr::from_arrow(3, &src, &dst, Some(&w));
653 assert_eq!(g.neighbours(0), &[1, 2]);
654 assert_eq!(g.weights[0], 3.0, "weight column streamed in");
655 let s = relax_cpu(&g, &LayoutParams::default(), 80);
657 assert!(s.pos.iter().all(|q| q[0].is_finite() && q[1].is_finite()));
658 }
659}