diff --git a/docs/many-body-force/README.md b/docs/many-body-force/README.md index 9746b0c6..813e5aaf 100644 --- a/docs/many-body-force/README.md +++ b/docs/many-body-force/README.md @@ -3,8 +3,18 @@ This is a walkthrough of the repulsion force introduced in [#240](https://github.com/cosmosgl/graph/pull/240) — what the old algorithm did, what the new one does instead, and why the change matters. The code lives in -`src/modules/ForceManyBody/`; the short engineering record is -`history/2026/2026-07-08-many-body-repulsion.md`. +`src/modules/ForceManyBody/`; the short engineering records are +`history/2026/2026-07-08-many-body-repulsion.md` and (for the sampling-noise follow-up +described near the end) `history/2026/2026-08-14-nearfield-jitter.md`. + +Since that follow-up the force has **two paths**: graphs of at most 4,096 points are computed +**exactly** — one all-pairs pass, no grid, no sampling (`force-allpairs.frag`, see +[Small graphs are exact](#small-graphs-are-exact-the-all-pairs-path)) — and everything below +describes the grid + Monte-Carlo machinery that larger graphs use. + +The problem-and-fix story behind that follow-up — the measurements, the captured trajectories, +and its own figures — is also a standalone page: [`jitter-fix.html`](./jitter-fix.html) +(self-contained; open it locally in a browser — GitHub shows only its source). ## The problem both algorithms solve @@ -96,10 +106,22 @@ computes an **unbiased estimate** instead: ![Depth peeling and Horvitz–Thompson weighting](c-depth-peeling.svg) **Sampling (`build-nearfield-slots.vert`):** every tick, each point gets a fresh pseudo-random -hash. Eight "depth peeling" passes then run over the finest grid; pass *k* selects, per cell, +hash. K "depth peeling" passes then run over the finest grid; pass *k* selects, per cell, the point with the smallest hash *not yet selected by passes 0..k−1* (the GPU depth test does -the per-cell minimum for free). After 8 passes, each cell's 8 **slot textures** hold a uniform -random 8-subset of its points — re-drawn from scratch every tick. +the per-cell minimum for free). After K passes, the K layers of a **slot array texture** +(`sampler2DArray`, one layer per pass) hold a uniform random K-subset of each cell's points — +re-drawn from scratch every tick. + +**K adapts to the graph size** (`getNearFieldSlotCount` in `index.ts`): 32 slots up to 16k +points, 16 up to 65k, 8 above. Peeling is inherently sequential — one render pass per slot — +so K is a direct trade between per-tick cost and sampling variance; big graphs keep the cheap +estimator (per-point noise is sub-pixel at that scale), smaller ones buy more slots so that +realistic hub-cell occupancies are covered exactly. Each pass ping-pongs between two plain 2D +targets (pass *k* must sample pass *k−1*'s output, and sampling one layer of a texture while +rendering to another layer of the same texture is a WebGL feedback loop) and its result is +copied into its array layer. + +![Peeling ping-pong into the slot array](h-pingpong-peel.svg) **The hash must be an integer hash.** The first version used the classic `fract(sin(index * 12.9898 + seed * 78.233) * 43758.5453)` one-liner, which is quietly broken at @@ -124,7 +146,7 @@ the next pass's comparison. Cost is a wash: ~8 integer ALU ops replace a special the true pairwise forces from the sampled slots (skipping itself), then scales the sum by ``` -others / sampled // e.g. cell holds 12 other points, 8 sampled → × 12/8 +others / sampled // e.g. cell holds 48 other points, 32 sampled → × 48/32 ``` This is the **Horvitz–Thompson estimator**: since each of the cell's `others` points had equal @@ -136,16 +158,16 @@ centroid term**, so the tangential force component survives: Two properties fall out for free: -- **Sparse cells are exact.** A cell with ≤ 8 points is sampled exhaustively +- **Sparse cells are exact.** A cell with ≤ K points is sampled exhaustively (`others == sampled`, weight = 1). With the finest grid at 2·√n per axis, the average cell holds ¼ point — so for typical graphs the near field *is* the exact all-pairs force, and the - approximation only kicks in inside genuinely dense hubs. This is also why the prototyped - separate brute-force path for small graphs was dropped: the grid path is already effectively - exact there, and faster. -- **The sampling noise is a feature.** The estimate is unbiased but noisy, and the noise is - re-rolled every tick. Because every force is scaled by `alpha`, the noise *anneals*: large - early, when clumps need breaking apart, shrinking to nothing as the layout settles. It is - precisely the jitter that lets stacked points find distinct directions to escape along. + approximation only kicks in inside genuinely dense hubs. +- **The sampling noise is (mostly) a feature.** The estimate is unbiased but noisy, and the + noise is re-rolled every tick. Because every force is scaled by `alpha`, the noise *anneals*: + large early, when clumps need breaking apart, shrinking to nothing as the layout settles. It + is precisely the jitter that lets stacked points find distinct directions to escape along. + The caveat — and what the 2026-08-14 follow-up fixed — is layouts where the density never + disperses: see [When the noise stops annealing](#when-the-noise-stops-annealing) below. ### Step 4 — two stability guards @@ -157,11 +179,11 @@ Both live in `force-nearfield.frag`, both born from real failure modes: ring" around the stack. Instead, each point kicks along its own per-point random vector, so a pile disperses. - **Per-tick velocity clamp (2 × cell size).** The `others/sampled` weight is unbiased but - high-variance: in a cell holding far more points than 8 slots, a couple of very close samples - can be multiplied into a huge one-tick kick — flinging points across the screen at startup - and ejecting points from dense cluster centers. The clamp caps the magnitude and keeps the - direction; genuine spreading kicks are far below the bound, and bulk expansion is driven by - the far-field levels anyway. + high-variance: in a cell holding far more points than sampling slots, a couple of very close + samples can be multiplied into a huge one-tick kick — flinging points across the screen at + startup and ejecting points from dense cluster centers. The clamp caps the magnitude and + keeps the direction; genuine spreading kicks are far below the bound, and bulk expansion is + driven by the far-field levels anyway. ### The per-tick pipeline @@ -170,7 +192,60 @@ Both live in `force-nearfield.frag`, both born from real failure modes: Orchestrated by `src/modules/ForceManyBody/index.ts`: `drawLevels()` → `drawNearFieldSlots()` → `drawForces()` (per-level force passes plus the near-field pass, all blending additively into the shared velocity texture). Integration into -positions is the same step every other force uses. +positions is the same step every other force uses. Small graphs replace all three with a +single `drawAllPairsForce()` pass — next section. + +## When the noise stops annealing + +The "noise is a feature" argument has a hole, found the hard way on a real graph (the +163-country border-adjacency network): it assumes the density that causes the sampling +variance *disperses*. Under pure repulsion it does — the clump expands, occupancy falls to ~K +within a second, and the noise dies with it, sub-pixel before anyone sees it. But when link +attraction or gravity holds a hub together *while alpha stays high* (a long-running layout, a +reheated one, `start()` on interaction), occupancy stays far above the slot count forever, and +the per-tick re-drawn sample turns into visible, permanent shimmer: measured on that country +graph, every point wandered ~0.5 units per tick with a ~92° mean direction change — a pure +random walk stacked on a settled layout, while an exact all-pairs reference under the same +integration was three orders of magnitude stiller. + +![Same cell, a fresh sample every tick](f-resample-jitter.svg) + +Two changes closed it (2026-08-14): + +1. **Small graphs skip the estimator entirely** — the all-pairs path below. This is the case + where sustained dense hubs are both most common and cheapest to compute exactly. +2. **K became adaptive** (32/16/8 — see step 3). Mid-size graphs get 2–4× more samples, which + both extends the exactly-covered occupancy range and shrinks the residual variance + (amplitude ∝ occupancy/K · 1/√K) — while the ≥ 65k tier keeps today's 8-slot cost + unchanged. + +The **Performance → Repulsion Jitter: Fixed vs Before** story shows this live: the real +graph that surfaced it, run side by side — today's exact path next to the pre-fix sampled +configuration (forced back on through a story-only internals patch) — each with alpha held +at 1, a sliding-window step/turn meter, and a trajectory trace of one dense-cell point. + +## Small graphs are exact: the all-pairs path + +At or below **4,096 points** (`ALL_PAIRS_MAX_POINTS`) the force runs as one full-screen pass +(`force-allpairs.frag`): each point loops over every other point and sums the same clamped +inverse-distance pairwise force the grid path uses, with the same coincident-point random +kick. No pyramid, no peeling, no sampling — the result is exact at *any* cell occupancy, so +there is no noise to anneal and nothing to shimmer. The two paths also bound the same thing: +pairs within the near-field scale (2 × the finest cell size the grid path would use) are +jittered and capped like the near-field pass's sum, so a coincident stack expands instead of +teleporting, while farther pairs pass through unbounded like the level passes. + +![The two paths](g-two-paths.svg) + +It is also simply faster there. Depth peeling costs one render pass per slot, ~0.1 ms of +fixed overhead each; at 2k points, 64 experimental slots measured ~6.4 ms/step while the +single all-pairs pass measures **~1.8 ms/step** — n² texel loops are trivial work at this +scale (4096² ≈ 17M pair evaluations). The threshold sits where that stops being true: the +next power of two would already cost several milliseconds. + +(#240 prototyped and dropped exactly this path, when the grid looked "effectively exact" for +small graphs. That held for dispersing layouts; the sustained-density case above is why it +came back.) ## Why it is better than the old one @@ -180,15 +255,16 @@ positions is the same step every other force uses. | Dense hubs | collapse into disks / petals | spread into natural clouds | | Stacked points | never separate (void rings) | random kick disperses them | | Coverage seams | depend on `theta` tuning | exact once-tiling, nothing to tune | -| Small / sparse graphs | always approximate | effectively **exact** (cells ≤ 8 points) | +| Small graphs (≤ 4,096 points) | always approximate | **exact** — dedicated all-pairs pass | +| Sparse cells (larger graphs) | always approximate | **exact** (cells ≤ K points, K = 8–32) | | Bias | systematic (centroid direction) | none — unbiased estimator; noise anneals with alpha | | `simulationRepulsionTheta` | required tuning | deprecated no-op (accepted, ignored) | | Speed | baseline | **~1.2–4× faster per step** across the practical range | -| Code paths | one, plus the theta special-casing | one, for every graph size | +| Code paths | one, plus the theta special-casing | two: exact ≤ 4k points, grid + sampling above | The speedup comes from the fixed 3×3/6×6 loop structure (compact, coherent texel fetches; -no data-dependent band walking) — measure it yourself with the **Misc → Repulsion Benchmark** -story, which steps the simulation directly and forces a readback so the numbers aren't capped +no data-dependent band walking) — measure it yourself with the **Performance → Repulsion +Benchmark** story, which steps the simulation directly and forces a readback so the numbers aren't capped by the display refresh rate. The before/after videos in the [PR description](https://github.com/cosmosgl/graph/pull/240) show the visual difference on a dense-hub graph. @@ -197,13 +273,16 @@ dense-hub graph. The old algorithm's error was a systematic *bias* — invisible per-frame, but it deformed layouts (disks, petals, void rings) and never went away. The new algorithm's error is -*variance* — visible as per-tick jitter in dense regions while `alpha` is high, but centered on -the exact answer and vanishing as the simulation cools. For a layout engine that trade is -clearly right: the final layout is what users keep, and the transient jitter is doing useful -work (annealing) on the way there. - -Cost side: 8 extra slot textures at the finest grid resolution (at the 512² cap that is -8 × 512² × 2 floats = 16 MB of GPU memory) and the 8 peeling passes per tick — both already +*variance* — per-tick jitter in dense regions while `alpha` is high, centered on the exact +answer and vanishing as the simulation cools. For a layout engine that trade is right: the +final layout is what users keep, and the transient jitter is doing useful work (annealing) on +the way there. Where the variance stopped vanishing — sustained-density layouts — it was +removed outright (exact ≤ 4k points) or shrunk (adaptive K); the residual is confined to +over-K-occupancy hub cells in graphs above 4k points while alpha is high. + +Cost side: the slot array texture plus two peel targets at the finest grid resolution (at the +512² cap with K = 8 that is 10 × 512² × 2 floats = 20 MB of GPU memory; smaller graphs have +more layers but a proportionally smaller grid) and the K peeling passes per tick — all already included in the benchmark numbers above. ## Glossary diff --git a/docs/many-body-force/c-depth-peeling.svg b/docs/many-body-force/c-depth-peeling.svg index 7cc1d2b3..2cf69a3e 100644 --- a/docs/many-body-force/c-depth-peeling.svg +++ b/docs/many-body-force/c-depth-peeling.svg @@ -5,5 +5,5 @@ -One finest-level cell, one tick13 points, each hashed with this tick’s random seed.62.77.56.97.91.81.94.74.64.20.93.12.82sampled (8 smallest hashes)left out8 depth-peeling passes → 8 slot texturespass k keeps the smallest hash not yet peeledslot 0point #11hash 0.12slot 1point #9hash 0.20slot 2point #2hash 0.56slot 3point #0hash 0.62slot 4point #8hash 0.64slot 5point #7hash 0.74slot 6point #1hash 0.77slot 7point #5hash 0.81Horvitz–Thompson weightingsampled sum = full sum, on averagecell has 12 other points,8 of them sampled →F ≈ (12 / 8) · Σ F(sampled pair)E[F] = exact all-pairs sum(unbiased, no centroid term)≤ 8 points in cell → exact.A fresh random subset every tick:the sampling noise shrinks with alphaand acts as annealing jitter. +One finest-level cell, one tick13 points, each hashed with this tick’s random seed.62.77.56.97.91.81.94.74.64.20.93.12.82sampled (the K smallest hashes)left outK depth-peeling passes → K slot array layers (K = 8 drawn)pass k keeps the smallest hash not yet peeledslot 0point #11hash 0.12slot 1point #9hash 0.20slot 2point #2hash 0.56slot 3point #0hash 0.62slot 4point #8hash 0.64slot 5point #7hash 0.74slot 6point #1hash 0.77slot 7point #5hash 0.81Horvitz–Thompson weightingsampled sum = full sum, on averagecell has 12 other points,8 of them sampled →F ≈ (12 / 8) · Σ F(sampled pair)E[F] = exact all-pairs sum(unbiased, no centroid term)≤ K points in cell → exact.A fresh random subset every tick: thenoise anneals with alpha while densitydisperses. Sustained density is why Kadapts (32/16/8) and ≤ 4k graphs go exact. \ No newline at end of file diff --git a/docs/many-body-force/d-gpu-pipeline.svg b/docs/many-body-force/d-gpu-pipeline.svg index e7279682..816b206c 100644 --- a/docs/many-body-force/d-gpu-pipeline.svg +++ b/docs/many-body-force/d-gpu-pipeline.svg @@ -1,9 +1,9 @@ - + - -One simulation tick of the repulsion force (all on the GPU)positionsone texel per point(x, y)1 · aggregate levelsdraw n points into each grid,additive blend accumulates[Σx, Σy, count] per cell2 · build near-field slots8 depth-peel passes over thefinest grid: a fresh random8-subset per cell, every tick3 · force passesper level: centroid repulsion+ near field: weighted pairs;all add into the velocity textureThe integration step (velocity → positions) is shared with all other forces and unchanged.Levels: 4², 8², … up to ≈ 2·√n per axis (capped 512²). Slot textures: 8 × finest grid, [point index, hash] each. + +One simulation tick of the repulsion force (all on the GPU)positionsone texel per point(x, y)1 · aggregate levelsdraw n points into each grid,additive blend accumulates[Σx, Σy, count] per cell2 · build near-field slotsK depth-peel passes over thefinest grid: a fresh randomK-subset per cell (K = 32/16/8)3 · force passesper level: centroid repulsion+ near field: weighted pairs;all add into the velocity textureThe integration step (velocity → positions) is shared with all other forces and unchanged.Levels: 4², 8², … up to ≈ 2·√n per axis (capped 512²). Slot array: K layers × finest grid, [point index, hash] each.Graphs of ≤ 4,096 points skip all three passes: a single exact all-pairs pass writes the velocity texture directly. \ No newline at end of file diff --git a/docs/many-body-force/f-resample-jitter.svg b/docs/many-body-force/f-resample-jitter.svg new file mode 100644 index 00000000..489b9801 --- /dev/null +++ b/docs/many-body-force/f-resample-jitter.svg @@ -0,0 +1,9 @@ + + + + + + + +Same cell, same positions — a fresh sample every ticktick t — sample A (amber)force on ⊙ = Σ over sample A × m/Ktick t+1 — sample B (teal)same positions, different draw → different answerthe point, tick after tick≈0.5 units/tick, ~92° mean turn — at equilibriumThe estimate is unbiased — averaged over ticks it equals the exact force — but each tick draws a fresh K-subset andweights it by m/K. While the cell stays dense (m ≫ K) that per-tick difference never shrinks: the noise is not annealing,it is a permanent random walk stacked on a settled layout. + \ No newline at end of file diff --git a/docs/many-body-force/g-two-paths.svg b/docs/many-body-force/g-two-paths.svg new file mode 100644 index 00000000..ff6cc9b2 --- /dev/null +++ b/docs/many-body-force/g-two-paths.svg @@ -0,0 +1,9 @@ + + + + + + + +Two paths — the point count picks onegraph ofn pointsn ≤ 4096usesAllPairsnoaggregate thegrid pyramiddepth-peel K slotsper finest cellK = 32 (≤16k) · 16 (≤65k) · 8 abovefar field per level +sampled near field ×m/Kyesforce-allpairs.fragone exact O(n²) pass — no grid, no samplingvelocitytextureadditivesingle writeSame pairwise falloff and coincident-point kick on both paths — the small-graph path just computes every pair instead of estimating. + \ No newline at end of file diff --git a/docs/many-body-force/gen-diagrams.cjs b/docs/many-body-force/gen-diagrams.cjs index 965017de..eaf7c83e 100644 --- a/docs/many-body-force/gen-diagrams.cjs +++ b/docs/many-body-force/gen-diagrams.cjs @@ -152,7 +152,7 @@ const arrow = (x1, y1, x2, y2, stroke, width = 2, dash = '') => } // ---------------------------------------------------------------- Diagram C -// Depth peeling K=8 random slots per cell + Horvitz–Thompson weighting. +// Depth peeling K random slots per cell (K = 8 drawn) + Horvitz–Thompson weighting. { const W = 980; const H = 430 let b = '' @@ -177,13 +177,13 @@ const arrow = (x1, y1, x2, y2, stroke, width = 2, dash = '') => b += `` b += txt(x, y + 24, h.toFixed(2).slice(1), { size: 10.5, fill: isP ? C.covered[0] : C.sub, anchor: 'middle' }) }) - b += txt(cellX, cellY + cellS + 28, 'sampled (8 smallest hashes)', { size: 12, fill: C.covered[0] }) + b += txt(cellX, cellY + cellS + 28, 'sampled (the K smallest hashes)', { size: 12, fill: C.covered[0] }) b += `` b += txt(cellX + 202, cellY + cellS + 28, 'left out', { size: 12, fill: C.sub }) // Peeling passes → slot textures const sx = 380; const sy = 78 - b += txt(sx, 40, '8 depth-peeling passes → 8 slot textures', { size: 15, weight: '600' }) + b += txt(sx, 40, 'K depth-peeling passes → K slot array layers (K = 8 drawn)', { size: 15, weight: '600' }) b += txt(sx, 58, 'pass k keeps the smallest hash not yet peeled', { size: 12, fill: C.sub }) order.slice(0, 8).forEach((o, k) => { const y = sy + k * 34 @@ -203,17 +203,18 @@ const arrow = (x1, y1, x2, y2, stroke, width = 2, dash = '') => b += txt(ex + 16, ey + 88, 'F ≈ (12 / 8) · Σ F(sampled pair)', { size: 14.5, weight: '600', fill: C.accent }) b += txt(ex + 16, ey + 122, 'E[F] = exact all-pairs sum', { size: 13.5 }) b += txt(ex + 16, ey + 142, '(unbiased, no centroid term)', { size: 12.5, fill: C.sub }) - b += txt(ex + 16, ey + 166, '≤ 8 points in cell → exact.', { size: 13, weight: '600' }) - b += txt(ex, ey + 210, 'A fresh random subset every tick:', { size: 12.5, fill: C.sub }) - b += txt(ex, ey + 228, 'the sampling noise shrinks with alpha', { size: 12.5, fill: C.sub }) - b += txt(ex, ey + 246, 'and acts as annealing jitter.', { size: 12.5, fill: C.sub }) + b += txt(ex + 16, ey + 166, '≤ K points in cell → exact.', { size: 13, weight: '600' }) + b += txt(ex, ey + 210, 'A fresh random subset every tick: the', { size: 12.5, fill: C.sub }) + b += txt(ex, ey + 228, 'noise anneals with alpha while density', { size: 12.5, fill: C.sub }) + b += txt(ex, ey + 246, 'disperses. Sustained density is why K', { size: 12.5, fill: C.sub }) + b += txt(ex, ey + 264, 'adapts (32/16/8) and ≤ 4k graphs go exact.', { size: 12.5, fill: C.sub }) fs.writeFileSync(path.join(OUT, 'c-depth-peeling.svg'), svgDoc(W, H, b)) } // ---------------------------------------------------------------- Diagram D // The per-tick GPU pipeline. { - const W = 980; const H = 300 + const W = 980; const H = 322 let b = '' const box = (x, y, w, h, title, lines, color, soft) => { let s = `` @@ -227,11 +228,12 @@ const arrow = (x1, y1, x2, y2, stroke, width = 2, dash = '') => b += arrow(170, midY + 55, 208, midY + 55, C.frame, 2.5) b += box(210, midY, 220, 110, '1 · aggregate levels', ['draw n points into each grid,', 'additive blend accumulates', '[Σx, Σy, count] per cell'], C.covered[1], C.coveredSoft[1]) b += arrow(430, midY + 55, 468, midY + 55, C.frame, 2.5) - b += box(470, midY, 220, 110, '2 · build near-field slots', ['8 depth-peel passes over the', 'finest grid: a fresh random', '8-subset per cell, every tick'], C.covered[0], C.coveredSoft[0]) + b += box(470, midY, 220, 110, '2 · build near-field slots', ['K depth-peel passes over the', 'finest grid: a fresh random', 'K-subset per cell (K = 32/16/8)'], C.covered[0], C.coveredSoft[0]) b += arrow(690, midY + 55, 728, midY + 55, C.frame, 2.5) b += box(730, midY, 230, 110, '3 · force passes', ['per level: centroid repulsion', '+ near field: weighted pairs;', 'all add into the velocity texture'], C.covered[2], C.coveredSoft[2]) b += txt(20, 250, 'The integration step (velocity → positions) is shared with all other forces and unchanged.', { size: 12.5, fill: C.sub }) - b += txt(20, 272, 'Levels: 4², 8², … up to ≈ 2·√n per axis (capped 512²). Slot textures: 8 × finest grid, [point index, hash] each.', { size: 12.5, fill: C.sub }) + b += txt(20, 272, 'Levels: 4², 8², … up to ≈ 2·√n per axis (capped 512²). Slot array: K layers × finest grid, [point index, hash] each.', { size: 12.5, fill: C.sub }) + b += txt(20, 294, 'Graphs of ≤ 4,096 points skip all three passes: a single exact all-pairs pass writes the velocity texture directly.', { size: 12.5, fill: C.sub }) fs.writeFileSync(path.join(OUT, 'd-gpu-pipeline.svg'), svgDoc(W, H, b)) } @@ -277,4 +279,121 @@ const arrow = (x1, y1, x2, y2, stroke, width = 2, dash = '') => fs.writeFileSync(path.join(OUT, 'e-old-theta-bands.svg'), svgDoc(W, H, b)) } +// ---------------------------------------------------------------- Diagram F +// The sustained-density failure mode: a fresh sample every tick means the +// force on a confined point is re-rolled noise, not annealing jitter. +{ + const W = 980; const H = 430 + let b = '' + b += txt(20, 34, 'Same cell, same positions — a fresh sample every tick', { size: 16, weight: '600' }) + // One shared dot layout (relative to a 200×200 cell): sample A one tick, + // sample B the next, one dot never drawn either tick, plus the observed point. + const setA = [[60, 50], [120, 55], [75, 80], [165, 100], [115, 105], [70, 140], [130, 145], [120, 165]] + const setB = [[90, 40], [150, 70], [45, 85], [135, 90], [55, 115], [145, 120], [100, 135], [160, 150]] + const extra = [[90, 175]] + const self = [105, 75] + const panel = (ox, label, sampled, unsampled, color, arrowTo) => { + let s = txt(ox, 66, label, { size: 13, fill: C.sub }) + s += `` + for (const [x, y] of [...unsampled, ...extra]) s += `` + for (const [x, y] of sampled) s += `` + s += `` + s += arrow(ox + self[0], 76 + self[1], ox + arrowTo[0], 76 + arrowTo[1], color, 2.5) + return s + } + b += panel(30, 'tick t — sample A (amber)', setA, setB, C.covered[2], [68, 32]) + b += txt(30, 300, 'force on ⊙ = Σ over sample A × m/K', { size: 12, fill: C.sub }) + b += panel(360, 'tick t+1 — sample B (teal)', setB, setA, C.covered[1], [144, 116]) + b += txt(360, 300, 'same positions, different draw → different answer', { size: 12, fill: C.sub }) + b += txt(690, 66, 'the point, tick after tick', { size: 13, fill: C.sub }) + b += `` + b += `` + b += `` + b += txt(690, 300, '≈0.5 units/tick, ~92° mean turn — at equilibrium', { size: 12, fill: C.sub }) + b += txt(20, 345, 'The estimate is unbiased — averaged over ticks it equals the exact force — but each tick draws a fresh K-subset and', { size: 13 }) + b += txt(20, 365, 'weights it by m/K. While the cell stays dense (m ≫ K) that per-tick difference never shrinks: the noise is not annealing,', { size: 13 }) + b += txt(20, 385, 'it is a permanent random walk stacked on a settled layout.', { size: 13 }) + fs.writeFileSync(path.join(OUT, 'f-resample-jitter.svg'), svgDoc(W, H, b)) +} + +// ---------------------------------------------------------------- Diagram G +// The two paths after the fix: exact all-pairs ≤ 4,096 points, grid + sampling above. +{ + const W = 980; const H = 330 + const MONO = 'ui-monospace, SFMono-Regular, Menlo, monospace' + let b = '' + b += txt(20, 34, 'Two paths — the point count picks one', { size: 16, weight: '600' }) + b += `` + b += txt(75, 152, 'graph of', { size: 13, anchor: 'middle' }) + b += txt(75, 170, 'n points', { size: 13, anchor: 'middle' }) + b += `` + b += `n ≤ 4096` + b += txt(240, 168, 'usesAllPairs', { size: 11, fill: C.sub, anchor: 'middle' }) + b += arrow(130, 156, 156, 156, C.text, 1.5) + b += arrow(266, 129, 352, 94, C.text, 1.5) + b += txt(300, 96, 'no', { size: 12, fill: C.sub }) + b += `` + b += txt(438, 84, 'aggregate the', { size: 13, anchor: 'middle' }) + b += txt(438, 102, 'grid pyramid', { size: 13, anchor: 'middle' }) + b += arrow(518, 88, 552, 88, C.text, 1.5) + b += `` + b += txt(636, 84, 'depth-peel K slots', { size: 13, anchor: 'middle' }) + b += txt(636, 102, 'per finest cell', { size: 13, anchor: 'middle' }) + b += `K = 32 (≤16k) · 16 (≤65k) · 8 above` + b += arrow(716, 88, 750, 88, C.text, 1.5) + b += `` + b += txt(854, 84, 'far field per level +', { size: 13, anchor: 'middle' }) + b += txt(854, 102, 'sampled near field ×m/K', { size: 13, anchor: 'middle' }) + b += arrow(266, 183, 424, 234, C.text, 1.5) + b += txt(318, 226, 'yes', { size: 12, fill: C.sub }) + b += `` + b += `force-allpairs.frag` + b += txt(555, 262, 'one exact O(n²) pass — no grid, no sampling', { size: 13, anchor: 'middle' }) + b += `` + b += txt(905, 222, 'velocity', { size: 13, anchor: 'middle' }) + b += txt(905, 240, 'texture', { size: 13, anchor: 'middle' }) + b += arrow(854, 116, 887, 196, C.text, 1.5) + b += txt(892, 156, 'additive', { size: 11, fill: C.sub }) + b += arrow(680, 242, 846, 230, C.text, 1.5) + b += txt(742, 222, 'single write', { size: 11, fill: C.sub }) + b += txt(20, 312, 'Same pairwise falloff and coincident-point kick on both paths — the small-graph path just computes every pair instead of estimating.', { size: 12, fill: C.sub }) + fs.writeFileSync(path.join(OUT, 'g-two-paths.svg'), svgDoc(W, H, b)) +} + +// ---------------------------------------------------------------- Diagram H +// Peeling ping-pong: two plain 2D targets, each pass's result copied into its +// layer of the slot array texture (rendering into a layer of a texture you are +// sampling another layer of is a WebGL feedback loop). +{ + const W = 980; const H = 430 + const MONO = 'ui-monospace, SFMono-Regular, Menlo, monospace' + let b = '' + b += txt(20, 34, 'Peeling ping-pong: two plain 2D targets, copied into array layers', { size: 16, weight: '600' }) + const cols = [ + { x: 100, label: 'pass 0', l1: 'writes target A', l2: 'smallest hash / cell', soft: C.coveredSoft[0] }, + { x: 320, label: 'pass 1', l1: 'reads A · writes B', l2: 'next-smallest hash', soft: C.coveredSoft[1] }, + { x: 540, label: 'pass 2', l1: 'reads B · writes A', l2: 'next-smallest hash', soft: C.coveredSoft[0] }, + { x: 760, label: 'pass 3 … K−1', l1: 'reads A · writes B', l2: '… alternating', soft: C.coveredSoft[1] }, + ] + cols.forEach((col, k) => { + b += txt(col.x + 80, 96, col.label, { size: 13, fill: C.sub, anchor: 'middle' }) + b += `` + b += txt(col.x + 80, 136, col.l1, { size: 13, anchor: 'middle' }) + b += txt(col.x + 80, 156, col.l2, { size: 12, fill: C.sub, anchor: 'middle' }) + if (k < cols.length - 1) { + b += arrow(col.x + 160, 142, col.x + 216, 142, C.text, 1.5) + b += txt(col.x + 188, 130, 'winner', { size: 11, fill: C.sub, anchor: 'middle' }) + } + b += arrow(col.x + 80, 174, col.x + 80, 266, C.covered[2], 2) + b += `` + b += txt(col.x + 80, 298, k < 3 ? `layer ${k}` : 'layer 3 … K−1', { size: 13, anchor: 'middle' }) + }) + b += `copyTextureToTexture` + b += `` + b += `slotsTexture : sampler2DArray — the force pass loops s = 0 … slotCount−1` + b += txt(20, 390, 'Why not render into layer k directly? Pass k must sample pass k−1\'s output — and sampling one layer of a texture while', { size: 13 }) + b += txt(20, 410, 'rendering to another layer of the same texture is a WebGL feedback loop. So passes render to plain 2D targets instead.', { size: 13 }) + fs.writeFileSync(path.join(OUT, 'h-pingpong-peel.svg'), svgDoc(W, H, b)) +} + console.log('generated:', fs.readdirSync(OUT).join(', ')) diff --git a/docs/many-body-force/h-pingpong-peel.svg b/docs/many-body-force/h-pingpong-peel.svg new file mode 100644 index 00000000..96fa3d58 --- /dev/null +++ b/docs/many-body-force/h-pingpong-peel.svg @@ -0,0 +1,9 @@ + + + + + + + +Peeling ping-pong: two plain 2D targets, copied into array layerspass 0writes target Asmallest hash / cellwinnerlayer 0pass 1reads A · writes Bnext-smallest hashwinnerlayer 1pass 2reads B · writes Anext-smallest hashwinnerlayer 2pass 3 … K−1reads A · writes B… alternatinglayer 3 … K−1copyTextureToTextureslotsTexture : sampler2DArray — the force pass loops s = 0 … slotCount−1Why not render into layer k directly? Pass k must sample pass k−1's output — and sampling one layer of a texture whilerendering to another layer of the same texture is a WebGL feedback loop. So passes render to plain 2D targets instead. + \ No newline at end of file diff --git a/docs/many-body-force/jitter-fix.html b/docs/many-body-force/jitter-fix.html new file mode 100644 index 00000000..75e1b21e --- /dev/null +++ b/docs/many-body-force/jitter-fix.html @@ -0,0 +1,448 @@ + + + + + +Near-Field Shimmer — the cosmos.gl many-body jitter fix + + + +
+
+

cosmos.gl · many-body force

+

Near-Field Shimmer

+

Why a settled layout wouldn't hold still, how the cause was pinned to + one line of sampling design, and how the fix split the repulsion force in two. For the full + algorithm walkthrough, see the many-body deep dive; this page + covers one problem and its solution.

+
+ +

The one-sentence version: the Monte-Carlo near field re-draws its random +sample every tick, which is harmless while dense clumps disperse but becomes permanent visible +shimmer when links or gravity keep a settled layout dense — so the force computes small graphs +(≤ 4,096 points) exactly with a dedicated all-pairs pass, and gives larger graphs +2–4× more sampling slots via an adaptive, sampler2DArray-backed slot count.

+ +
+

01 · Background

+

How the repulsion force computes

+

Repulsion is O(n²) if computed literally, so the engine uses a P3M + scheme: a grid pyramid lumps far-away mass into cell centroids, and the only region that + needs real pairwise forces — the 3×3 finest-cell neighborhood around each point — gets a + Monte-Carlo estimate of them.

+ +
+
+ + + + + + +One simulation tick of the repulsion force (all on the GPU)positionsone texel per point(x, y)1 · aggregate levelsdraw n points into each grid,additive blend accumulates[Σx, Σy, count] per cell2 · build near-field slotsK depth-peel passes over thefinest grid: a fresh randomK-subset per cell (K = 32/16/8)3 · force passesper level: centroid repulsion+ near field: weighted pairs;all add into the velocity textureThe integration step (velocity → positions) is shared with all other forces and unchanged.Levels: 4², 8², … up to ≈ 2·√n per axis (capped 512²). Slot array: K layers × finest grid, [point index, hash] each.Graphs of ≤ 4,096 points skip all three passes: a single exact all-pairs pass writes the velocity texture directly. +
+
One simulation tick, entirely on the GPU: aggregate the pyramid, peel the + near-field sample, then blend every force pass additively into the velocity texture.
+
+ +
+
+ + + + + + +Level 0: 4 × 4 cellswhole space minus the 3×3 shellcentroid force at this level3×3 → deferred to the next levelLevel 1: 8 × 8 cellsprevious 3×3 refined, minus its own 3×3centroid force at this level3×3 → deferred to the next levelLevel 2: 16 × 16 cellsprevious 3×3 refined, minus its own 3×3centroid force at this level3×3 → near field (Monte-Carlo pairs)Every region of space is charged to exactly one pass — no gaps, no double counting. The finest grid adapts to n (≈ 2·√n cells per axis, 8²…512²). +
+
The pyramid tiles space exactly once: each level covers its 6×6 child block + minus its own 3×3, and the leftover finest 3×3 neighborhood is the near field.
+
+ +

The near field is the part this page is about. Every tick: each point gets a fresh + pseudo-random hash; K depth-peeling passes select, per finest cell, the K + points with the smallest hashes — a uniform random K-subset, re-drawn from scratch every + tick; then the force pass sums true pairwise forces from those K points and scales the sum + by others / sampled — the Horvitz–Thompson estimator:

+ +

src/modules/ForceManyBody/force-nearfield.frag

+ 
vec2 pairSum = vec2(0.0);
+float sampled = 0.0;
+for (int s = 0; s < slots; s += 1) {
+  vec2 slot = texelFetch(slotsTexture, ivec3(cell, s), 0).rg;
+  if (slot.x < 0.0) break;                    // cell exhausted
+  pairSum += slotVelocity(slot, position, selfIndex, random.rg, sampled);
+}
+// Each of the cell's `others` points had equal probability sampled/others
+// of being drawn; dividing by it makes E[estimate] = the exact all-pairs sum.
+if (sampled > 0.0) velocity += (others / sampled) * pairSum;
+ +
+
+ + + + + + +One finest-level cell, one tick13 points, each hashed with this tick’s random seed.62.77.56.97.91.81.94.74.64.20.93.12.82sampled (the K smallest hashes)left outK depth-peeling passes → K slot array layers (K = 8 drawn)pass k keeps the smallest hash not yet peeledslot 0point #11hash 0.12slot 1point #9hash 0.20slot 2point #2hash 0.56slot 3point #0hash 0.62slot 4point #8hash 0.64slot 5point #7hash 0.74slot 6point #1hash 0.77slot 7point #5hash 0.81Horvitz–Thompson weightingsampled sum = full sum, on averagecell has 12 other points,8 of them sampled →F ≈ (12 / 8) · Σ F(sampled pair)E[F] = exact all-pairs sum(unbiased, no centroid term)≤ K points in cell → exact.A fresh random subset every tick: thenoise anneals with alpha while densitydisperses. Sustained density is why Kadapts (32/16/8) and ≤ 4k graphs go exact. +
+
Depth peeling extracts each cell's K smallest-hash points; the + Horvitz–Thompson weight makes the sampled sum an unbiased estimate of the exact one.
+
+ +

Two properties matter here:

+
    +
  • A cell with ≤ K points is sampled exhaustively (others == sampled, + weight 1) — its near field is exact. With the finest grid at ~2·√n cells per axis the + average cell holds ¼ point, so for typical graphs the whole near field is exact.
    • +
  • A cell with more than K points gets an unbiased but noisy + estimate — and the noise is re-rolled every tick, because the sample is.
    • +
+
+ +
+

02 · The problem

+

Sampling noise that never anneals

+

The design treated the per-tick noise as a feature — annealing jitter. Every force is + scaled by alpha, so the noise is large early, when clumps genuinely need random + kicks to break apart, and fades as the simulation cools. Under pure repulsion + the argument is airtight twice over: the dense clump that causes the variance + disperses within ~60 ticks, occupancy falls toward K, and the noise dies with it — sub-pixel + before anyone sees it.

+

The hole: the argument assumes density disperses. When link attraction or gravity + holds a hub together while alpha stays high — a long-running layout, a reheated one, + start() on interaction — cell occupancy stays far above K forever, and + the re-drawn sample turns into permanent, visible shimmer.

+ +
+
+ + + + + + +Same cell, same positions — a fresh sample every ticktick t — sample A (amber)force on ⊙ = Σ over sample A × m/Ktick t+1 — sample B (teal)same positions, different draw → different answerthe point, tick after tick≈0.5 units/tick, ~92° mean turn — at equilibriumThe estimate is unbiased — averaged over ticks it equals the exact force — but each tick draws a fresh K-subset andweights it by m/K. While the cell stays dense (m ≫ K) that per-tick difference never shrinks: the noise is not annealing,it is a permanent random walk stacked on a settled layout. +
+
The mechanism of the bug: with the cell holding m ≈ 45 others and K = 8 slots, + each tick multiplies a different random 8 by 45/8 ≈ 5.6. Unbiased in expectation — + but at equilibrium the tick-to-tick swing is pure noise, annealing nothing.
+
+ +

The graph that surfaced it

+

The 163-country border-adjacency network (642 links) — small, real, and exactly the + sustained-density shape: link attraction plus gravity pull it into a clump a few finest cells + wide, occupancy ~45 per cell against 8 slots. Measured at equilibrium with alpha held at 1, + every point wandered ~0.5 units per tick with a ~92° mean direction + change between consecutive steps. A ~90° mean turn is the signature of a pure random + walk: consecutive steps are uncorrelated. A settled layout should be still; this one shimmered + indefinitely.

+ +

How the cause was pinned

+

A CPU replica of the force differing only in the per-tick K = 8 re-sampling + reproduced the GPU numbers digit-for-digit, while a CPU exact all-pairs reference + under identical falloff, friction, and alpha was ~1000× stiller. That rules out integration + and damping — the re-sampling is the jitter.

+ +

A settled layout is still. Sampling noise may only exist + where it is annealing something, never at equilibrium of a small or medium graph.

+ +

The trajectory, measured

+

The Performance → Repulsion Jitter: Fixed vs Before Storybook story runs + the country graph twice with identical data and seed — once on today's exact path, once with + the pre-fix sampled configuration (K = 8) forced back on — and traces one point from the + densest finest cell. These are its trajectory panels, captured live, same point in both + configurations:

+ +
+
+ Fixed: trajectory is a single smooth arc spanning about 3 units over 360 ticks +

Fixed (exact all-pairs path): step 0.02 u/tick · + turn 0.6° — the arc is the layout's slow coherent drift under alpha held + at 1, not noise.

+
+
+ Before: trajectory is a dense random-walk tangle spanning 6.76 units over 360 ticks +

Before the fix (sampled near field, K = 8): + step 0.46 u/tick · turn 94° — a random walk stacked on a settled + layout.

+
+
+ +

The shapes say more than the spans: smooth flow versus a scribble whose direction re-rolls + every tick. That's why the turn angle is the story's key metric — it separates flow (0.6°) + from noise (94°) even when both cover ground.

+
+ +
+

03 · The fix, part one

+

Small graphs are exact

+

src/modules/ForceManyBody/index.ts

+ 
const ALL_PAIRS_MAX_POINTS = 4096
+
+/** Small graphs skip the grid + Monte-Carlo machinery entirely. */
+private get usesAllPairs (): boolean {
+  return (this.data.pointsNumber ?? 0) <= ALL_PAIRS_MAX_POINTS
+}
+
+public run (): void {
+  // …
+  if (this.usesAllPairs) {
+    this.drawAllPairsForce()   // one pass, done
+    return
+  }
+  this.drawLevels()            // grid pyramid
+  this.drawNearFieldSlots()    // K peeling passes
+  this.drawForces()            // far field + sampled near field
+}
+ +
+
+ + + + + + +Two paths — the point count picks onegraph ofn pointsn ≤ 4096usesAllPairsnoaggregate thegrid pyramiddepth-peel K slotsper finest cellK = 32 (≤16k) · 16 (≤65k) · 8 abovefar field per level +sampled near field ×m/Kyesforce-allpairs.fragone exact O(n²) pass — no grid, no samplingvelocitytextureadditivesingle writeSame pairwise falloff and coincident-point kick on both paths — the small-graph path just computes every pair instead of estimating. +
+
The point count picks the path. Both write the same velocity texture with the + same pairwise falloff — the small-graph path just computes every pair instead of estimating.
+
+ +

At or below 4,096 points the whole force becomes one full-screen pass: + each fragment is one point, looping over every other point.

+ +

src/modules/ForceManyBody/force-allpairs.frag

+ 
for (int i = 0; i < count; i += 1) {
+  if (i == selfIndex) continue;
+  ivec2 texel = ivec2(i % size, i / size);
+  if (texelFetch(exitTexture, texel, 0).g > 0.5) continue;   // absent point
+  vec2 otherPosition = texelFetch(positionsTexture, texel, 0).rg;
+  // … same clamped inverse-distance falloff as the grid path,
+  //   split into near/far sums at the near-field radius …
+}
+ +

The pairwise falloff is identical to the grid path's (the shader comment + pins it: "must stay identical"), including the coincident-point random kick — so the physics + is unchanged, only computed exactly. No grid, no peeling, no sampling → exact at any + cell occupancy → zero noise, nothing to shimmer. And the two paths bound the same thing: + pairs within the grid path's near-field scale are jittered and capped at 2 × the finest cell + size it would use (so a coincident stack expands instead of teleporting), while farther + pairs pass through unbounded, as the level passes do — dynamics stay continuous across the + threshold.

+ +

Why exact is also faster here

+

Counterintuitive but measured: depth peeling is inherently sequential — one render pass per + slot, ~0.1 ms of fixed cost each — while the O(n²) texel loop is trivial work at this scale + (4096² ≈ 17M pair evaluations, around a millisecond on modest GPUs). At 2k points, a 64-slot + peel measured ~6.4 ms/step; the single all-pairs pass measures + ~1.8 ms/step. The 4,096 threshold sits where this stops being true — the next + power of two would already cost several milliseconds of n² work.

+
+ +
+

04 · The fix, part two

+

Adaptive K above the threshold

+

Above 4,096 points, n² is genuinely unaffordable (100k² = 10¹⁰ pairs), so the sampler + stays — but the noise amplitude scales like occupancy/K · 1/√K, so raising K + buys a lot: going 8 → 32 slots cuts amplitude ~8× and widens the exactly-covered + range to cells of ≤ 32 points.

+ +

src/modules/ForceManyBody/index.ts

+ 
const getNearFieldSlotCount = (pointsNumber: number): number => {
+  if (pointsNumber <= 16384) return 32
+  if (pointsNumber <= 65536) return 16
+  return 8   // the ≥65k tier keeps exactly its old cost
+}
+ +

Peel cost is K sequential passes, so K is a direct cost/variance dial: mid-size graphs — + where a single dense hub is still common and per-point noise still visible — buy more slots; + huge graphs keep the cheap 8-slot estimator, where per-point noise is sub-pixel anyway and the + peel cost would dominate.

+ +

The enabler: sampler2DArray

+

Previously K was hard-wired at 8 as eight hand-unrolled samplers, because + GLSL ES 3.0 can't index a sampler2D[] with a loop variable:

+ +

before — force-nearfield.frag

+ 
uniform sampler2D slotTexture0;
+uniform sampler2D slotTexture1;
+// … slotTexture2 … slotTexture7, and eight matching unrolled reads.
+// Changing K meant editing the shader (and the bindings, and the constant).
+ +

after — one array texture, K is a runtime uniform

+ 
uniform highp sampler2DArray slotsTexture;
+…
+int slots = int(slotCount);
+for (int s = 0; s < slots; s += 1) {
+  vec2 slot = texelFetch(slotsTexture, ivec3(cell, s), 0).rg;
+  if (slot.x < 0.0) break;   // exhausted cell peels empty (−1) slots
+  pairSum += slotVelocity(slot, position, selfIndex, random.rg, sampled);
+}
+ +

That swap is what makes K a plain runtime value instead of a shader edit — and it lifts the + old practical ceiling of the WebGL2 texture-unit budget (8 slot samplers plus the other bound + textures was pressing against the 16-unit floor).

+ +

The ping-pong wrinkle

+

Peeling pass k must sample pass k−1's output (to know which + points are already taken), and sampling one layer of a texture while rendering to another + layer of the same texture is a WebGL feedback loop. So the passes ping-pong between + two plain 2D targets, and each pass's result is copied into its array layer:

+ +
+
+ + + + + + +Peeling ping-pong: two plain 2D targets, copied into array layerspass 0writes target Asmallest hash / cellwinnerlayer 0pass 1reads A · writes Bnext-smallest hashwinnerlayer 1pass 2reads B · writes Anext-smallest hashwinnerlayer 2pass 3 … K−1reads A · writes B… alternatinglayer 3 … K−1copyTextureToTextureslotsTexture : sampler2DArray — the force pass loops s = 0 … slotCount−1Why not render into layer k directly? Pass k must sample pass k−1's output — and sampling one layer of a texture whilerendering to another layer of the same texture is a WebGL feedback loop. So passes render to plain 2D targets instead. +
+
Two render targets alternate as write/read; every pass's winners are copied + into their own layer of the array texture the force pass loops over.
+
+
+ +
+

05 · Results

+

Before and after

+
+ + + + + + + +
MeasurementBeforeAfter
Country graph, mean step at equilibrium~0.46 units/tick0.02 units/tick
Country graph, mean turn angle~92° (random walk)0.6° (still)
Synthetic 1,024-point single-cell clump vs CPU exact referencenoisydirection noise 0.0°, path efficiency 1.000 — indistinguishable, incl. a gravity-confined 260-points/cell-forever case
2k points~6.4 ms/step (64-slot experiment)1.81 ms/step (exact)
≥ 100k pointsbaselineunchanged — 100k 6.63, 200k 13.79 ms/step
+
+

Full benchmark row after the fix: 2k 1.81 · 5k 3.80 · 20k 2.29 · 50k 3.90 · 100k 6.63 · + 200k 13.79 ms/step. The 5k reading sitting above 20k is the adaptive K at work: 5k is just + over the all-pairs threshold, so it pays 32 sequential peel passes; 20k runs K = 16.

+

Nothing changed at the API surface: no new config, and simulationRepulsion + behaves as before. See it live in Performance → Repulsion Jitter: Fixed vs + Before, and the full algorithm in the deep dive.

+
+
+ + diff --git a/history/2026/2026-08-14-nearfield-jitter.md b/history/2026/2026-08-14-nearfield-jitter.md new file mode 100644 index 00000000..8ae08404 --- /dev/null +++ b/history/2026/2026-08-14-nearfield-jitter.md @@ -0,0 +1,78 @@ + + +# Near-field sampling jitter: exact small-graph path + adaptive slot count + +**Date:** 2026-08-14 +**Commits:** `feat(stories): country-borders story reproducing near-field sampling jitter` (`776d15a`), `fix(force): exact all-pairs repulsion below 4k points; adaptive near-field sampling above` (`ea88026`), `feat(stories): replace the jitter repro with a fixed-vs-before comparison — both paths visible at HEAD` (`2fe05a6`) + +## Why + +The Monte-Carlo near field (`2026-07-08-many-body-repulsion.md`) re-draws each finest cell's +K = 8 sample every tick and weights it by `count/sampled`. The design treated the resulting +per-tick variance as annealing jitter — correct for pure repulsion, where the dense clump that +causes the variance disperses within ~60 ticks and the noise dies with it, sub-pixel before +anyone sees it. But when link attraction or gravity holds density up *while alpha stays high* +(long-running or reheated layouts), cell occupancy never falls toward K and the re-sampling +noise becomes permanent visible shimmer. Surfaced on a real 163-country border graph: at +equilibrium every point wandered ~0.5 units/tick with ~92° mean direction change between +consecutive ticks — a pure random walk on top of a settled layout. A CPU all-pairs reference +under identical falloff/friction/alpha was ~1000× stiller, and a CPU replica differing *only* +in the per-tick K = 8 re-sampling reproduced the GPU numbers digit-for-digit — pinning the +cause to the sampling, not integration or damping. + +## What changed + +Two complementary mechanisms in `src/modules/ForceManyBody/`: + +- **Exact all-pairs path for small graphs** (`force-allpairs.frag`): at + `pointsNumber ≤ ALL_PAIRS_MAX_POINTS` (4,096) the whole force is one O(n²) full-screen pass — + same clamped inverse-distance falloff and coincident-point kick as the grid path, absent + (NaN) points skipped, no pyramid or slot allocation at all. Exact at any occupancy, so zero + sampling noise. Also *faster* there: depth peeling is one sequential render pass per slot + (~0.1 ms fixed cost each), while the n² texel loop is trivial at this scale — measured + ~1.8 ms/step at 2k points vs ~6.4 ms for a 64-slot peel. (#240 prototyped and dropped this + path when the grid looked "effectively exact" for small graphs; sustained-density layouts are + why it returned.) +- **Adaptive slot count via `sampler2DArray`**: the near-field slots moved from 8 hand-unrolled + `sampler2D`s (WebGL2 texture-unit ceiling, and unrollable only by editing the shader) to one + array texture looped with a `slotCount` uniform. `getNearFieldSlotCount` now returns + 32 (≤ 16k points) / 16 (≤ 65k) / 8 (above) — more slots extend the exactly-covered occupancy + range and shrink residual variance (amplitude ∝ occupancy/K · 1/√K), while the ≥ 65k tier is + cost-identical to before. Peeling ping-pongs between two plain 2D targets and copies each + pass's result into its array layer — pass k must sample pass k−1's output, and sampling one + layer of a texture while rendering to another layer of the same texture is a WebGL feedback + loop. + +## Results + +- Country borders graph (163 points, 642 links, alpha held at 1): step 0.46 → 0.02 units/tick, + mean turn 92.5° → 0.8° — visually still. +- A synthetic dense clump measured during the investigation (1,024 points in one finest cell): + direction noise 0.0° in every tick window, path efficiency 1.000 — indistinguishable from + the CPU exact reference, including the pathological gravity-confined case (260 points/cell + forever). +- Repulsion benchmark: 2k **1.81 ms/step** (exact), 5k 3.80, 20k 2.29, 50k 3.90, + 100k 6.63, 200k 13.79 — the ≥ 100k path byte-identical in cost to before the change. + +## Notes + +- No config or public-API change; `simulationRepulsion` behaves as before. Not a breaking + change — `migration-notes.md` intentionally untouched. +- The deep dive (`docs/many-body-force/README.md`) now documents the two-path structure, the + sustained-density failure mode, and the adaptive K. + +## Example + +- **Repulsion Jitter: Fixed vs Before** (`src/stories/performance/country-borders-comparison.ts`, + Storybook *Performance*): the real graph that surfaced the bug, run twice side by side with + identical data and seed — left today's exact path (settled and still), right the pre-fix + configuration (sampled near field, K = 8) forced back on through a story-only patch of + ForceManyBody internals via the repo's src alias (per-instance config marker; fails loudly if + the internals move; does not resolve against the published package). Each side has a + step/turn meter and a trajectory panel tracing one dense-cell point over 360 ticks — a + smooth drift arc today vs a random-walk tangle before. +- This story replaced the original single-pane repro + (`feat(stories): country-borders story reproducing near-field sampling jitter`): after the + fix, that story ran on the exact path and could no longer show live the shimmer it + documented. It was committed *before* the fix precisely so the shimmer is observable at its + own commit — that bisectable evidence remains in git history. diff --git a/src/modules/ForceManyBody/force-allpairs.frag b/src/modules/ForceManyBody/force-allpairs.frag new file mode 100644 index 00000000..aea36edc --- /dev/null +++ b/src/modules/ForceManyBody/force-allpairs.frag @@ -0,0 +1,115 @@ +#version 300 es +precision highp float; +// Fragment shaders default int to mediump, guaranteed only to 32767 — +// point indices go far higher. +precision highp int; + +// Exact all-pairs repulsion for small graphs. One fragment per point, looping +// over every other point — O(n²) total, but below the brute-force threshold a +// single pass is both cheaper than the grid pyramid's sequential depth-peeling +// passes and exact at any cell occupancy: no Monte-Carlo sampling, hence none +// of the per-tick re-sampling noise that shows up as shimmer in dense layouts +// (see force-nearfield.frag for the sampled path used above the threshold). + +uniform sampler2D positionsTexture; +uniform sampler2D randomValues; +uniform sampler2D exitTexture; + +#ifdef USE_UNIFORM_BUFFERS +layout(std140) uniform forceAllPairsUniforms { + float pointsTextureSize; + float pointsNumber; + float alpha; + float repulsion; + float maxStep; +} forceAllPairs; + +#define pointsTextureSize forceAllPairs.pointsTextureSize +#define pointsNumber forceAllPairs.pointsNumber +#define alpha forceAllPairs.alpha +#define repulsion forceAllPairs.repulsion +#define maxStep forceAllPairs.maxStep +#else +uniform float pointsTextureSize; +uniform float pointsNumber; +uniform float alpha; +uniform float repulsion; +uniform float maxStep; +#endif + +out vec4 fragColor; + +// Same clamped inverse-distance falloff as the grid-path shaders (must stay identical). +vec2 pairwiseVelocity(vec2 position, vec2 otherPosition, vec2 randomDir) { + vec2 distVector = position - otherPosition; + float l = dot(distVector, distVector); + if (l <= 0.0) { + // Exactly coincident points have no separation direction, so an + // inverse-distance force is undefined and they would stay stacked forever. + // Kick along this point's random vector instead (each point has a + // different one, so a pile disperses). + distVector = randomDir; + l = dot(distVector, distVector); + if (l <= 0.0) return vec2(0.0); + } + float distanceMin2 = 1.0; + if (l < distanceMin2) l = sqrt(distanceMin2 * l); + float addV = alpha * repulsion / sqrt(l); + return addV * normalize(distVector); +} + +void main() { + ivec2 pointTexel = ivec2(gl_FragCoord.xy); + int size = int(pointsTextureSize); + int selfIndex = pointTexel.y * size + pointTexel.x; + int count = int(pointsNumber); + + // Fragments beyond the point count are unused texture pixels. + if (selfIndex >= count) { + fragColor = vec4(0.0); + return; + } + + // An absent point must neither move nor repel (its position is NaN). + vec4 selfExit = texelFetch(exitTexture, pointTexel, 0); + if (selfExit.g > 0.5) { + fragColor = vec4(0.0); + return; + } + + vec2 position = texelFetch(positionsTexture, pointTexel, 0).rg; + vec4 random = texelFetch(randomValues, pointTexel, 0); + + // Pairs are split at the grid path's near-field scale (maxStep = 2 × the + // finest cell size it would use at this point count): closer pairs + // correspond to its 3×3 near-field pass — jittered and bounded below, the + // way that pass bounds its own sum — farther pairs to its level passes, + // which it leaves unbounded. The same split keeps the dynamics continuous + // across the point-count threshold between the two paths. + vec2 nearVelocity = vec2(0.0); + vec2 farVelocity = vec2(0.0); + float nearRadius2 = maxStep * maxStep; + for (int i = 0; i < count; i += 1) { + if (i == selfIndex) continue; + ivec2 texel = ivec2(i % size, i / size); + if (texelFetch(exitTexture, texel, 0).g > 0.5) continue; + vec2 otherPosition = texelFetch(positionsTexture, texel, 0).rg; + vec2 distVector = position - otherPosition; + vec2 pairVelocity = pairwiseVelocity(position, otherPosition, random.rg); + if (dot(distVector, distVector) < nearRadius2) nearVelocity += pairVelocity; + else farVelocity += pairVelocity; + } + + // Random jitter proportional to the near velocity, to keep points from + // sticking (same as the near-field pass). + nearVelocity += nearVelocity * random.rg; + + // Bound the per-tick near kick. Exactness alone does not bound it: the + // falloff still diverges at near-zero separations, and a coincident stack + // sums n−1 same-direction random kicks into one fling. The clamp caps the + // magnitude and keeps the direction. + float speed = length(nearVelocity); + if (speed > maxStep) nearVelocity *= maxStep / speed; + + fragColor = vec4(farVelocity + nearVelocity, 0.0, 0.0); +} diff --git a/src/modules/ForceManyBody/force-nearfield.frag b/src/modules/ForceManyBody/force-nearfield.frag index 73d62d4d..ed3073da 100644 --- a/src/modules/ForceManyBody/force-nearfield.frag +++ b/src/modules/ForceManyBody/force-nearfield.frag @@ -24,17 +24,10 @@ precision highp int; uniform sampler2D positionsTexture; uniform sampler2D levelTexture; uniform sampler2D randomValues; -// One sampler per near-field slot. We list them out instead of using an array -// because WebGL2's GLSL won't let you index a sampler array in a loop. Keep this -// list the same length as NEAR_FIELD_SLOTS in index.ts. -uniform sampler2D slotTexture0; -uniform sampler2D slotTexture1; -uniform sampler2D slotTexture2; -uniform sampler2D slotTexture3; -uniform sampler2D slotTexture4; -uniform sampler2D slotTexture5; -uniform sampler2D slotTexture6; -uniform sampler2D slotTexture7; +// All near-field slots live in one array texture (one layer per depth-peeling +// pass), so the slot count is a runtime uniform instead of a hard-wired list +// of sampler2Ds. Float data — highp, the default sampler precision is lowp. +uniform highp sampler2DArray slotsTexture; #ifdef USE_UNIFORM_BUFFERS layout(std140) uniform forceNearFieldUniforms { @@ -43,6 +36,7 @@ layout(std140) uniform forceNearFieldUniforms { float cellSize; float alpha; float repulsion; + float slotCount; } forceNearField; #define pointsTextureSize forceNearField.pointsTextureSize @@ -50,12 +44,14 @@ layout(std140) uniform forceNearFieldUniforms { #define cellSize forceNearField.cellSize #define alpha forceNearField.alpha #define repulsion forceNearField.repulsion +#define slotCount forceNearField.slotCount #else uniform float pointsTextureSize; uniform float levelGridSize; uniform float cellSize; uniform float alpha; uniform float repulsion; +uniform float slotCount; #endif out vec4 fragColor; @@ -126,18 +122,15 @@ void main() { vec2 pairSum = vec2(0.0); float sampled = 0.0; - // Same story as the sampler list above: no looping over samplers in - // WebGL2, so we read each slot on its own line. This has to match - // NEAR_FIELD_SLOTS too (and the samplers above, and the bindings in - // index.ts). - pairSum += slotVelocity(texelFetch(slotTexture0, cell, 0).rg, position, selfIndex, random.rg, sampled); - pairSum += slotVelocity(texelFetch(slotTexture1, cell, 0).rg, position, selfIndex, random.rg, sampled); - pairSum += slotVelocity(texelFetch(slotTexture2, cell, 0).rg, position, selfIndex, random.rg, sampled); - pairSum += slotVelocity(texelFetch(slotTexture3, cell, 0).rg, position, selfIndex, random.rg, sampled); - pairSum += slotVelocity(texelFetch(slotTexture4, cell, 0).rg, position, selfIndex, random.rg, sampled); - pairSum += slotVelocity(texelFetch(slotTexture5, cell, 0).rg, position, selfIndex, random.rg, sampled); - pairSum += slotVelocity(texelFetch(slotTexture6, cell, 0).rg, position, selfIndex, random.rg, sampled); - pairSum += slotVelocity(texelFetch(slotTexture7, cell, 0).rg, position, selfIndex, random.rg, sampled); + // One layer per depth-peeling pass. An exhausted cell peels empty slots + // (index -1) for the remaining layers; the early break skips them — + // empty layers can't be followed by occupied ones within a cell. + int slots = int(slotCount); + for (int s = 0; s < slots; s += 1) { + vec2 slot = texelFetch(slotsTexture, ivec3(cell, s), 0).rg; + if (slot.x < 0.0) break; + pairSum += slotVelocity(slot, position, selfIndex, random.rg, sampled); + } // Horvitz–Thompson weighting: the sample is uniform among the cell's // other points (conditioned on whether the point itself was peeled), diff --git a/src/modules/ForceManyBody/index.ts b/src/modules/ForceManyBody/index.ts index c3af6ab6..b2d0b7e6 100644 --- a/src/modules/ForceManyBody/index.ts +++ b/src/modules/ForceManyBody/index.ts @@ -6,6 +6,7 @@ import calculateLevelFrag from '@/graph/modules/ForceManyBody/calculate-level.fr import calculateLevelPreciseVert from '@/graph/modules/ForceManyBody/calculate-level.vert?raw' import forceLevelPreciseFrag from '@/graph/modules/ForceManyBody/force-level.frag?raw' import forceNearFieldFrag from '@/graph/modules/ForceManyBody/force-nearfield.frag?raw' +import forceAllPairsFrag from '@/graph/modules/ForceManyBody/force-allpairs.frag?raw' import buildNearFieldSlotsVert from '@/graph/modules/ForceManyBody/build-nearfield-slots.vert?raw' import buildNearFieldSlotsFrag from '@/graph/modules/ForceManyBody/build-nearfield-slots.frag?raw' import { createIndexesForBuffer } from '@/graph/modules/Shared/buffer' @@ -18,21 +19,60 @@ import updateVert from '@/graph/modules/Shared/quad.vert?raw' */ const MAX_GRID_SIZE = 512 +/** + * Finest grid resolution per axis for a point count: ~2·√n, floored at 8², + * capped at MAX_GRID_SIZE. Shared by the pyramid allocation and the all-pairs + * pass's per-tick velocity clamp, which must bound with the same cell size the + * grid path would use at the same count. + */ +const getFinestGridSize = (pointsNumber: number): number => + Math.min(MAX_GRID_SIZE, Math.max(8, Math.pow(2, Math.ceil(Math.log2(2 * Math.sqrt(pointsNumber)))))) + /** * How many points per finest-level cell get exact pairwise repulsion each tick. - * We pick a fresh random subset every tick, so over time every point in a busy - * cell takes its turn being treated exactly. + * A cell holding at most this many points is sampled exhaustively — its near + * field is exact. Above it, a fresh random subset is drawn every tick and + * Horvitz–Thompson weighted; unbiased, but the per-tick re-drawing makes the + * force estimate noisy in proportion to occupancy/slots. In layouts where + * something keeps density up (link attraction into hubs, gravity) while alpha + * stays high, that noise is visible as per-point shimmer. + * + * So the slot count scales down as the graph grows: small graphs get enough + * slots that realistic cell occupancies are covered exactly (the country-scale + * graph that surfaced the jitter peaks around ~50 points per cell), while large + * graphs keep the cheap 8-slot estimator — at that scale per-point noise is + * sub-pixel and the peel cost (slots × points per tick) dominates instead. + * The slots live in one sampler2DArray layer each, so this is a plain runtime + * value — no shader changes needed to retune it. + * + * Graphs at or below ALL_PAIRS_MAX_POINTS never reach this path at all — they + * take the exact all-pairs pass instead, so the tiers start above it. + */ +const getNearFieldSlotCount = (pointsNumber: number): number => { + if (pointsNumber <= 16384) return 32 + if (pointsNumber <= 65536) return 16 + return 8 +} + +/** + * At or below this point count the whole force is computed exactly: one + * all-pairs O(n²) pass (force-allpairs.frag) replaces the grid pyramid and the + * Monte-Carlo near field. Two reasons it wins there: * - * Heads up — this isn't a knob you can just turn. The JS side (allocation, the - * peel loop) follows this number automatically, but the shader can't: WebGL2's - * GLSL won't let you loop over a sampler array, so force-nearfield.frag spells - * out every slot by hand. So if you change this, you also have to update - * force-nearfield.frag (the sampler list and the unrolled reads) and the - * bindings in drawForces — each of those spots has a matching note. If you'd - * rather make it truly tunable, switch the slots to a sampler2DArray and all the - * hand-syncing goes away. + * - Zero sampling noise at any cell occupancy. The sampled near field is only + * exact while a cell holds ≤ slot-count points; a small dense graph (hubs + * held tight by links or gravity) can concentrate hundreds of points in one + * finest cell, and the per-tick re-sampled estimate then jitters visibly. + * - It's cheaper. Depth peeling is inherently sequential — one render pass per + * slot — and at small point counts that fixed per-pass cost dominates the + * actual work (measured ~6ms/step for 64 slots at 2k points, vs ~1ms for the + * single all-pairs pass whose n² texel loop is trivial at this scale). + * + * The crossover is set by the n² fragment work: 4096² ≈ 17M pair evaluations + * per step stays around a millisecond on modest GPUs, while the next power of + * two would already cost several. */ -const NEAR_FIELD_SLOTS = 8 +const ALL_PAIRS_MAX_POINTS = 4096 /** A grid-level aggregation target ([sum(x), sum(y), count, 0] per cell). */ type LevelTarget = { @@ -42,15 +82,21 @@ type LevelTarget = { gridSize: number; } -/** A near-field depth-peeling slot target ([point index, hash] per cell). */ +/** A near-field depth-peeling render target ([point index, hash] per cell). */ type SlotTarget = { texture: Texture; fbo: Framebuffer; } +/** Ping-pong pair: each peel pass writes one and reads the other. */ +const PEEL_TARGETS = 2 + /** * GPU many-body (repulsion) force. * + * Graphs at or below ALL_PAIRS_MAX_POINTS are computed exactly in a single + * all-pairs pass (see that constant for why). Above it: + * * A Barnes-Hut-style grid pyramid (each level covers its aligned 6×6 child block * minus the Chebyshev-1 shell) whose finest 3×3 neighborhood is closed by an * unbiased Monte-Carlo near field: a per-tick depth-peeled random subset of each @@ -67,17 +113,29 @@ export class ForceManyBody extends CoreModule { /** Grid level count; `0` until create() allocates the pyramid. */ private levels = 0 private levelTargets = new Map() + /** Near-field slot count for the current point count (getNearFieldSlotCount). */ + private nearFieldSlots = 0 /** - * Near-field point slots: NEAR_FIELD_SLOTS textures sharing the finest - * level's grid layout, each holding [point index, hash] per cell — built by - * depth peeling every tick (see build-nearfield-slots.vert). + * Near-field point slots: one sampler2DArray layer per depth-peeling pass, + * sharing the finest level's grid layout, each holding [point index, hash] + * per cell — rebuilt every tick (see build-nearfield-slots.vert). */ - private nearFieldSlotTargets: SlotTarget[] = [] + private slotsArrayTexture: Texture | undefined + /** + * The two ping-pong peel render targets: pass k draws into k % 2 while + * reading the previous pass's result from (k + 1) % 2, then the result is + * copied into layer k of slotsArrayTexture. Peeling can't render into the + * array layers directly — pass k needs to sample pass k−1's output, and + * sampling one layer of a texture while rendering to another is a WebGL + * feedback loop. + */ + private peelTargets: SlotTarget[] = [] private calculateLevelsCommand: Model | undefined private forceLevelCommand: Model | undefined private buildNearFieldSlotsCommand: Model | undefined private forceNearFieldCommand: Model | undefined + private forceAllPairsCommand: Model | undefined private forceVertexCoordBuffer: Buffer | undefined @@ -115,19 +173,42 @@ export class ForceManyBody extends CoreModule { cellSize: number; alpha: number; repulsion: number; + slotCount: number; + }; + }> | undefined + + private forceAllPairsUniformStore: UniformStore<{ + forceAllPairsUniforms: { + pointsTextureSize: number; + pointsNumber: number; + alpha: number; + repulsion: number; + maxStep: number; }; }> | undefined private previousPointsTextureSize: number | undefined private previousPointsNumber: number | undefined + /** Small graphs skip the grid + Monte-Carlo machinery entirely (see ALL_PAIRS_MAX_POINTS). */ + private get usesAllPairs (): boolean { + return (this.data.pointsNumber ?? 0) <= ALL_PAIRS_MAX_POINTS + } + public create (): void { const { device, store } = this if (!store.pointsTextureSize) return // (Re)allocate the grid pyramid + near-field slots for the current point // count (resizing levels and dropping any that the pyramid no longer needs). - this.createLevels() + // Small graphs take the exact all-pairs pass and don't need any of it — + // drop whatever a previously larger graph left behind. + if (this.usesAllPairs) { + this.destroyLevelTargets() + this.levels = 0 + } else { + this.createLevels() + } // Random jitter texture to prevent sticking const totalPixels = store.pointsTextureSize * store.pointsTextureSize @@ -351,6 +432,7 @@ export class ForceManyBody extends CoreModule { cellSize: 'f32', alpha: 'f32', repulsion: 'f32', + slotCount: 'f32', }, defaultUniforms: { pointsTextureSize: store.pointsTextureSize, @@ -358,10 +440,57 @@ export class ForceManyBody extends CoreModule { cellSize: 0, alpha: store.alpha, repulsion: this.config.simulationRepulsion, + slotCount: 0, }, }, }) + // Exact all-pairs command (fullscreen quad — the small-graph path) + this.forceAllPairsUniformStore ||= new UniformStore(device, { + forceAllPairsUniforms: { + uniformTypes: { + // Order MUST match shader declaration order (std140 layout) + pointsTextureSize: 'f32', + pointsNumber: 'f32', + alpha: 'f32', + repulsion: 'f32', + maxStep: 'f32', + }, + defaultUniforms: { + pointsTextureSize: store.pointsTextureSize, + pointsNumber: data.pointsNumber, + alpha: store.alpha, + repulsion: this.config.simulationRepulsion, + maxStep: 0, + }, + }, + }) + + this.forceAllPairsCommand ||= new Model(device, { + fs: forceAllPairsFrag, + vs: updateVert, + topology: 'triangle-strip', + vertexCount: 4, + attributes: { + vertexCoord: this.forceVertexCoordBuffer, + }, + bufferLayout: [ + { name: 'vertexCoord', format: 'float32x2' }, + ], + defines: { + USE_UNIFORM_BUFFERS: true, + }, + bindings: { + forceAllPairsUniforms: this.forceAllPairsUniformStore.getManagedUniformBuffer('forceAllPairsUniforms'), + // All texture bindings will be set dynamically in drawAllPairsForce() method + }, + parameters: { + blend: false, + depthWriteEnabled: false, + depthCompare: 'always', + }, + }) + this.forceNearFieldCommand ||= new Model(device, { fs: forceNearFieldFrag, vs: updateVert, @@ -407,9 +536,15 @@ export class ForceManyBody extends CoreModule { return } + // Small graphs: one exact all-pairs pass, no grid, no sampling. + if (this.usesAllPairs) { + this.drawAllPairsForce() + return + } + // Nothing to do until the grid pyramid and near-field slots are allocated // (create() builds them; this guards a partial/failed allocation). - if (this.levelTargets.size === 0 || this.nearFieldSlotTargets.length !== NEAR_FIELD_SLOTS) return + if (this.levelTargets.size === 0 || this.peelTargets.length !== PEEL_TARGETS || !this.slotsArrayTexture) return this.drawLevels() this.drawNearFieldSlots() @@ -430,6 +565,8 @@ export class ForceManyBody extends CoreModule { this.buildNearFieldSlotsCommand = undefined this.forceNearFieldCommand?.destroy() this.forceNearFieldCommand = undefined + this.forceAllPairsCommand?.destroy() + this.forceAllPairsCommand = undefined // 2. Destroy Framebuffers + 3. Textures (grid targets destroy their FBOs // before their textures internally) @@ -448,6 +585,8 @@ export class ForceManyBody extends CoreModule { this.buildNearFieldSlotsUniformStore = undefined this.forceNearFieldUniformStore?.destroy() this.forceNearFieldUniformStore = undefined + this.forceAllPairsUniformStore?.destroy() + this.forceAllPairsUniformStore = undefined // 5. Destroy Buffers (passed via attributes - NOT owned by Models, must destroy manually) if (this.pointIndices && !this.pointIndices.destroyed) { @@ -460,6 +599,47 @@ export class ForceManyBody extends CoreModule { this.forceVertexCoordBuffer = undefined } + /** + * The small-graph path: a single exact all-pairs pass into the velocity FBO. + * Replaces the pyramid + near-field passes below ALL_PAIRS_MAX_POINTS. + */ + private drawAllPairsForce (): void { + const { device, store, data, points } = this + if (!points) return + if (!this.forceAllPairsCommand || !this.forceAllPairsUniformStore) return + if (!points.previousPositionTexture || points.previousPositionTexture.destroyed) return + if (!points.exitTexture || points.exitTexture.destroyed) return + if (!this.randomValuesTexture || this.randomValuesTexture.destroyed) return + if (!points.velocityFbo || points.velocityFbo.destroyed) return + if (!data.pointsNumber) return + + this.forceAllPairsUniformStore.setUniforms({ + forceAllPairsUniforms: { + pointsTextureSize: store.pointsTextureSize ?? 0, + pointsNumber: data.pointsNumber, + alpha: store.alpha, + repulsion: this.config.simulationRepulsion, + // The near-field pass's per-tick bound and the shader's near/far split + // radius, computed from the finest cell size the grid path would use + // at this point count. + maxStep: 2 * (store.adjustedSpaceSize / getFinestGridSize(data.pointsNumber)), + }, + }) + + this.forceAllPairsCommand.setBindings({ + positionsTexture: points.previousPositionTexture, + randomValues: this.randomValuesTexture, + exitTexture: points.exitTexture, + }) + + const drawPass = device.beginRenderPass({ + framebuffer: points.velocityFbo, + clearColor: [0, 0, 0, 0], + }) + this.forceAllPairsCommand.draw(drawPass) + drawPass.end() + } + /** Aggregates points into every grid level texture. */ private drawLevels (): void { const { device, store, data, points } = this @@ -502,11 +682,13 @@ export class ForceManyBody extends CoreModule { } /** - * Rebuilds the near-field point slots for this tick: NEAR_FIELD_SLOTS + * Rebuilds the near-field point slots for this tick: `nearFieldSlots` * depth-peeling passes over the points, each capturing the eligible point with * the smallest per-tick random hash per finest-level cell (see * build-nearfield-slots.vert). Re-seeded every tick so dense cells rotate all - * their points through exact pairwise treatment. + * their points through exact pairwise treatment. Each pass ping-pongs between + * the two peel targets (reading the previous pass's output), then its result + * is copied into its layer of the slot array texture. */ private drawNearFieldSlots (): void { const { device, store, data, points } = this @@ -515,14 +697,16 @@ export class ForceManyBody extends CoreModule { if (!points.previousPositionTexture || points.previousPositionTexture.destroyed) return if (!points.exitTexture || points.exitTexture.destroyed) return if (!data.pointsNumber || !this.pointIndices) return + if (!this.slotsArrayTexture || this.slotsArrayTexture.destroyed) return const finest = this.levelTargets.get(this.levels - 1) if (!finest || finest.texture.destroyed) return const randomSeed = store.getRandomFloat(0, 1) - for (let slot = 0; slot < this.nearFieldSlotTargets.length; slot += 1) { - const target = this.nearFieldSlotTargets[slot] - if (!target || target.fbo.destroyed) continue + for (let slot = 0; slot < this.nearFieldSlots; slot += 1) { + const target = this.peelTargets[slot % PEEL_TARGETS] + const previous = this.peelTargets[(slot + 1) % PEEL_TARGETS] + if (!target || target.fbo.destroyed || !previous || previous.texture.destroyed) continue this.buildNearFieldSlotsUniformStore.setUniforms({ buildNearFieldSlotsUniforms: { @@ -544,7 +728,7 @@ export class ForceManyBody extends CoreModule { // draw to run — any texture that is not the render target works. previousSlot: slot === 0 ? points.previousPositionTexture - : this.nearFieldSlotTargets[slot - 1]!.texture, + : previous.texture, }) const slotPass = device.beginRenderPass({ @@ -555,6 +739,19 @@ export class ForceManyBody extends CoreModule { }) this.buildNearFieldSlotsCommand.draw(slotPass) slotPass.end() + + // Publish this pass's result as layer `slot` of the array texture that + // the near-field force pass samples. + const commandEncoder = device.createCommandEncoder() + commandEncoder.copyTextureToTexture({ + sourceTexture: target.texture, + destinationTexture: this.slotsArrayTexture, + destinationOrigin: [0, 0, slot], + width: finest.gridSize, + height: finest.gridSize, + }) + // finish() destroys the encoder itself and returns the command buffer. + device.submit(commandEncoder.finish()) } } @@ -567,7 +764,8 @@ export class ForceManyBody extends CoreModule { if (!points) return if (!this.forceLevelCommand || !this.forceLevelUniformStore) return if (!this.forceNearFieldCommand || !this.forceNearFieldUniformStore) return - if (this.nearFieldSlotTargets.length !== NEAR_FIELD_SLOTS) return + if (this.peelTargets.length !== PEEL_TARGETS) return + if (!this.slotsArrayTexture || this.slotsArrayTexture.destroyed) return if (!points.previousPositionTexture || points.previousPositionTexture.destroyed) return if (!this.randomValuesTexture || this.randomValuesTexture.destroyed) return if (!points.velocityFbo || points.velocityFbo.destroyed) return @@ -610,25 +808,15 @@ export class ForceManyBody extends CoreModule { cellSize, alpha: store.alpha, repulsion: this.config.simulationRepulsion, + slotCount: this.nearFieldSlots, }, }) - // One binding per slot, listed out to match the samplers over in - // force-nearfield.frag — both lists have to stay NEAR_FIELD_SLOTS long. - // We only reach here after confirming we have a full set of slots (check - // at the top of this method), so the `!` on each one is safe. this.forceNearFieldCommand.setBindings({ positionsTexture: points.previousPositionTexture, levelTexture: target.texture, randomValues: this.randomValuesTexture, - slotTexture0: this.nearFieldSlotTargets[0]!.texture, - slotTexture1: this.nearFieldSlotTargets[1]!.texture, - slotTexture2: this.nearFieldSlotTargets[2]!.texture, - slotTexture3: this.nearFieldSlotTargets[3]!.texture, - slotTexture4: this.nearFieldSlotTargets[4]!.texture, - slotTexture5: this.nearFieldSlotTargets[5]!.texture, - slotTexture6: this.nearFieldSlotTargets[6]!.texture, - slotTexture7: this.nearFieldSlotTargets[7]!.texture, + slotsTexture: this.slotsArrayTexture, }) this.forceNearFieldCommand.draw(drawPass) } @@ -647,11 +835,7 @@ export class ForceManyBody extends CoreModule { const { device } = this const pointsNumber = this.data.pointsNumber ?? 0 - const targetGridSize = 2 * Math.sqrt(pointsNumber) - const finestGridSize = Math.min( - MAX_GRID_SIZE, - Math.max(8, Math.pow(2, Math.ceil(Math.log2(targetGridSize)))) - ) + const finestGridSize = getFinestGridSize(pointsNumber) this.levels = Math.log2(finestGridSize) - 1 for (let level = 0; level < this.levels; level += 1) { @@ -694,33 +878,38 @@ export class ForceManyBody extends CoreModule { } /** - * Allocates the depth-peeling slot targets ([point index, hash] per cell) plus - * a depth attachment each for the peel's smallest-hash selection. + * Allocates the near-field sampling resources: the two ping-pong depth-peeling + * targets ([point index, hash] per cell, with a depth attachment each for the + * peel's smallest-hash selection) and the slot array texture (one layer per + * peeling pass) that the force pass samples. */ private createNearFieldSlotTargets (finest: LevelTarget): void { const { device } = this - // These slots follow the finest level's grid, and that grid does change size - // as the graph grows or shrinks (it snaps to powers of two). So unlike the - // level targets, we really might need to resize here: if what we already have - // matches the finest grid and we've got all NEAR_FIELD_SLOTS, keep it; - // otherwise throw it away and rebuild at the new size. Every slot is the same - // size, so checking slot 0 tells us about the whole set. - const existing = this.nearFieldSlotTargets[0] + const slots = getNearFieldSlotCount(this.data.pointsNumber ?? 0) + // These targets follow the finest level's grid, and that grid does change + // size as the graph grows or shrinks (it snaps to powers of two) — and the + // slot count changes with the point count too. If everything we already + // have matches, keep it; otherwise throw it away and rebuild. + const existing = this.peelTargets[0] if ( existing && !existing.texture.destroyed && existing.texture.width === finest.gridSize && existing.texture.height === finest.gridSize && - this.nearFieldSlotTargets.length === NEAR_FIELD_SLOTS + this.peelTargets.length === PEEL_TARGETS && + this.slotsArrayTexture && + !this.slotsArrayTexture.destroyed && + this.nearFieldSlots === slots ) return this.destroyNearFieldSlotTargets() - for (let slot = 0; slot < NEAR_FIELD_SLOTS; slot += 1) { + this.nearFieldSlots = slots + for (let target = 0; target < PEEL_TARGETS; target += 1) { const texture = device.createTexture({ width: finest.gridSize, height: finest.gridSize, format: 'rg32float', - usage: Texture.SAMPLE | Texture.RENDER, + usage: Texture.SAMPLE | Texture.RENDER | Texture.COPY_SRC, }) const fbo = device.createFramebuffer({ width: finest.gridSize, @@ -733,16 +922,29 @@ export class ForceManyBody extends CoreModule { // smallest-hash point from the whole tick's sample. depthStencilAttachment: 'depth24plus', }) - this.nearFieldSlotTargets.push({ texture, fbo }) + this.peelTargets.push({ texture, fbo }) } + this.slotsArrayTexture = device.createTexture({ + dimension: '2d-array', + width: finest.gridSize, + height: finest.gridSize, + depth: slots, + format: 'rg32float', + usage: Texture.SAMPLE | Texture.COPY_DST, + }) } private destroyNearFieldSlotTargets (): void { - for (const target of this.nearFieldSlotTargets) { + for (const target of this.peelTargets) { if (!target.fbo.destroyed) target.fbo.destroy() if (!target.texture.destroyed) target.texture.destroy() } - this.nearFieldSlotTargets = [] + this.peelTargets = [] + if (this.slotsArrayTexture && !this.slotsArrayTexture.destroyed) { + this.slotsArrayTexture.destroy() + } + this.slotsArrayTexture = undefined + this.nearFieldSlots = 0 } private destroyLevelTargets (): void { diff --git a/src/stories/performance.stories.ts b/src/stories/performance.stories.ts index b98cc9d3..b4e27e89 100644 --- a/src/stories/performance.stories.ts +++ b/src/stories/performance.stories.ts @@ -7,6 +7,7 @@ import { collisionStressTest } from './performance/collision-stress-test' import { pointOcclusionCulling } from './performance/point-occlusion-culling' import { onDemandRendering } from './performance/on-demand-rendering' import { repulsionBenchmark } from './performance/repulsion-benchmark' +import { countryBordersComparison } from './performance/country-borders-comparison' import createCosmosRaw from './create-cosmos?raw' import generateMeshDataRaw from './generate-mesh-data?raw' @@ -16,6 +17,8 @@ import collisionStressTestRaw from './performance/collision-stress-test?raw' import pointOcclusionCullingRaw from './performance/point-occlusion-culling?raw' import onDemandRenderingRaw from './performance/on-demand-rendering?raw' import repulsionBenchmarkRaw from './performance/repulsion-benchmark?raw' +import countryBordersComparisonRaw from './performance/country-borders-comparison?raw' +import countryBordersDataRaw from './performance/country-borders-data?raw' // These exist to show a cost or a limit rather than a feature. Most run an FPS // monitor; several will be slow to start on weak hardware. @@ -87,5 +90,17 @@ export const RepulsionBenchmark: Story = { }, } +export const CountryBordersComparison: Story = { + ...createStory(countryBordersComparison), + name: 'Repulsion Jitter: Fixed vs Before', + tags: ['perf', 'advanced', 'interactive'], + parameters: { + sourceCode: [ + { name: 'Story', code: countryBordersComparisonRaw }, + { name: 'Data', code: countryBordersDataRaw }, + ], + }, +} + // eslint-disable-next-line import/no-default-export export default meta diff --git a/src/stories/performance/country-borders-comparison.ts b/src/stories/performance/country-borders-comparison.ts new file mode 100644 index 00000000..d7f0a754 --- /dev/null +++ b/src/stories/performance/country-borders-comparison.ts @@ -0,0 +1,398 @@ +import { Graph, GraphConfig } from '@cosmos.gl/graph' +import { ForceManyBody } from '@/graph/modules/ForceManyBody' +import { COUNTRY_BORDER_LINKS } from './country-borders-data' + +// What the many-body jitter fix changed, live: the country borders network +// (163 points, 642 links) run twice with identical data, seed, and simulation +// settings. Left is cosmos.gl today — small graphs take the exact all-pairs +// path, and the settled layout is still. Right is how every graph ran before +// the fix — the grid + Monte-Carlo near field with K = 8 sampling slots, +// forced back on here purely for comparison. Each side reports step/turn over +// a sliding window and traces one point from its densest finest cell: a +// smooth drift arc today, a random-walk tangle before. +// +// ⚠ HOW THE FORCING WORKS — READ BEFORE COPYING. The engine deliberately does +// not expose the path choice or the slot count; this story reaches into +// unexported internals (`ForceManyBody`) through the repo's `@/graph` source +// alias and patches two private members at module load. That import only +// resolves inside this repository — none of the patch code works against the +// published package, and any refactor of those internals may break it (the +// patch throws loudly if the members disappear). It exists so the comparison +// is reproducible without engine changes; it is not a supported API. + +const SEED = 42 +const SPACE = 4096 +const METER_WINDOW = 120 // ticks +const TRAIL_WINDOW = 360 // ticks of path history in the trajectory panel +const TRAIL_PICK_TICK = 180 // pick the traced point after the clump has formed + +// ── Internals patch (repo-only, see header) ───────────────────────────────── + +// Per-instance marker smuggled through the config: applyConfig copies unknown +// keys verbatim onto the merged config object, and every force module holds a +// reference to that object — so the patched members below can tell the two +// twin graphs apart at call time. +const FORCE_SAMPLED_KEY = '__storyForceSampledRepulsion' +// Symbol.for survives HMR re-execution, keeping the patch single-layered. +const PATCH_FLAG = Symbol.for('cosmos.stories.force-sampled-repulsion-patch') + +type ForceManyBodyPrivate = { + config: Record; + data: { pointsNumber?: number }; + nearFieldSlots?: number; +} + +const patchForceManyBody = (): void => { + const proto = ForceManyBody.prototype as unknown as Record + if (proto[PATCH_FLAG]) return + + const usesAllPairs = Object.getOwnPropertyDescriptor(proto, 'usesAllPairs') + const createSlotTargets = Object.getOwnPropertyDescriptor(proto, 'createNearFieldSlotTargets') + if (typeof usesAllPairs?.get !== 'function' || typeof createSlotTargets?.value !== 'function') { + throw new Error('country-borders-comparison: ForceManyBody internals moved — update this story') + } + const originalUsesAllPairs = usesAllPairs.get + const originalCreateSlotTargets = createSlotTargets.value as (...args: unknown[]) => unknown + + // Marked instances always take the grid + sampled near-field path. + Object.defineProperty(proto, 'usesAllPairs', { + configurable: true, + get (this: ForceManyBodyPrivate): boolean { + if (this.config[FORCE_SAMPLED_KEY]) return false + return originalUsesAllPairs.call(this) as boolean + }, + }) + + // getNearFieldSlotCount is module-private and can't be patched directly; it + // reads data.pointsNumber only inside this method, so shadow the getter with + // a count from its 8-slot tier for the duration of the call. + Object.defineProperty(proto, 'createNearFieldSlotTargets', { + configurable: true, + writable: true, + value (this: ForceManyBodyPrivate, ...args: unknown[]): unknown { + if (!this.config[FORCE_SAMPLED_KEY]) return originalCreateSlotTargets.apply(this, args) + Object.defineProperty(this.data, 'pointsNumber', { value: 100_000, configurable: true }) + try { + const result = originalCreateSlotTargets.apply(this, args) + // The pane's "K = 8" header must stay true: fail loudly if the slot + // tiers are retuned so 100k points no longer maps to 8 slots, or the + // count stops flowing through data.pointsNumber. + if (this.nearFieldSlots !== 8) { + throw new Error(`country-borders-comparison: expected 8 near-field slots, got ${this.nearFieldSlots} — update this story`) + } + return result + } finally { + delete this.data.pointsNumber + } + }, + }) + + proto[PATCH_FLAG] = true +} + +// ── Shared data ───────────────────────────────────────────────────────────── + +const mulberry32 = (seed: number): (() => number) => { + let a = seed >>> 0 + return () => { + a = (a + 0x6D2B79F5) >>> 0 + let t = a + t = Math.imul(t ^ (t >>> 15), t | 1) + t ^= t + Math.imul(t ^ (t >>> 7), t | 61) + return ((t ^ (t >>> 14)) >>> 0) / 4294967296 + } +} + +const parseLinks = (): { count: number; links: Float32Array } => { + const indexOf = new Map() + const pairs = COUNTRY_BORDER_LINKS.trim().split(/\s+/) + const links = new Float32Array(pairs.length * 2) + for (const [i, pair] of pairs.entries()) { + const [source, target] = pair!.split('-') as [string, string] + for (const code of [source, target]) { + if (!indexOf.has(code)) indexOf.set(code, indexOf.size) + } + links[i * 2] = indexOf.get(source)! + links[i * 2 + 1] = indexOf.get(target)! + } + return { count: indexOf.size, links } +} + +// ── One instrumented pane (graph + meter + trajectory panel) ──────────────── + +type Pane = { + element: HTMLDivElement; + graph: Graph; + destroy: () => void; +} + +const createPane = (title: string, accent: string, forceSampled: boolean): Pane => { + const pane = document.createElement('div') + pane.style.cssText = 'flex:1;min-width:0;display:flex;flex-direction:column;overflow:hidden;' + + const header = document.createElement('div') + header.style.cssText = 'padding:10px 14px;font-size:12px;line-height:1.6;flex:none;' + const titleLine = document.createElement('div') + titleLine.style.cssText = `font-weight:bold;color:${accent};` + titleLine.textContent = title + const meter = document.createElement('div') + meter.textContent = 'Warming up…' + header.appendChild(titleLine) + header.appendChild(meter) + pane.appendChild(header) + + const stage = document.createElement('div') + stage.style.cssText = 'flex:1;min-height:200px;position:relative;' + pane.appendChild(stage) + + const host = document.createElement('div') + host.style.cssText = 'height:100%;' + stage.appendChild(host) + + const trailPanel = document.createElement('div') + trailPanel.style.cssText = + 'position:absolute;top:10px;right:10px;width:200px;background:rgba(26,29,35,0.88);border:1px solid #3a3f47;padding:8px;pointer-events:none;' + const trailCanvas = document.createElement('canvas') + trailCanvas.style.cssText = 'display:block;width:184px;height:184px;' + const trailLabel = document.createElement('div') + trailLabel.style.cssText = 'margin-top:6px;font-size:11px;line-height:1.4;color:#9aa3ad;' + trailLabel.textContent = 'Trajectory: waiting for the clump to form…' + trailPanel.appendChild(trailCanvas) + trailPanel.appendChild(trailLabel) + stage.appendChild(trailPanel) + + const { count: n, links } = parseLinks() + const rng = mulberry32(SEED) + const positions = new Float32Array(n * 2) + for (let i = 0; i < n; i += 1) { + positions[i * 2] = SPACE * (0.25 + rng() * 0.5) + positions[i * 2 + 1] = SPACE * (0.25 + rng() * 0.5) + } + + // Finest many-body grid for this point count (mirrors ForceManyBody.createLevels). + const finestGrid = Math.min(512, Math.max(8, 2 ** Math.ceil(Math.log2(2 * Math.sqrt(n))))) + const cell = SPACE / finestGrid + const cellOf = (x: number, y: number): number => { + const cx = Math.min(finestGrid - 1, Math.max(0, Math.floor(x / cell))) + const cy = Math.min(finestGrid - 1, Math.max(0, Math.floor(y / cell))) + return cy * finestGrid + cx + } + + let previous: number[] | null = null + let previousStep: Float32Array | null = null + const stepWindow: number[] = [] + const turnWindow: number[] = [] + let tick = 0 + let tracedIndex = -1 + const trail: number[] = [] + + const pickFromDensestCell = (tracked: number[]): number => { + const counts = new Map() + for (let i = 0; i < n; i += 1) { + const key = cellOf(tracked[i * 2]!, tracked[i * 2 + 1]!) + counts.set(key, (counts.get(key) ?? 0) + 1) + } + let bestKey = -1 + let bestCount = 0 + for (const [key, count] of counts) { + if (count > bestCount) { + bestCount = count + bestKey = key + } + } + for (let i = 0; i < n; i += 1) { + if (cellOf(tracked[i * 2]!, tracked[i * 2 + 1]!) === bestKey) return i + } + return 0 + } + + const niceUnits = (target: number): number => { + const pow = 10 ** Math.floor(Math.log10(target)) + for (const mult of [5, 2, 1]) { + if (mult * pow <= target) return mult * pow + } + return pow + } + + const drawTrail = (): void => { + const ctx = trailCanvas.getContext('2d') + if (!ctx || trail.length < 4) return + const size = 184 + const dpr = window.devicePixelRatio || 1 + if (trailCanvas.width !== size * dpr) { + trailCanvas.width = size * dpr + trailCanvas.height = size * dpr + } + ctx.setTransform(dpr, 0, 0, dpr, 0, 0) + ctx.clearRect(0, 0, size, size) + + let minX = Infinity + let maxX = -Infinity + let minY = Infinity + let maxY = -Infinity + for (let i = 0; i < trail.length; i += 2) { + minX = Math.min(minX, trail[i]!) + maxX = Math.max(maxX, trail[i]!) + minY = Math.min(minY, trail[i + 1]!) + maxY = Math.max(maxY, trail[i + 1]!) + } + const span = Math.max(maxX - minX, maxY - minY) + // Auto-fit, but never zoom past 1 unit per panel — a still point must read + // as a dot, not be inflated into a false tangle. + const fitSpan = Math.max(span, 1) + const scale = (size - 24) / fitSpan + const midX = (minX + maxX) / 2 + const midY = (minY + maxY) / 2 + const toX = (x: number): number => size / 2 + (x - midX) * scale + const toY = (y: number): number => size / 2 + (y - midY) * scale + + ctx.lineWidth = 1.5 + for (let i = 2; i < trail.length; i += 2) { + ctx.strokeStyle = `rgba(233,161,59,${(0.15 + 0.85 * (i / trail.length)).toFixed(3)})` + ctx.beginPath() + ctx.moveTo(toX(trail[i - 2]!), toY(trail[i - 1]!)) + ctx.lineTo(toX(trail[i]!), toY(trail[i + 1]!)) + ctx.stroke() + } + ctx.fillStyle = '#f3c063' + ctx.beginPath() + ctx.arc(toX(trail[trail.length - 2]!), toY(trail[trail.length - 1]!), 3, 0, Math.PI * 2) + ctx.fill() + + const barUnits = niceUnits((0.4 * size) / scale) + ctx.strokeStyle = '#9aa3ad' + ctx.lineWidth = 1 + ctx.beginPath() + ctx.moveTo(10, size - 10.5) + ctx.lineTo(10 + barUnits * scale, size - 10.5) + ctx.stroke() + ctx.fillStyle = '#9aa3ad' + ctx.font = '10px monospace' + ctx.fillText(`${barUnits} u`, 12, size - 15) + + trailLabel.textContent = + `Trajectory: point ${tracedIndex} (densest cell), last ${trail.length / 2} ticks — path spans ${span.toFixed(2)} u` + } + + const config: GraphConfig = { + spaceSize: SPACE, + fitViewOnInit: true, + fitViewPadding: 0.3, + rescalePositions: false, + randomSeed: SEED, + pointDefaultSize: 6, + linkDefaultWidth: 1, + renderLinks: true, + // Hold alpha ≈ 1 so the sampled side's shimmer persists instead of + // annealing away — what a long-running / reheated layout experiences. + simulationDecay: 1e12, + onSimulationTick: (): void => { + tick += 1 + const tracked = graph.getTrackedPointPositionsArray() + if (tracked.length !== n * 2) return + if (previous) { + const steps = new Float32Array(n * 2) + let stepSum = 0 + let turnSum = 0 + let turnCount = 0 + for (let i = 0; i < n; i += 1) { + const dx = tracked[i * 2]! - previous[i * 2]! + const dy = tracked[i * 2 + 1]! - previous[i * 2 + 1]! + steps[i * 2] = dx + steps[i * 2 + 1] = dy + const len = Math.hypot(dx, dy) + stepSum += len + if (previousStep) { + const px = previousStep[i * 2]! + const py = previousStep[i * 2 + 1]! + const plen = Math.hypot(px, py) + if (len > 1e-9 && plen > 1e-9) { + const cos = Math.min(1, Math.max(-1, (dx * px + dy * py) / (len * plen))) + turnSum += Math.acos(cos) * (180 / Math.PI) + turnCount += 1 + } + } + } + stepWindow.push(stepSum / n) + if (turnCount > 0) turnWindow.push(turnSum / turnCount) + if (stepWindow.length > METER_WINDOW) stepWindow.shift() + if (turnWindow.length > METER_WINDOW) turnWindow.shift() + previousStep = steps + } + previous = tracked + + if (tracedIndex < 0 && tick >= TRAIL_PICK_TICK) tracedIndex = pickFromDensestCell(tracked) + if (tracedIndex >= 0) { + trail.push(tracked[tracedIndex * 2]!, tracked[tracedIndex * 2 + 1]!) + if (trail.length > TRAIL_WINDOW * 2) trail.splice(0, trail.length - TRAIL_WINDOW * 2) + drawTrail() + } + + if (tick % 15 === 0 && stepWindow.length > 0) { + const meanStep = stepWindow.reduce((a, b) => a + b, 0) / stepWindow.length + const meanTurn = turnWindow.length > 0 ? turnWindow.reduce((a, b) => a + b, 0) / turnWindow.length : 0 + meter.textContent = + `last ${stepWindow.length} ticks: step ${meanStep.toFixed(2)} u/tick turn ${meanTurn.toFixed(1)}°` + } + }, + } + // The per-instance marker read by the patched internals (see the header). + if (forceSampled) (config as Record)[FORCE_SAMPLED_KEY] = true + + const graph = new Graph(host, config) + + let destroyed = false + const setup = async (): Promise => { + await graph.ready + if (destroyed) return + graph.setPointPositions(positions, true) + graph.setLinks(links) + graph.render() + graph.trackPointPositionsByIndices(Array.from({ length: n }, (_, i) => i)) + graph.start(1) + // Re-fit while the layout finds its equilibrium shape, then leave the + // camera alone. + for (const delay of [1500, 3500, 6000]) { + setTimeout(() => { + if (!destroyed) graph.fitView(400) + }, delay) + } + } + // eslint-disable-next-line no-console + setup().catch((error) => console.error('[country-borders-comparison] failed', error)) + + return { + element: pane, + graph, + destroy: (): void => { + destroyed = true + }, + } +} + +// ── The story: both paths, same data, same tick ───────────────────────────── + +export const countryBordersComparison = (): { graph: Graph; div: HTMLDivElement; destroy?: () => void } => { + patchForceManyBody() + + const outer = document.createElement('div') + outer.style.cssText = 'height:100vh;width:100%;background:#1a1d23;color:#e0e0e0;font-family:monospace;display:flex;overflow:hidden;' + + const left = createPane('Today: exact all-pairs repulsion — the current engine', '#7fd1c0', false) + const right = createPane('Before the fix: sampled near field, K = 8 — forced for comparison', '#e9a13b', true) + right.element.style.borderLeft = '1px solid #3a3f47' + outer.appendChild(left.element) + outer.appendChild(right.element) + + return { + graph: left.graph, + div: outer, + destroy: (): void => { + // The story contract tears down `graph` (the left pane) itself; the + // right pane's graph is this story's own cleanup responsibility. + left.destroy() + right.destroy() + right.graph.destroy() + }, + } +} diff --git a/src/stories/performance/country-borders-data.ts b/src/stories/performance/country-borders-data.ts new file mode 100644 index 00000000..9d5641e1 --- /dev/null +++ b/src/stories/performance/country-borders-data.ts @@ -0,0 +1,45 @@ +// Real-world border adjacencies between countries (ISO 3166-1 alpha-2 codes), +// exported 2026-08-15. Each token is one directed source-target link; most +// borders appear in both directions, exactly as the source dataset lists them. +export const COUNTRY_BORDER_LINKS: string = + 'AZ-TR DE-CH GE-TR HR-SI KG-UZ KR-KP MO-CN SN-MR SZ-ZA UA-SK DZ-EH GY-VE HU-UA KW-SA LT-RU MC-FR ' + + 'MR-EH MW-ZM NO-SE SD-SS SM-IT SV-HN TN-LY BE-NL CA-US CZ-SK DK-DE MK-RS QA-SA AL-RS AT-CH CM-NG ' + + 'ER-SD GF-SR GT-MX PL-UA PY-BR SA-YE BO-PE BW-ZW EC-PE GB-IE GH-TG GQ-GA ID-TL LB-SY MM-TH NL-DE ' + + 'SE-NO SL-LR TZ-ZM VN-LA FI-SE KP-RU KZ-UZ LA-VN LS-ZA ME-RS MF-SX NI-HN UY-BR ZM-ZW AF-UZ CG-GA ' + + 'CH-LI CL-PE CN-VN CR-PA GI-ES IL-SY IN-PK MY-TH NE-NG RO-UA RW-UG US-MX AD-ES AM-TR AR-UY BN-MY ' + + 'CF-SD ET-SD HN-NI JO-SY LV-RU SI-IT VA-IT AE-SA AO-ZM BG-TR BY-UA HK-CN IE-GB LI-CH LY-TN MN-RU ' + + 'TD-SD TH-MM UZ-TM BD-MM EG-SD GA-GQ IQ-TR ML-SN NP-IN RS-RO UG-TZ CI-ML EH-MA GM-SN HT-DO MX-US ' + + 'PE-EC PS-JO YE-SA BT-IN BZ-MX CD-ZM GW-SN IR-TM KH-VN MZ-ZW OM-YE TL-ID TM-UZ BF-TG BI-TZ BJ-TG ' + + 'EE-RU FR-CH GR-TR PT-ES SK-UA SO-KE SR-GY TJ-UZ TR-SY VE-GY CO-VE DJ-SO IT-VA LR-SL MA-EH MD-UA ' + + 'NA-ZM NG-NE PG-ID RU-UA SS-UG TG-GH BR-VE DO-HT KE-UG LU-FR PA-CR PK-IR SY-TR ZA-ZW BA-RS ES-PT ' + + 'GN-SL SX-MF ZW-ZM AZ-RU DE-PL GE-RU HR-RS KG-TJ SN-ML SZ-MZ UA-RU DZ-TN GY-SR HU-SI KW-IQ LT-PL ' + + 'MR-SN MW-TZ NO-RU SD-LY SV-GT TN-DZ BE-LU CZ-PL MK-GR AL-MK AT-SI CM-GA ER-ET GF-BR GT-HN PL-SK ' + + 'PY-BO SA-AE BO-PY BW-ZM EC-CO GH-CI GQ-CM ID-PG LB-IL MM-LA NL-BE SE-FI SL-GN TZ-UG VN-CN FI-RU ' + + 'KP-KR KZ-TM LA-TH ME-HR NI-CR UY-AR ZM-TZ AF-TM CG-CD CH-IT CL-BO CN-TJ CR-NI IL-PS IN-NP MY-ID ' + + 'NE-ML RO-RS RW-TZ US-CA AD-FR AM-IR AR-PY CF-SS ET-SS HN-GT JO-SA LV-LT SI-HU AE-OM AO-NA BG-RS ' + + 'BY-RU LI-AT LY-SD MN-CN TD-NG TH-MY UZ-TJ BD-IN EG-PS GA-CG IQ-SY ML-NE NP-CN RS-MK UG-SS CI-LR ' + + 'EH-MR MX-GT PE-CO PS-IL YE-OM BT-CN BZ-GT CD-UG GW-GN IR-TR KH-TH MZ-ZM OM-SA TM-KZ BF-NE BI-RW ' + + 'BJ-NG EE-LV FR-ES GR-MK SK-PL SO-ET SR-GF TJ-KG TR-IQ VE-CO CO-PE DJ-ET IT-CH LR-GN MA-ES MD-RO ' + + 'NA-ZA NG-TD RU-PL SS-SD TG-BF BR-UY KE-TZ LU-DE PA-CO PK-IN SY-LB ZA-SZ BA-ME ES-MA GN-SN ZW-ZA ' + + 'AZ-IR DE-NL GE-AZ HR-ME KG-KZ SN-GW UA-RO DZ-NE GY-BR HU-SK LT-LV MR-ML MW-MZ NO-FI SD-ER BE-DE ' + + 'CZ-DE MK-BG AL-ME AT-SK CM-GQ ER-DJ GT-SV PL-RU PY-AR SA-QA BO-CL BW-ZA GH-BF ID-MY MM-IN TZ-RW ' + + 'VN-KH FI-NO KP-CN KZ-RU LA-MM ME-BA ZM-NA AF-TJ CG-CF CH-DE CL-AR CN-RU IL-LB IN-MM MY-BN NE-LY ' + + 'RO-MD RW-CD AM-GE AR-CL CF-CD ET-SO HN-SV JO-PS LV-EE SI-HR AO-CD BG-RO BY-PL LY-NE TD-NE TH-LA ' + + 'UZ-KG EG-LY GA-CM IQ-SA ML-MR RS-ME UG-RW CI-GN EH-DZ MX-BZ PE-CL PS-EG CD-TZ IR-PK KH-LA MZ-TZ ' + + 'OM-AE TM-IR BF-ML BI-CD BJ-NE FR-MC GR-BG SK-HU SO-DJ SR-BR TJ-CN TR-IR VE-BR CO-PA DJ-ER IT-SI ' + + 'LR-CI MA-DZ NA-BW NG-CM RU-NO SS-KE TG-BJ BR-SR KE-SS LU-BE PK-CN SY-JO ZA-NA BA-HR ES-GI GN-ML ' + + 'ZW-MZ AZ-GE DE-LU GE-AM HR-HU KG-CN SN-GN UA-PL DZ-MA HU-RS LT-BY MR-DZ SD-ET BE-FR CZ-AT MK-AL ' + + 'AL-GR AT-LI CM-CG GT-BZ PL-LT SA-OM BO-BR BW-NA MM-CN TZ-MZ KZ-KG LA-KH ME-AL ZM-MZ AF-PK CG-CM ' + + 'CH-FR CN-PK IL-JO IN-CN NE-TD RO-HU RW-BI AM-AZ AR-BR CF-CG ET-KE JO-IL LV-BY SI-AT AO-CG BG-MK ' + + 'BY-LT LY-EG TD-LY TH-KH UZ-KZ EG-IL IQ-KW ML-GN RS-HU UG-KE CI-GH PE-BR CD-SS IR-IQ MZ-ZA TM-AF ' + + 'BF-GH BJ-BF FR-LU GR-AL SK-CZ TJ-AF TR-GR CO-EC IT-SM NA-AO NG-BJ RU-MN SS-ET BR-PE KE-SO PK-AF ' + + 'SY-IL ZA-MZ ES-FR GN-LR ZW-BW AZ-AM DE-FR HR-BA SN-GM UA-MD DZ-MR HU-RO SD-EG AT-IT CM-TD PL-DE ' + + 'SA-KW BO-AR MM-BD TZ-MW KZ-CN LA-CN ZM-MW AF-IR CG-AO CH-AT CN-NP IL-EG IN-BT NE-BF RO-BG AR-BO ' + + 'CF-TD ET-ER JO-IQ BG-GR BY-LV LY-TD TD-CF UZ-AF IQ-JO ML-CI RS-HR UG-CD CI-BF PE-BO CD-RW IR-AZ ' + + 'MZ-SZ BF-CI FR-IT SK-AT TR-GE CO-BR IT-FR RU-LT SS-CD BR-PY KE-ET SY-IQ ZA-LS ES-AD GN-GW DE-DK ' + + 'UA-HU DZ-ML HU-HR SD-TD AT-HU CM-CF PL-CZ SA-JO TZ-KE ZM-CD AF-CN CN-MM IN-BD NE-BJ CF-CM ET-DJ ' + + 'LY-DZ TD-CM IQ-IR ML-BF RS-BG CD-CG IR-AM MZ-MW BF-BJ FR-DE TR-BG IT-AT RU-LV SS-CF BR-GY ZA-BW ' + + 'GN-CI DE-CZ UA-BY DZ-LY HU-AT SD-CF AT-DE PL-BY SA-IQ TZ-CD ZM-BW CN-MN NE-DZ ML-DZ RS-BA CD-CF ' + + 'IR-AF FR-BE TR-AZ RU-KP BR-GF DE-BE AT-CZ TZ-BI ZM-AO CN-MO RS-AL CD-BI FR-AD TR-AM RU-KZ BR-CO ' + + 'DE-AT CN-LA CD-AO RU-GE BR-BO CN-KG RU-FI BR-AR CN-KP RU-EE CN-KZ RU-CN CN-IN RU-BY CN-HK RU-AZ ' + + 'CN-BT CN-AF'