facett-graphview/src/layout_grid.wgsl
// **GPU force layout, grid-pyramid repulsion** — the O(n · log n) twin of `layout.wgsl`.
//
// `cs_layout` in `layout.wgsl` sums repulsion over ALL pairs: n² per step. At 100 000
// nodes that is 10^10 pair terms per step; a discrete GPU pays ~18 ms for it, an
// integrated one pays seconds. This file replaces the all-pairs loop with a regular
// quad-tree built fresh every step (the grid pyramid): nodes are binned into a
// G × G grid (G a power of two), coarser levels are sums of four children, and each
// node takes repulsion from
//
// * the 3 × 3 finest cells around its own (own cell with itself subtracted), and
// * at every level, the up-to-27 cells that are children of its parent's 3 × 3
// neighbourhood but outside its own 3 × 3 neighbourhood (the classic
// interaction list — every far cell is at least one cell width away, θ ≤ 1).
//
// A cell acts as its mass at its centroid. The spring, gravity and integration terms
// are line-for-line the ones in `cs_layout`, so the two kernels differ ONLY in how
// repulsion is summed.
//
// Passes per step (encoded by `GpuLayout::submit_steps`):
// 1. clear `fine` and `bounds` (buffer clears, no shader)
// 2. cs_bounds — workgroup min/max of positions → 4 atomics
// 3. cs_bin — each node adds (mass, offset·mass) into its finest cell
// 4. cs_fine — finest cells → level K of `cells` (mass-weighted sums)
// 5. cs_reduce ×K — level l = sum of four children at level l+1 (dynamic offset)
// 6. cs_force — repulsion from the pyramid + springs + gravity + integrate
struct LNode {
a: vec4<f32>,
b: vec4<f32>,
};
struct U {
// x = n, y = repulsion, z = attraction, w = ideal_len
a: vec4<f32>,
// x = gravity, y = damping, z = dt, w = min_dist
b: vec4<f32>,
// x = G (finest side), y = K (= log2 G), z, w = unused
c: vec4<f32>,
};
struct L {
// x = the level this cs_reduce dispatch WRITES (its children are level x + 1)
level: vec4<u32>,
};
@group(0) @binding(0) var<uniform> P: U;
@group(0) @binding(1) var<storage, read> src: array<LNode>;
@group(0) @binding(2) var<storage, read_write> dst: array<LNode>;
@group(0) @binding(3) var<storage, read> offsets: array<u32>;
@group(0) @binding(4) var<storage, read> targets: array<u32>;
// 4 atomics: max of ordered(-x), ordered(-y), ordered(x), ordered(y). Cleared to 0.
@group(0) @binding(5) var<storage, read_write> bounds: array<atomic<u32>, 4>;
// 3 per finest cell: mass·16, Σ ux·mass·16·1024, Σ uy·mass·16·1024 (u = offset in cell).
@group(0) @binding(6) var<storage, read_write> fine: array<atomic<u32>>;
// The pyramid, level l at offset (4^l − 1) / 3: (Σ m·x, Σ m·y, Σ m, 0).
@group(0) @binding(7) var<storage, read_write> cells: array<vec4<f32>>;
@group(1) @binding(0) var<uniform> LV: L;
const MASS_Q: f32 = 16.0;
const OFF_Q: f32 = 1024.0;
// Order-preserving f32 → u32: larger float ⇒ larger u32, so atomicMax works on floats.
fn ordered(x: f32) -> u32 {
let b = bitcast<u32>(x);
if ((b & 0x80000000u) != 0u) { return ~b; }
return b | 0x80000000u;
}
fn unordered(u: u32) -> f32 {
if ((u & 0x80000000u) != 0u) { return bitcast<f32>(u & 0x7fffffffu); }
return bitcast<f32>(~u);
}
fn level_base(l: u32) -> u32 {
// (4^l − 1) / 3
return ((1u << (2u * l)) - 1u) / 3u;
}
// The square domain every pass agrees on this step: (min.x, min.y, side).
fn domain() -> vec3<f32> {
let mnx = -unordered(atomicLoad(&bounds[0]));
let mny = -unordered(atomicLoad(&bounds[1]));
let mxx = unordered(atomicLoad(&bounds[2]));
let mxy = unordered(atomicLoad(&bounds[3]));
let side = max(max(mxx - mnx, mxy - mny), 1e-3) * 1.0001;
return vec3<f32>(mnx, mny, side);
}
var<workgroup> wg: array<vec4<f32>, 256>;
@compute @workgroup_size(256)
fn cs_bounds(@builtin(global_invocation_id) gid: vec3<u32>, @builtin(local_invocation_index) li: u32) {
let n = u32(P.a.x);
var v = vec4<f32>(-3.4e38, -3.4e38, -3.4e38, -3.4e38);
if (gid.x < n) {
let p = src[gid.x].a.xy;
v = vec4<f32>(-p.x, -p.y, p.x, p.y);
}
wg[li] = v;
workgroupBarrier();
var s = 128u;
loop {
if (s == 0u) { break; }
if (li < s) { wg[li] = max(wg[li], wg[li + s]); }
workgroupBarrier();
s = s >> 1u;
}
if (li == 0u) {
let r = wg[0];
atomicMax(&bounds[0], ordered(r.x));
atomicMax(&bounds[1], ordered(r.y));
atomicMax(&bounds[2], ordered(r.z));
atomicMax(&bounds[3], ordered(r.w));
}
}
fn cell_of(p: vec2<f32>, dom: vec3<f32>, side_cells: u32) -> vec2<u32> {
let f = (p - dom.xy) / dom.z * f32(side_cells);
let m = f32(side_cells - 1u);
return vec2<u32>(u32(clamp(f.x, 0.0, m)), u32(clamp(f.y, 0.0, m)));
}
@compute @workgroup_size(256)
fn cs_bin(@builtin(global_invocation_id) gid: vec3<u32>) {
let n = u32(P.a.x);
let i = gid.x;
if (i >= n) { return; }
let g = u32(P.c.x);
let dom = domain();
let me = src[i];
let p = me.a.xy;
let c = cell_of(p, dom, g);
let cs = dom.z / f32(g);
let u = clamp((p - (dom.xy + vec2<f32>(c) * cs)) / cs, vec2<f32>(0.0), vec2<f32>(0.999));
let mq = max(u32(round(me.b.x * MASS_Q)), 1u);
let k = (c.y * g + c.x) * 3u;
atomicAdd(&fine[k], mq);
atomicAdd(&fine[k + 1u], u32(u.x * OFF_Q) * mq);
atomicAdd(&fine[k + 2u], u32(u.y * OFF_Q) * mq);
}
@compute @workgroup_size(256)
fn cs_fine(@builtin(global_invocation_id) gid: vec3<u32>) {
let g = u32(P.c.x);
let kk = u32(P.c.y);
let idx = gid.x;
if (idx >= g * g) { return; }
let mq = atomicLoad(&fine[idx * 3u]);
var out = vec4<f32>(0.0);
if (mq > 0u) {
let dom = domain();
let cs = dom.z / f32(g);
let cx = idx % g;
let cy = idx / g;
let m = f32(mq) / MASS_Q;
let ux = f32(atomicLoad(&fine[idx * 3u + 1u])) / (f32(mq) * OFF_Q);
let uy = f32(atomicLoad(&fine[idx * 3u + 2u])) / (f32(mq) * OFF_Q);
let c = dom.xy + (vec2<f32>(f32(cx), f32(cy)) + vec2<f32>(ux, uy)) * cs;
out = vec4<f32>(c * m, m, 0.0);
}
cells[level_base(kk) + idx] = out;
}
@compute @workgroup_size(256)
fn cs_reduce(@builtin(global_invocation_id) gid: vec3<u32>) {
let l = LV.level.x;
let side = 1u << l;
let idx = gid.x;
if (idx >= side * side) { return; }
let x = idx % side;
let y = idx / side;
let cside = side * 2u;
let cb = level_base(l + 1u);
let c00 = cells[cb + (2u * y) * cside + 2u * x];
let c10 = cells[cb + (2u * y) * cside + 2u * x + 1u];
let c01 = cells[cb + (2u * y + 1u) * cside + 2u * x];
let c11 = cells[cb + (2u * y + 1u) * cside + 2u * x + 1u];
cells[level_base(l) + idx] = c00 + c10 + c01 + c11;
}
fn push(pos_i: vec2<f32>, cell: vec4<f32>, k: f32, min2: f32) -> vec2<f32> {
if (cell.z <= 0.0) { return vec2<f32>(0.0); }
let c = cell.xy / cell.z;
let d = pos_i - c;
let dist2 = max(dot(d, d), min2);
return d * (k * cell.z / dist2);
}
@compute @workgroup_size(64)
fn cs_force(@builtin(global_invocation_id) gid: vec3<u32>) {
let i = gid.x;
let n = u32(P.a.x);
if (i >= n) { return; }
let rep = P.a.y;
let attr = P.a.z;
let ideal = P.a.w;
let grav = P.b.x;
let damp = P.b.y;
let dt = P.b.z;
let mind = P.b.w;
let min2 = mind * mind;
let g = u32(P.c.x);
let kk = u32(P.c.y);
let me = src[i];
let pos_i = me.a.xy;
let wi = me.b.x;
let k = rep * wi;
let dom = domain();
var force = vec2<f32>(0.0, 0.0);
// ── far field: the interaction list at every level ──
for (var l = 1u; l <= kk; l = l + 1u) {
let side = 1u << l;
let mine = cell_of(pos_i, dom, side);
let par = mine / 2u;
let base = level_base(l);
let px0 = max(i32(par.x) - 1, 0);
let py0 = max(i32(par.y) - 1, 0);
let px1 = min(i32(par.x) + 1, i32(side / 2u) - 1);
let py1 = min(i32(par.y) + 1, i32(side / 2u) - 1);
for (var y = py0 * 2; y <= py1 * 2 + 1; y = y + 1) {
for (var x = px0 * 2; x <= px1 * 2 + 1; x = x + 1) {
if (abs(x - i32(mine.x)) <= 1 && abs(y - i32(mine.y)) <= 1) { continue; }
force = force + push(pos_i, cells[base + u32(y) * side + u32(x)], k, min2);
}
}
}
// ── near field: the finest 3 × 3, own cell without itself ──
let mine = cell_of(pos_i, dom, g);
let base = level_base(kk);
for (var y = max(i32(mine.y) - 1, 0); y <= min(i32(mine.y) + 1, i32(g) - 1); y = y + 1) {
for (var x = max(i32(mine.x) - 1, 0); x <= min(i32(mine.x) + 1, i32(g) - 1); x = x + 1) {
var cell = cells[base + u32(y) * g + u32(x)];
if (x == i32(mine.x) && y == i32(mine.y)) {
// Subtract this node's own (quantised) contribution.
let cs = dom.z / f32(g);
let u = clamp((pos_i - (dom.xy + vec2<f32>(mine) * cs)) / cs, vec2<f32>(0.0), vec2<f32>(0.999));
let mq = f32(max(u32(round(wi * MASS_Q)), 1u)) / MASS_Q;
let own = dom.xy + (vec2<f32>(mine) + floor(u * OFF_Q) / OFF_Q) * cs;
cell = cell - vec4<f32>(own * mq, mq, 0.0);
if (cell.z < 0.5 / MASS_Q) { continue; }
}
force = force + push(pos_i, cell, k, min2);
}
}
// ── spring attraction along CSR neighbours (as in cs_layout) ──
let start = offsets[i];
let end = offsets[i + 1u];
for (var e = start; e < end; e = e + 1u) {
let t = targets[e];
let d = src[t].a.xy - pos_i;
let dist = max(length(d), 1e-4);
force = force + d * (attr * (dist - ideal) / dist);
}
// ── centring gravity ──
force = force - pos_i * grav;
// ── integrate (mass = wi) ──
let inv_m = 1.0 / max(wi, 1e-4);
let vel = (me.a.zw + force * inv_m * dt) * damp;
let pos = pos_i + vel * dt;
dst[i].a = vec4<f32>(pos.x, pos.y, vel.x, vel.y);
dst[i].b = me.b;
}