diff --git a/__snapshots__/test/snapshots.test.ts.js b/__snapshots__/test/snapshots.test.ts.js index 2498f17cd0..0ab84c0b2f 100644 --- a/__snapshots__/test/snapshots.test.ts.js +++ b/__snapshots__/test/snapshots.test.ts.js @@ -193,12 +193,12 @@ exports['Physics snapshots Barnes Hut defaults binary tree 1'] = { exports['Physics snapshots Force Atlas 2 defaults 2 connected nodes 1'] = { "1": { - "x": -11, - "y": 58 + "x": -14, + "y": 64 }, "2": { - "x": -8, - "y": -68 + "x": -10, + "y": -74 } } @@ -215,172 +215,172 @@ exports['Physics snapshots Force Atlas 2 defaults 2 disconnected nodes 1'] = { exports['Physics snapshots Force Atlas 2 defaults 2 doubly connected nodes 1'] = { "1": { - "x": -18, - "y": 64 + "x": -16, + "y": 71 }, "2": { - "x": 3, - "y": -65 + "x": 0, + "y": -73 } } exports['Physics snapshots Force Atlas 2 defaults 6 node circle 1'] = { "1": { - "x": -23, - "y": -151 + "x": -31, + "y": -191 }, "2": { - "x": -133, - "y": -48 + "x": -173, + "y": -61 }, "3": { - "x": -111, - "y": 102 + "x": -141, + "y": 129 }, "4": { - "x": 34, - "y": 147 + "x": 41, + "y": 187 }, "5": { - "x": 145, - "y": 45 + "x": 184, + "y": 58 }, "6": { - "x": 120, - "y": -104 + "x": 151, + "y": -131 } } exports['Physics snapshots Force Atlas 2 defaults 6 node complete graph 1'] = { "1": { - "x": -160, - "y": -106 + "x": -211, + "y": -144 }, "2": { - "x": -171, - "y": 94 + "x": -234, + "y": 117 }, "3": { - "x": 175, - "y": 84 + "x": 220, + "y": 137 }, "4": { - "x": 11, - "y": 197 + "x": -8, + "y": 264 }, "5": { - "x": 39, - "y": -151 + "x": 22, + "y": -247 }, "6": { - "x": 110, - "y": -109 + "x": 220, + "y": -116 } } exports['Physics snapshots Force Atlas 2 defaults 6 node star 1'] = { "1": { - "x": 12, - "y": -2 + "x": 9, + "y": -1 }, "2": { - "x": -146, - "y": -12 + "x": -196, + "y": -16 }, "3": { - "x": -40, - "y": -150 + "x": -51, + "y": -197 }, "4": { - "x": 154, - "y": -73 + "x": 193, + "y": -97 }, "5": { - "x": 121, - "y": 116 + "x": 155, + "y": 145 }, "6": { - "x": -48, - "y": 146 + "x": -68, + "y": 192 } } exports['Physics snapshots Force Atlas 2 defaults binary tree 1'] = { "1": { - "x": 52, - "y": 48 + "x": 73, + "y": 104 }, "2": { - "x": -97, - "y": 87 + "x": -129, + "y": 123 }, "3": { - "x": 206, - "y": 34 + "x": 269, + "y": 41 }, "4": { - "x": -113, - "y": -86 + "x": -166, + "y": -124 }, "5": { - "x": -135, - "y": 245 + "x": -202, + "y": 345 }, "6": { - "x": 211, - "y": -132 + "x": 327, + "y": -187 }, "7": { - "x": 199, - "y": 204 + "x": 312, + "y": 276 }, "8": { - "x": -278, - "y": -65 + "x": -396, + "y": -106 }, "9": { - "x": -26, - "y": -228 + "x": -61, + "y": -342 }, "10": { - "x": -279, - "y": 227 + "x": -382, + "y": 329 }, "11": { - "x": -49, - "y": 358 + "x": -86, + "y": 488 }, "12": { - "x": 356, - "y": -104 + "x": 501, + "y": -161 }, "13": { - "x": 213, - "y": -277 + "x": 296, + "y": -368 }, "14": { - "x": 340, - "y": 156 + "x": 484, + "y": 221 }, "15": { - "x": 134, - "y": 333 + "x": 215, + "y": 439 }, "16": { - "x": -304, - "y": -209 + "x": -437, + "y": -282 }, "17": { - "x": -360, - "y": 55 + "x": -499, + "y": 49 }, "18": { - "x": 63, - "y": -344 + "x": 81, + "y": -461 }, "19": { - "x": -138, - "y": -327 + "x": -182, + "y": -474 } } diff --git a/lib/network/modules/components/physics/FA2BasedRepulsionSolver.js b/lib/network/modules/components/physics/FA2BasedRepulsionSolver.js index 16fe1dbbfe..9aef81d36a 100644 --- a/lib/network/modules/components/physics/FA2BasedRepulsionSolver.js +++ b/lib/network/modules/components/physics/FA2BasedRepulsionSolver.js @@ -16,6 +16,57 @@ class ForceAtlas2BasedRepulsionSolver extends BarnesHutSolver { this._rng = Alea("FORCE ATLAS 2 BASED REPULSION SOLVER"); } + /** + * The repulsion "charge" of a node: its mass weighted by its degree + * (`edges.length + 1`). ForceAtlas2 scales repulsion by node degree so that + * highly interconnected hubs are pushed to the periphery. + * + * The degree must enter the force symmetrically. ForceAtlas2's repulsion is + * `k * (deg(a)+1) * (deg(b)+1) / distance` (Jacomy et al. 2014), i.e. it + * depends on the product of both degrees, so the force on `a` from `b` equals + * the force on `b` from `a` (Newton's third law). Applying the degree to only + * the receiving node — as was done previously — makes the pairwise force + * non-reciprocal and injects a net torque every tick. Velocity damping caps + * but never removes it, so a converged layout spins about its center of mass + * forever (issue #2193). Folding the degree into the mass used both to build + * the Barnes-Hut tree and to evaluate the force keeps it reciprocal. + * @param {Node} node + * @returns {number} the degree-weighted mass of the node + * @private + */ + _getMass(node) { + return node.options.mass * (node.edges.length + 1); + } + + /** + * Update the mass of a branch using the degree-weighted mass (see _getMass), + * so the accumulated branch mass and its center of mass are consistent with + * the degree-weighted repulsion applied in _calculateForces. + * @param {object} parentBranch + * @param {Node} node + * @private + */ + _updateBranchMass(parentBranch, node) { + const mass = this._getMass(node); + const centerOfMass = parentBranch.centerOfMass; + const totalMass = parentBranch.mass + mass; + const totalMassInv = 1 / totalMass; + + centerOfMass.x = centerOfMass.x * parentBranch.mass + node.x * mass; + centerOfMass.x *= totalMassInv; + + centerOfMass.y = centerOfMass.y * parentBranch.mass + node.y * mass; + centerOfMass.y *= totalMassInv; + + parentBranch.mass = totalMass; + const biggestSize = Math.max( + Math.max(node.height, node.radius), + node.width, + ); + parentBranch.maxWidth = + parentBranch.maxWidth < biggestSize ? biggestSize : parentBranch.maxWidth; + } + /** * Calculate the forces based on the distance. * @param {number} distance @@ -38,14 +89,15 @@ class ForceAtlas2BasedRepulsionSolver extends BarnesHutSolver { ); } - const degree = node.edges.length + 1; - // the dividing by the distance cubed instead of squared allows us to get the fx and fy components without sines and cosines - // it is shorthand for gravityforce with distance squared and fx = dx/distance * gravityForce + // dividing by the distance squared (rather than the cubed used by + // BarnesHutSolver) yields ForceAtlas2's linear 1/distance repulsion once the + // dx/dy components are factored in. Both parentBranch.mass and the node's + // mass are degree-weighted (see _getMass), so the repulsion is proportional + // to (degree+1) on both sides and stays reciprocal. const gravityForce = (this.options.gravitationalConstant * parentBranch.mass * - node.options.mass * - degree) / + this._getMass(node)) / Math.pow(distance, 2); const fx = dx * gravityForce; const fy = dy * gravityForce; diff --git a/test/physics/forceatlas2-rotation.test.ts b/test/physics/forceatlas2-rotation.test.ts new file mode 100644 index 0000000000..01e4c664ae --- /dev/null +++ b/test/physics/forceatlas2-rotation.test.ts @@ -0,0 +1,197 @@ +import { expect } from "chai"; +import Network from "../../lib/network/Network.js"; +import { canvasMockify } from "../canvas-mock.js"; + +/** + * Regression test for issue #2193: with the `forceAtlas2Based` solver an + * asymmetric graph converges to a stable *shape* but the whole layout then + * keeps rotating about its center of mass forever. + * + * Root cause: `FA2BasedRepulsionSolver` scaled the pairwise repulsion by the + * *receiving* node's degree only, so the repulsion between two nodes of + * different degree was non-reciprocal -- it violated Newton's third law. + * ForceAtlas2's repulsion is defined as `k * (deg(a)+1) * (deg(b)+1) / dist` + * (Jacomy et al. 2014): symmetric in the two degrees, hence reciprocal. The + * one-sided factor injects a small net torque every tick; velocity damping + * removes relative motion but not a rigid rotation, so the converged layout + * spins at a terminal rate. The fix folds `(degree+1)` into the Barnes-Hut + * effective mass (both when building the tree and when evaluating the force), + * which keeps the degree-weighted "hubs to the periphery" behaviour while + * making the force reciprocal again. The bug only shows on a + * degree-heterogeneous, asymmetric graph -- on a symmetric one the per-pair + * torques cancel, which is why it stayed undiagnosed. + * + * What this test measures: the *actual rotation of the node positions*, not a + * velocity proxy. It drives the engine with `minVelocity: 0` so the solver + * never reports "stabilized" (the regime in which the live network keeps + * ticking and the rotation is visible), converges the shape, then measures the + * best-fit rigid rotation angle of the whole configuration over a long window + * via the orthogonal Procrustes formula + * angle = atan2( Σ(x0*y1 - y0*x1), Σ(x0*x1 + y0*y1) ) + * where (x0,y0) and (x1,y1) are node positions relative to the centroid at the + * start and end of the window. + * + * On the buggy code the layout rotates ~122 deg over the 20000-step window and + * never settles; with the fix it rotates < 0.5 deg and settles. The threshold + * below (5 deg) sits more than an order of magnitude clear of both. + */ +describe("ForceAtlas2 rotation (issue #2193)", function (): void { + beforeEach(function () { + this.clearJSDOM = canvasMockify("
"); + this.container = document.getElementById("mynetwork"); + }); + + afterEach(function () { + this.clearJSDOM(); + + delete this.clearJSDOM; + delete this.container; + }); + + it("a converged asymmetric layout does not keep rotating", function (): void { + this.timeout(120000); + + // --- Build a deterministic, asymmetric, degree-heterogeneous graph. ------ + // A symmetric graph would NOT exhibit the bug (its per-pair torques cancel), + // so the heterogeneity of the hub degrees below is essential. Initial + // positions come from a seeded PRNG so the run is fully reproducible. + const rng = mulberry32(0xc0ffee); + const nodes: { id: number; x: number; y: number }[] = []; + const edges: { from: number; to: number }[] = []; + + const place = (id: number): void => { + nodes.push({ id, x: (rng() - 0.5) * 600, y: (rng() - 0.5) * 600 }); + }; + + const hubs = [0, 1, 2, 3, 4]; + hubs.forEach((h): void => place(h)); + // Connect the hubs in a cycle. + for (let i = 0; i < hubs.length; i++) { + edges.push({ from: hubs[i], to: hubs[(i + 1) % hubs.length] }); + } + // Give each hub a very different number of leaves -> heterogeneous degrees. + let nextId = hubs.length; + const leavesPerHub = [40, 20, 10, 4, 1]; + hubs.forEach((hub, idx): void => { + for (let k = 0; k < leavesPerHub[idx]; k++) { + const id = nextId++; + place(id); + edges.push({ from: hub, to: id }); + } + }); + + const network = new Network( + this.container, + { nodes, edges }, + { + physics: { + solver: "forceAtlas2Based", + // Constant timestep -> deterministic, comparable measurements. Drive + // the engine ourselves below instead of relying on the asynchronous, + // render-loop-based stabilization. + adaptiveTimestep: false, + stabilization: { enabled: false }, + // Never report "stabilized" so physicsTick keeps integrating; this is + // the regime in which the perpetual rotation is observable (a settled + // engine would simply freeze and hide it). + minVelocity: 0, + }, + }, + ); + + const engine = network.physics; + engine.updatePhysicsData(); + + // Converge the shape. + for (let i = 0; i < 6000; i++) { + engine.physicsTick(); + } + + // Snapshot positions relative to the centroid, advance a long window, then + // measure the net rigid rotation between the two snapshots. + const before = relativePositions(network); + const windowSteps = 20000; + for (let i = 0; i < windowSteps; i++) { + engine.physicsTick(); + } + const after = relativePositions(network); + + const rotationDeg = Math.abs( + (bestFitRotation(before, after) * 180) / Math.PI, + ); + + // Stop the render loop before the JSDOM is torn down (the canvas mock + // shims requestAnimationFrame onto setTimeout, which would otherwise fire + // after teardown and throw "window is not defined"). + network.destroy(); + + // Buggy code: ~122 deg (and still climbing -- it never settles). + // Fixed code: < 0.5 deg. The threshold sits well clear of both. + expect(rotationDeg).to.be.below(5); + }); +}); + +/** + * Node positions relative to the configuration's centroid. + * @param network - the vis Network to inspect + * @returns parallel arrays of x and y offsets from the centroid + */ +function relativePositions(network: any): { x: number[]; y: number[] } { + const nodes = network.body.nodes; + const ids = network.physics.physicsBody.physicsNodeIndices; + + let cx = 0; + let cy = 0; + for (const id of ids) { + cx += nodes[id].x; + cy += nodes[id].y; + } + cx /= ids.length; + cy /= ids.length; + + const x: number[] = []; + const y: number[] = []; + for (const id of ids) { + x.push(nodes[id].x - cx); + y.push(nodes[id].y - cy); + } + return { x, y }; +} + +/** + * Best-fit rigid rotation angle taking configuration `a` onto configuration `b` + * (orthogonal Procrustes, 2D): angle = atan2(Σ a×b, Σ a·b). Both configurations + * must be centered (see relativePositions) and index-aligned. + * @param a - start positions relative to centroid + * @param b - end positions relative to centroid + * @returns the signed best-fit rotation angle in radians + */ +function bestFitRotation( + a: { x: number[]; y: number[] }, + b: { x: number[]; y: number[] }, +): number { + let cross = 0; // Σ (ax*by - ay*bx) + let dot = 0; // Σ (ax*bx + ay*by) + for (let i = 0; i < a.x.length; i++) { + cross += a.x[i] * b.y[i] - a.y[i] * b.x[i]; + dot += a.x[i] * b.x[i] + a.y[i] * b.y[i]; + } + return Math.atan2(cross, dot); +} + +/** + * Small seeded PRNG (mulberry32) so initial node positions are deterministic + * and the test is fully reproducible. + * @param seed - 32-bit seed + * @returns a function returning the next pseudo-random number in [0, 1) + */ +function mulberry32(seed: number): () => number { + let a = seed >>> 0; + return function (): number { + a |= 0; + a = (a + 0x6d2b79f5) | 0; + let t = Math.imul(a ^ (a >>> 15), 1 | a); + t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t; + return ((t ^ (t >>> 14)) >>> 0) / 4294967296; + }; +}