From 2256af7448bcf5490a4d3c1d2fbb089cd1518c9f Mon Sep 17 00:00:00 2001 From: Pitchaya Boonsarngsuk <2285135b@student.gla.ac.uk> Date: Thu, 22 Mar 2018 16:02:45 +0000 Subject: [PATCH] =?UTF-8?q?=E0=B9=81=E0=B8=81=E0=B9=89=20coding=20style?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/barnesHut.js | 310 ++++++----- src/hybridSimulation.js | 406 +++++++------- src/interpolation/helpers.js | 7 +- src/interpolation/interpBruteForce.js | 2 +- src/interpolation/interpCommon.js | 18 +- src/link.js | 216 ++++---- src/neighbourSampling.js | 2 +- src/t-sne.js | 752 +++++++++++++------------- src/xy.js | 26 +- 9 files changed, 866 insertions(+), 873 deletions(-) diff --git a/src/barnesHut.js b/src/barnesHut.js index fed7583..24f3a5b 100644 --- a/src/barnesHut.js +++ b/src/barnesHut.js @@ -1,158 +1,152 @@ -import constant from "./constant"; -import jiggle from "./jiggle"; -import {x, y} from "./xy"; -import {quadtree} from "d3-quadtree"; - -/** - * The refinement of the existing Barnes-Hut implementation in D3 - * to fit the use case of the project. Previously the algorithm stored - * strength as internal node, now the random child is stored as internal - * node and the force calculations are done between the node and that internal - * object if they are sufficiently far away. - * The check to see if the nodes are far away was also changed to the one described in original Barnes-Hut paper. - * @return {force} calculated forces. - */ -export default function() { - var nodes, - node, - alpha, - distance = constant(300), - theta = 0.5; - - /** - * Constructs a quadtree at every iteration and apply the forces by visiting - * each node in a tree. - * @param {number} _ - controls the stopping of the - * particle simulations. - */ - function force(_) { - var i, n = nodes.length, tree = quadtree(nodes, x, y).visitAfter(accumulate); - for (alpha = _, i = 0; i < n; ++i) { - node = nodes[i], tree.visit(apply); - } - } - - /** - * Function used during the tree construction to fill out the nodes with - * correct data. Internal nodes acquire the random child while the leaf - * nodes accumulate forces from coincident quadrants. - * @param {quadrant} quad - node representing the quadrant in quadtree. - */ - function accumulate(quad) { - var q, d, children = []; - - // For internal nodes, accumulate forces from child quadrants. - if (quad.length) { - for (var i = 0; i < 4; ++i) { - if ((q = quad[i]) && (d = q.data)) { - children.push(d); - } - } - // Choose a random child. - quad.data = children[Math.floor(Math.random() * children.length)]; - quad.x = quad.data.x; - quad.y = quad.data.y; - } - - // For leaf nodes, accumulate forces from coincident quadrants. - else { - q = quad; - q.x = q.data.x; - q.y = q.data.y; - } - } - - /** - * Function that applies the forces for each node. If the objects are - * far away, the approximation is made. Otherwise, forces are calculated - * directly between the nodes. - * @param {quadrant} quad - node representing the quadrant in quadtree. - * @param {number} x1 - lower x bound of the node. - * @param {number} _ - lower y bound of the node. - * @param {number} x2 - upper x bound of the node. - * @return {boolean} - true if the approximation was applied. - */ - function apply(quad, x1, _, x2) { - - var x = quad.data.x + quad.data.vx - node.x - node.vx, - y = quad.data.y + quad.data.vy - node.y - node.vy, - w = x2 - x1, - l = Math.sqrt(x * x + y * y); - - // Apply the Barnes-Hut approximation if possible. - // Limit forces for very close nodes; randomize direction if coincident. - if (w / l < theta) { - if (x === 0) x = jiggle(), l += x * x; - if (y === 0) y = jiggle(), l += y * y; - if (quad.data) { - l = (l - +distance(node, quad.data)) / l * alpha; - x *= l, y *= l; - quad.data.vx -= x; - quad.data.vy -= y; - node.vx += x; - node.vy += y; - } - return true; - } - - // Otherwise, process points directly. - else if (quad.length) return; - - // Limit forces for very close nodes; randomize direction if coincident. - if (quad.data !== node || quad.next) { - if (x === 0) x = jiggle(), l += x * x; - if (y === 0) y = jiggle(), l += y * y; - } - - do if (quad.data !== node) { - l = (l - +distance(node, quad.data)) / l * alpha; - x *= l, y *= l; - quad.data.vx -= x; - quad.data.vy -= y; - node.vx += x; - node.vy += y; - } while (quad = quad.next); - } - - /** - * Calculates the stress. Basically, it computes the difference between - * high dimensional distance and real distance. The lower the stress is, - * the better layout. - * @return {number} - stress of the layout. - */ - function getStress() { - var totalDiffSq = 0, totalHighDistSq = 0; - for (var i = 0, source, target, realDist, highDist; i < nodes.length; i++) { - for (var j = 0; j < nodes.length; j++) { - if (i !== j) { - source = nodes[i], target = nodes[j]; - realDist = Math.hypot(target.x-source.x, target.y-source.y); - highDist = +distance(nodes[i], nodes[j]); - totalDiffSq += Math.pow(realDist-highDist, 2); - totalHighDistSq += highDist * highDist; - } - } - } - return Math.sqrt(totalDiffSq/totalHighDistSq); - } - - // API for initializing the algorithm, setting parameters and querying - // metrics. - force.initialize = function(_) { - nodes = _; - }; - - force.distance = function(_) { - return arguments.length ? (distance = typeof _ === "function" ? _ : constant(+_), force) : distance; - }; - - force.theta = function(_) { - return arguments.length ? (theta = _, force) : theta; - }; - - force.stress = function() { - return getStress(); - }; - - return force; -} +import constant from "./constant"; +import jiggle from "./jiggle"; +import {x, y} from "./xy"; +import {quadtree} from "d3-quadtree"; + +/** + * The refinement of the existing Barnes-Hut implementation in D3 + * to fit the use case of the project. Previously the algorithm stored + * strength as internal node, now the random child is stored as internal + * node and the force calculations are done between the node and that internal + * object if they are sufficiently far away. + * The check to see if the nodes are far away was also changed to the one described in original Barnes-Hut paper. + * @return {force} calculated forces. + */ +export default function() { + var nodes, + node, + alpha, + distance = constant(300), + theta = 0.5; + + /** + * Constructs a quadtree at every iteration and apply the forces by visiting + * each node in a tree. + * @param {number} _ - controls the stopping of the + * particle simulations. + */ + function force(_) { + var i, n = nodes.length, tree = quadtree(nodes, x, y).visitAfter(accumulate); + for (alpha = _, i = 0; i < n; ++i) { + node = nodes[i], tree.visit(apply); + } + } + + /** + * Function used during the tree construction to fill out the nodes with + * correct data. Internal nodes acquire the random child while the leaf + * nodes accumulate forces from coincident quadrants. + * @param {quadrant} quad - node representing the quadrant in quadtree. + */ + function accumulate(quad) { + var q, d, children = []; + + // For internal nodes, accumulate forces from child quadrants. + if (quad.length) { + for (var i = 0; i < 4; ++i) { + if ((q = quad[i]) && (d = q.data)) { + children.push(d); + } + } + // Choose a random child. + quad.data = children[Math.floor(Math.random() * children.length)]; + quad.x = quad.data.x; + quad.y = quad.data.y; + } else { // For leaf nodes, accumulate forces from coincident quadrants. + q = quad; + q.x = q.data.x; + q.y = q.data.y; + } + } + + /** + * Function that applies the forces for each node. If the objects are + * far away, the approximation is made. Otherwise, forces are calculated + * directly between the nodes. + * @param {quadrant} quad - node representing the quadrant in quadtree. + * @param {number} x1 - lower x bound of the node. + * @param {number} _ - lower y bound of the node. + * @param {number} x2 - upper x bound of the node. + * @return {boolean} - true if the approximation was applied. + */ + function apply(quad, x1, _, x2) { + + var x = quad.data.x + quad.data.vx - node.x - node.vx, + y = quad.data.y + quad.data.vy - node.y - node.vy, + w = x2 - x1, + l = Math.sqrt(x * x + y * y); + + // Apply the Barnes-Hut approximation if possible. + // Limit forces for very close nodes; randomize direction if coincident. + if (w / l < theta) { + if (x === 0) x = jiggle(), l += x * x; + if (y === 0) y = jiggle(), l += y * y; + if (quad.data) { + l = (l - +distance(node, quad.data)) / l * alpha; + x *= l, y *= l; + quad.data.vx -= x; + quad.data.vy -= y; + node.vx += x; + node.vy += y; + } + return true; + } else if (quad.length) return; // Otherwise, process points directly. + + // Limit forces for very close nodes; randomize direction if coincident. + if (quad.data !== node || quad.next) { + if (x === 0) x = jiggle(), l += x * x; + if (y === 0) y = jiggle(), l += y * y; + } + + do if (quad.data !== node) { + l = (l - +distance(node, quad.data)) / l * alpha; + x *= l, y *= l; + quad.data.vx -= x; + quad.data.vy -= y; + node.vx += x; + node.vy += y; + } while (quad = quad.next); + } + + /** + * Calculates the stress. Basically, it computes the difference between + * high dimensional distance and real distance. The lower the stress is, + * the better layout. + * @return {number} - stress of the layout. + */ + function getStress() { + var totalDiffSq = 0, totalHighDistSq = 0; + for (var i = 0, source, target, realDist, highDist; i < nodes.length; i++) { + for (var j = 0; j < nodes.length; j++) { + if (i !== j) { + source = nodes[i], target = nodes[j]; + realDist = Math.hypot(target.x-source.x, target.y-source.y); + highDist = +distance(nodes[i], nodes[j]); + totalDiffSq += Math.pow(realDist-highDist, 2); + totalHighDistSq += highDist * highDist; + } + } + } + return Math.sqrt(totalDiffSq/totalHighDistSq); + } + + // API for initializing the algorithm, setting parameters and querying + // metrics. + force.initialize = function(_) { + nodes = _; + }; + + force.distance = function(_) { + return arguments.length ? (distance = typeof _ === "function" ? _ : constant(+_), force) : distance; + }; + + force.theta = function(_) { + return arguments.length ? (theta = _, force) : theta; + }; + + force.stress = function() { + return getStress(); + }; + + return force; +} diff --git a/src/hybridSimulation.js b/src/hybridSimulation.js index 77d6258..12ab577 100644 --- a/src/hybridSimulation.js +++ b/src/hybridSimulation.js @@ -1,203 +1,203 @@ -import {dispatch} from "d3-dispatch"; -import constant from "./constant"; -import interpBruteForce from "./interpolation/interpBruteForce"; -import interpolationPivots from "./interpolation/interpolationPivots"; -import {takeSampleFrom} from "./interpolation/helpers"; - -/** - * An implementation of Chalmers, Morrison, and Ross' 2002 hybrid layout - * algorithm with an option to use the 2003 pivot-based near neighbour searching - * method. - * It performs the 1996 neighbour sampling spring simulation model, on a - * "sample set", sqrt(n) samples of the data. - * Other data points are then interpolated into the model. - * Finally, another spring simulation may be performed on the entire dataset to - * clean up the model. - * @param {object} sim - D3 Simulation object - * @param {object} forceS - Pre-configured D3 force object for the sample set. - The ending handler will be re-configured. - Neighbour sampling force is expected, but other D3 - forces may also work. - * @param {object} forceF - Pre-configured D3 force object for the simultion ran - on the entire dataset at the end. - Neighbour sampling force is expected, but other D3 - forces may also work. - The force should not have any ending condition. - */ -export default function (sim, forceS, forceF) { - var - SAMPLE_ITERATIONS = 300, - FULL_ITERATIONS = 20, - interpDistanceFn, - NUM_PIVOTS = 0, - INTERP_FINE_ITS = 20, - sample = [], - remainder = [], - simulation = sim, - forceSample = forceS, - forceFull = forceF, - event = d3.dispatch("sampleTick", "fullTick", "startInterp", "end"), - initAlready = false, - nodes, - alreadyRanIterations, - hybrid; - - if(simulation != undefined) initSimulation(); - if(forceS != undefined || forceF != undefined) initForces(); - - // Performed on first run - function initialize() { - initAlready = true; - alreadyRanIterations = 0; - simulation - .on("tick", sampleTick) - .on("end", sampleEnded) - .nodes(sample) - .force("Sample force", forceSample); - console.log("Initialized Simulation for Hybrid"); - } - - function initForces(){ - if (forceSample.onStableVelo) { - forceSample.onStableVelo(sampleEnded); - } - - if (forceFull.onStableVelo) { - forceFull.onStableVelo(fullEnded); - } - - // Set default value for interpDistanceFn if not been specified yet - if(interpDistanceFn === undefined) { - if(forceFull.distance == 'function') - interpDistanceFn = forceFull.distance(); - else - interpDistanceFn = constant(300); - } - } - - function initSimulation(){ - nodes = simulation.nodes(); - simulation - .stop() - .alphaDecay(0) - .alpha(1) - - let sets = takeSampleFrom(nodes, Math.sqrt(nodes.length)); - sample = sets.sample; - remainder = sets.remainder; - } - - // Sample simulation ticked 1 frame, keep track of number of iterations here. - function sampleTick() { - event.call("sampleTick"); - if(alreadyRanIterations++ >= SAMPLE_ITERATIONS){ - sampleEnded(); - } - } - - // Full simulation ticked 1 frame, keep track of number of iterations here. - function fullTick() { - event.call("fullTick"); - if(alreadyRanIterations++ >= FULL_ITERATIONS){ - fullEnded(); - } - } - - function fullEnded() { - simulation.stop(); - initAlready = false; - simulation.force("Full force", null); - event.call("end"); - } - - function sampleEnded() { - simulation.stop(); - simulation.force("Sample force", null); - // Reset velocity of all nodes - for (let i=sample.length-1; i>=0; i--){ - sample[i].vx=0; - sample[i].vy=0; - } - - event.call("startInterp"); - if (NUM_PIVOTS>=1) { - interpolationPivots(sample, remainder, NUM_PIVOTS, interpDistanceFn, INTERP_FINE_ITS); - } else { - interpBruteForce(sample, remainder, interpDistanceFn, INTERP_FINE_ITS); - } - - event.call("fullTick"); - alreadyRanIterations = 0; - simulation - .on("tick", null) - .on("end", null) // The ending condition should be iterations count - .nodes(nodes); - - if (FULL_ITERATIONS<1 || forceF === undefined || forceF === null) { - event.call("end"); - return; - } - simulation - .on("tick", fullTick) - .force("Full force", forceFull) - .restart(); - } - - return hybrid = { - restart: function () { - if(!initAlready) initialize(); - simulation.restart(); - return hybrid; - }, - - stop: function () { - simulation.stop(); - return hybrid; - }, - - numPivots: function (_) { - return arguments.length ? (NUM_PIVOTS = +_, hybrid) : NUM_PIVOTS; - }, - - sampleIterations: function (_) { - return arguments.length ? (SAMPLE_ITERATIONS = +_, hybrid) : SAMPLE_ITERATIONS; - }, - - fullIterations: function (_) { - return arguments.length ? (FULL_ITERATIONS = +_, hybrid) : FULL_ITERATIONS; - }, - - interpFindTuneIts: function (_) { - return arguments.length ? (INTERP_FINE_ITS = +_, hybrid) : INTERP_FINE_ITS; - }, - - on: function (name, _) { - return arguments.length > 1 ? (event.on(name, _), hybrid) : event.on(name); - }, - - subSet: function (_) { - return arguments.length ? (sample = _, hybrid) : sample; - }, - - nonSubSet: function (_) { - return arguments.length ? (remainder = _, hybrid) : remainder; - }, - - interpDistanceFn: function (_) { - return arguments.length ? (interpDistanceFn = typeof _ === "function" ? _ : constant(+_), hybrid) : interpDistanceFn; - }, - - simulation: function (_) { - return arguments.length ? (initAlready = false, simulation = _, hybrid) : simulation; - }, - - forceSample: function (_) { - return arguments.length ? (forceSample = _, initForces(), hybrid) : forceSample; - }, - - forceFull: function (_) { - return arguments.length ? (forceFull = _, initForces(), hybrid) : forceFull; - }, - - }; -} +import {dispatch} from "d3-dispatch"; +import constant from "./constant"; +import interpBruteForce from "./interpolation/interpBruteForce"; +import interpolationPivots from "./interpolation/interpolationPivots"; +import {takeSampleFrom} from "./interpolation/helpers"; + +/** + * An implementation of Chalmers, Morrison, and Ross' 2002 hybrid layout + * algorithm with an option to use the 2003 pivot-based near neighbour searching + * method. + * It performs the 1996 neighbour sampling spring simulation model, on a + * "sample set", sqrt(n) samples of the data. + * Other data points are then interpolated into the model. + * Finally, another spring simulation may be performed on the entire dataset to + * clean up the model. + * @param {object} sim - D3 Simulation object + * @param {object} forceS - Pre-configured D3 force object for the sample set. + The ending handler will be re-configured. + Neighbour sampling force is expected, but other D3 + forces may also work. + * @param {object} forceF - Pre-configured D3 force object for the simultion ran + on the entire dataset at the end. + Neighbour sampling force is expected, but other D3 + forces may also work. + The force should not have any ending condition. + */ +export default function (sim, forceS, forceF) { + var + SAMPLE_ITERATIONS = 300, + FULL_ITERATIONS = 20, + interpDistanceFn, + NUM_PIVOTS = 0, + INTERP_FINE_ITS = 20, + sample = [], + remainder = [], + simulation = sim, + forceSample = forceS, + forceFull = forceF, + event = d3.dispatch("sampleTick", "fullTick", "startInterp", "end"), + initAlready = false, + nodes, + alreadyRanIterations, + hybrid; + + if(simulation != undefined) initSimulation(); + if(forceS != undefined || forceF != undefined) initForces(); + + // Performed on first run + function initialize() { + initAlready = true; + alreadyRanIterations = 0; + simulation + .on("tick", sampleTick) + .on("end", sampleEnded) + .nodes(sample) + .force("Sample force", forceSample); + console.log("Initialized Simulation for Hybrid"); + } + + function initForces(){ + if (forceSample.onStableVelo) { + forceSample.onStableVelo(sampleEnded); + } + + if (forceFull.onStableVelo) { + forceFull.onStableVelo(fullEnded); + } + + // Set default value for interpDistanceFn if not been specified yet + if(interpDistanceFn === undefined) { + if(forceFull.distance == 'function') + interpDistanceFn = forceFull.distance(); + else + interpDistanceFn = constant(300); + } + } + + function initSimulation(){ + nodes = simulation.nodes(); + simulation + .stop() + .alphaDecay(0) + .alpha(1) + + let sets = takeSampleFrom(nodes, Math.sqrt(nodes.length)); + sample = sets.sample; + remainder = sets.remainder; + } + + // Sample simulation ticked 1 frame, keep track of number of iterations here. + function sampleTick() { + event.call("sampleTick"); + if(alreadyRanIterations++ >= SAMPLE_ITERATIONS){ + sampleEnded(); + } + } + + // Full simulation ticked 1 frame, keep track of number of iterations here. + function fullTick() { + event.call("fullTick"); + if(alreadyRanIterations++ >= FULL_ITERATIONS){ + fullEnded(); + } + } + + function fullEnded() { + simulation.stop(); + initAlready = false; + simulation.force("Full force", null); + event.call("end"); + } + + function sampleEnded() { + simulation.stop(); + simulation.force("Sample force", null); + // Reset velocity of all nodes + for (let i=sample.length-1; i>=0; i--){ + sample[i].vx=0; + sample[i].vy=0; + } + + event.call("startInterp"); + if (NUM_PIVOTS>=1) { + interpolationPivots(sample, remainder, NUM_PIVOTS, interpDistanceFn, INTERP_FINE_ITS); + } else { + interpBruteForce(sample, remainder, interpDistanceFn, INTERP_FINE_ITS); + } + + event.call("fullTick"); + alreadyRanIterations = 0; + simulation + .on("tick", null) + .on("end", null) // The ending condition should be iterations count + .nodes(nodes); + + if (FULL_ITERATIONS<1 || forceF === undefined || forceF === null) { + event.call("end"); + return; + } + simulation + .on("tick", fullTick) + .force("Full force", forceFull) + .restart(); + } + + return hybrid = { + restart: function () { + if(!initAlready) initialize(); + simulation.restart(); + return hybrid; + }, + + stop: function () { + simulation.stop(); + return hybrid; + }, + + numPivots: function (_) { + return arguments.length ? (NUM_PIVOTS = +_, hybrid) : NUM_PIVOTS; + }, + + sampleIterations: function (_) { + return arguments.length ? (SAMPLE_ITERATIONS = +_, hybrid) : SAMPLE_ITERATIONS; + }, + + fullIterations: function (_) { + return arguments.length ? (FULL_ITERATIONS = +_, hybrid) : FULL_ITERATIONS; + }, + + interpFindTuneIts: function (_) { + return arguments.length ? (INTERP_FINE_ITS = +_, hybrid) : INTERP_FINE_ITS; + }, + + on: function (name, _) { + return arguments.length > 1 ? (event.on(name, _), hybrid) : event.on(name); + }, + + subSet: function (_) { + return arguments.length ? (sample = _, hybrid) : sample; + }, + + nonSubSet: function (_) { + return arguments.length ? (remainder = _, hybrid) : remainder; + }, + + interpDistanceFn: function (_) { + return arguments.length ? (interpDistanceFn = typeof _ === "function" ? _ : constant(+_), hybrid) : interpDistanceFn; + }, + + simulation: function (_) { + return arguments.length ? (initAlready = false, simulation = _, hybrid) : simulation; + }, + + forceSample: function (_) { + return arguments.length ? (forceSample = _, initForces(), hybrid) : forceSample; + }, + + forceFull: function (_) { + return arguments.length ? (forceFull = _, initForces(), hybrid) : forceFull; + }, + + }; +} diff --git a/src/interpolation/helpers.js b/src/interpolation/helpers.js index 4d97a4e..6e2f5b2 100644 --- a/src/interpolation/helpers.js +++ b/src/interpolation/helpers.js @@ -10,8 +10,8 @@ */ export function takeSampleFrom(sourceList, amount) { let randElements = [], - max = sourceList.length, - swap = false; + max = sourceList.length, + swap = false; if (amount >= max) { return {sample: sourceList, remainder: {}}; @@ -40,8 +40,7 @@ export function takeSampleFrom(sourceList, amount) { sample: remainder, remainder: randElements }; - } - else { + } else { return { sample: randElements, remainder: remainder diff --git a/src/interpolation/interpBruteForce.js b/src/interpolation/interpBruteForce.js index 14231ed..c4acd92 100644 --- a/src/interpolation/interpBruteForce.js +++ b/src/interpolation/interpBruteForce.js @@ -23,7 +23,7 @@ export default function(sampleSet, remainderSet, distanceFn, endingIts) { sampleSubset = takeSampleFrom(sampleSet, Math.sqrt(sampleSet.length)).sample, sampleSubsetDistanceCache = []; - // For each datapoint "node" to be interpolated + // For each datapoint "node" to be interpolated for (let i = remainderSet.length-1; i>=0; i--) { let node = remainderSet[i], diff --git a/src/interpolation/interpCommon.js b/src/interpolation/interpCommon.js index 4b5f480..860c311 100644 --- a/src/interpolation/interpCommon.js +++ b/src/interpolation/interpCommon.js @@ -36,7 +36,7 @@ export function placeNearToNearestNeighbour(node, nearNeighbour, radius, sampleS lowBound = 0.0, highBound = 0.0; - // Determine the closest quadrant + // Determine the closest quadrant if (dist0 == dist180) { if (dist90 > dist270) lowBound = highBound = 270; @@ -55,14 +55,12 @@ export function placeNearToNearestNeighbour(node, nearNeighbour, radius, sampleS lowBound = 90; highBound = 180; } + } else if (dist90 > dist270) { + lowBound = 270; + highBound = 360; } else { - if (dist90 > dist270) { - lowBound = 270; - highBound = 360; - } else { - lowBound = 0; - highBound = 90; - } + lowBound = 0; + highBound = 90; } // Determine the angle @@ -84,8 +82,8 @@ export function placeNearToNearestNeighbour(node, nearNeighbour, radius, sampleS function sumForcesToSample(node, samples, sampleCache) { let nodeVx = 0, - nodeVy = 0, - x, y, l, i, sample; + nodeVy = 0, + x, y, l, i, sample; for (i = samples.length-1; i >=0 ; i--) { sample = samples[i]; diff --git a/src/link.js b/src/link.js index 9bae6f3..3f046e8 100644 --- a/src/link.js +++ b/src/link.js @@ -1,108 +1,108 @@ -import constant from "./constant"; -import jiggle from "./jiggle"; - -/** - * Modified link force algorithm - * - simplify calculations for parameters locked for spring model - * - replace the use of links {} with loop. greatly reduce memory usage - * - removed other unused functions - * Alpha should be constant 1 for accurate simulation - */ -export default function() { - var dataSizeFactor, - distance = constant(30), - distances = [], - nodes, - stableVelocity = 0, - stableVeloHandler = null, - latestVelocityDiff = 0, - iterations = 1; - - function force(alpha) { - let n = nodes.length; - // Cache old velocity for comparison later - if (stableVeloHandler!==null && stableVelocity>=0) { - for (let i = n-1, node; i>=0; i--) { - node = nodes[i]; - node.oldvx = node.vx; - node.oldvy = node.vy; - } - } - - // Each iteration in a tick - for (var k = 0, source, target, i, j, x, y, l; k < iterations; ++k) { - // For each link - for (i = 1; i < n; i++) for (j = 0; j < i; j++) { - // jiggle so l won't be zero and divide by zero error after this - source = nodes[i]; - target = nodes[j]; - x = target.x + target.vx - source.x - source.vx || jiggle(); - y = target.y + target.vy - source.y - source.vy || jiggle(); - l = Math.sqrt(x * x + y * y); - l = (l - distances[i*(i-1)/2+j]) / l * dataSizeFactor * alpha; - x *= l, y *= l; - target.vx -= x; - target.vy -= y; - source.vx += x; - source.vy += y; - } - } - - // Calculate velocity changes, aka force applied. - if (stableVeloHandler!==null && stableVelocity>=0) { - let velocityDiff = 0; - for (let i = n-1, node; i>=0; i--) { - node = nodes[i]; - velocityDiff += Math.abs(Math.hypot(node.vx-node.oldvx, node.vy-node.oldvy)); - } - velocityDiff /= n; - latestVelocityDiff = velocityDiff; - - if(velocityDiff=0) { + for (let i = n-1, node; i>=0; i--) { + node = nodes[i]; + node.oldvx = node.vx; + node.oldvy = node.vy; + } + } + + // Each iteration in a tick + for (var k = 0, source, target, i, j, x, y, l; k < iterations; ++k) { + // For each link + for (i = 1; i < n; i++) for (j = 0; j < i; j++) { + // jiggle so l won't be zero and divide by zero error after this + source = nodes[i]; + target = nodes[j]; + x = target.x + target.vx - source.x - source.vx || jiggle(); + y = target.y + target.vy - source.y - source.vy || jiggle(); + l = Math.sqrt(x * x + y * y); + l = (l - distances[i*(i-1)/2+j]) / l * dataSizeFactor * alpha; + x *= l, y *= l; + target.vx -= x; + target.vy -= y; + source.vx += x; + source.vy += y; + } + } + + // Calculate velocity changes, aka force applied. + if (stableVeloHandler!==null && stableVelocity>=0) { + let velocityDiff = 0; + for (let i = n-1, node; i>=0; i--) { + node = nodes[i]; + velocityDiff += Math.abs(Math.hypot(node.vx-node.oldvx, node.vy-node.oldvy)); + } + velocityDiff /= n; + latestVelocityDiff = velocityDiff; + + if(velocityDiff 1) {return gaussRandom();} - return u * Math.sqrt(-2 * Math.log(r) / r); - } - - /** - * Return the normalized number. - * @return {number} normalized random number from Gaussian distribution. - */ - function randomN() { - return gaussRandom() * 1e-4; - } - - function sign(x) { - return x > 0 ? 1 : x < 0 ? -1 : 0; - } - - /** - * Create an array of length n filled with zeros. - * @param {number} n - length of array. - * @return {Float64Array} - array of zeros with length n. - */ - function zeros(n) { - if (typeof n === 'undefined' || isNaN(n)) { - return []; - } - return new Float64Array(n); // typed arrays are faster - } - - // Returns a 2d array of random numbers - /** - * Creates a 2d array filled with random numbers. - * @param {number} n - rows. - * @param {number} d - columns. - * @return {array} - 2d array - */ - function random2d(n, d) { - var x = []; - for (var i = 0; i < n; i++) { - var y = []; - for (var j = 0; j < d; j++) { - y.push(randomN()); - } - x.push(y); - } - return x; - } - - /** - * Compute the probability matrix using the provided data. - * @param {array} data - nodes. - * @param {number} perplexity - used to calculate entropy of distribution. - * @param {number} tol - limit for entropy difference. - * @return {2d array} - 2d matrix containing probabilities. - */ - function d2p(data, perplexity, tol) { - N = Math.floor(data.length); - var Htarget = Math.log(perplexity); // target entropy of distribution. - var P1 = zeros(N * N); // temporary probability matrix. - - var prow = zeros(N); // a temporary storage compartment. - for (var i = 0; i < N; i++) { - var betamin = -Infinity; - var betamax = Infinity; - var beta = 1; // initial value of precision. - var done = false; - var maxtries = 50; - - // Perform binary search to find a suitable precision beta - // so that the entropy of the distribution is appropriate. - var num = 0; - while (!done) { - // Compute entropy and kernel row with beta precision. - var psum = 0.0; - for (var j = 0; j < N; j++) { - var pj = Math.exp(-distance(data[i], data[j]) * beta); - // Ignore the diagonals - if (i === j) { - pj = 0; - } - prow[j] = pj; - psum += pj; - } - // Normalize p and compute entropy. - var Hhere = 0.0; - for (j = 0; j < N; j++) { - if (psum == 0) { - pj = 0; - } else { - pj = prow[j] / psum; - } - prow[j] = pj; - if (pj > 1e-7) { - Hhere -= pj * Math.log(pj); - } - } - - // Adjust beta based on result. - if (Hhere > Htarget) { - // Entropy was too high (distribution too diffuse) - // so we need to increase the precision for more peaky distribution. - betamin = beta; // move up the bounds. - if (betamax === Infinity) { - beta = beta * 2; - } else { - beta = (beta + betamax) / 2; - } - - } else { - // Converse case. Make distrubtion less peaky. - betamax = beta; - if (betamin === -Infinity) { - beta = beta / 2; - } else { - beta = (beta + betamin) / 2; - } - } - - // Stopping conditions: too many tries or got a good precision. - num++; - if (Math.abs(Hhere - Htarget) < tol || num >= maxtries) { - done = true; - } - } - - // Copy over the final prow to P1 at row i - for (j = 0; j < N; j++) { - P1[i * N + j] = prow[j]; - } - - } - - // Symmetrize P and normalize it to sum to 1 over all ij - var Pout = zeros(N * N); - var N2 = N * 2; - for (i = 0; i < N; i++) { - for (j = 0; j < N; j++) { - Pout[i * N + j] = Math.max((P1[i * N + j] + P1[j * N + i]) / N2, 1e-100); - } - } - return Pout; - } - - /** - * Initialize a starting (random) solution. - */ - function initSolution() { - Y = random2d(N, dim); - // Step gains to accelerate progress in unchanging directions. - gains = random2d(N, dim, 1.0); - // Momentum accumulator. - ystep = random2d(N, dim, 0.0); - iteration = 0; - } - - /** - * @return {2d array} the solution. - */ - function getSolution() { - return Y; - } - - /** - * Do a single step (iteration) for the layout. - * @return {number} the current cost. - */ - function step() { - iteration += 1; - - var cg = costGrad(Y); // Evaluate gradient. - var cost = cg.cost; - var grad = cg.grad; - - // Perform gradient step. - var ymean = zeros(dim); - for (var i = 0; i < N; i++) { - for (var d = 0; d < dim; d++) { - var gid = grad[i][d]; - var sid = ystep[i][d]; - var gainid = gains[i][d]; - - // Compute gain update. - var newgain = sign(gid) === sign(sid) ? gainid * 0.8 : gainid + 0.2; - if (newgain < 0.01) { - newgain = 0.01; - } - gains[i][d] = newgain; - - // Compute momentum step direction. - var momval = iteration < 250 ? 0.5 : 0.8; - var newsid = momval * sid - learningRate * newgain * grad[i][d]; - ystep[i][d] = newsid; - - // Do the step. - Y[i][d] += newsid; - - // Accumulate mean so that we can center later. - ymean[d] += Y[i][d]; - } - } - - // Reproject Y to have the zero mean. - for (i = 0; i < N; i++) { - for (d = 0; d < dim; d++) { - Y[i][d] -= ymean[d] / N; - } - } - return cost; - } - - /** - * Calculate the cost and the gradient. - * @param {2d array} Y - the current solution to evaluate. - * @return {object} that contains a cost and a gradient. - */ - function costGrad(Y) { - - var pmul = iteration < 100 ? 4 : 1; - - // Compute current Q distribution, unnormalized first. - var Qu = zeros(N * N); - var qsum = 0.0; - for (var i = 0; i < N; i++) { - for (var j = i + 1; j < N; j++) { - var dsum = 0.0; - for (var d = 0; d < dim; d++) { - var dhere = Y[i][d] - Y[j][d]; - dsum += dhere * dhere; - } - var qu = 1.0 / (1.0 + dsum); // Student t-distribution. - Qu[i * N + j] = qu; - Qu[j * N + i] = qu; - qsum += 2 * qu; - } - } - // Normalize Q distribution to sum to 1. - var NN = N * N; - var Q = zeros(NN); - for (var q = 0; q < NN; q++) { - Q[q] = Math.max(Qu[q] / qsum, 1e-100); - } - - var cost = 0.0; - var grad = []; - for (i = 0; i < N; i++) { - var gsum = new Array(dim); // Initialize gradiet for point i. - for (d = 0; d < dim; d++) { - gsum[d] = 0.0; - } - for (j = 0; j < N; j++) { - // Accumulate the cost. - cost += -P[i * N + j] * Math.log(Q[i * N + j]); - var premult = 4 * (pmul * P[i * N + j] - Q[i * N + j]) * Qu[i * N + j]; - for (d = 0; d < dim; d++) { - gsum[d] += premult * (Y[i][d] - Y[j][d]); - } - } - grad.push(gsum); - } - - return { - cost: cost, - grad: grad - }; - } - - /** - * Calculates the stress. Basically, it computes the difference between - * high dimensional distance and real distance. The lower the stress is, - * the better layout. - * @return {number} - stress of the layout. - */ - function getStress() { - var totalDiffSq = 0, - totalHighDistSq = 0; - for (var i = 0, source, target, realDist, highDist; i < nodes.length; i++) { - for (var j = 0; j < nodes.length; j++) { - if (i !== j) { - source = nodes[i], target = nodes[j]; - realDist = Math.hypot(target.x - source.x, target.y - source.y); - highDist = +distance(nodes[i], nodes[j]); - totalDiffSq += Math.pow(realDist - highDist, 2); - totalHighDistSq += highDist * highDist; - } - } - } - return Math.sqrt(totalDiffSq / totalHighDistSq); - } - - // API for initializing the algorithm, setting parameters and querying - // metrics. - force.initialize = function(_) { - nodes = _; - N = nodes.length; - // Initialize the probability matrix. - P = d2p(nodes, perplexity, 1e-4); - initSolution(); - }; - - force.id = function(_) { - return arguments.length ? (id = _, force) : id; - }; - - force.distance = function(_) { - return arguments.length ? (distance = typeof _ === "function" ? _ : constant(+_), force) : distance; - }; - - force.stress = function() { - return getStress(); - }; - - force.learningRate = function(_) { - return arguments.length ? (learningRate = +_, force) : learningRate; - }; - - force.perplexity = function(_) { - return arguments.length ? (perplexity = +_, force) : perplexity; - }; - - return force; -} +/* eslint-disable block-scoped-var */ +import constant from "./constant"; + +/** + * Set the node id accessor to the specified i. + * @param {node} d - node. + * @param {accessor} i - id accessor. + * @return {accessor} - node id accessor. + */ +function index(d, i) { + return i; +} + +/** + * t-SNE implementation in D3 by using the code existing in tsnejs + * (https://github.com/karpathy/tsnejs) to compute the solution. + */ +export default function() { + var id = index, + distance = constant(300), + nodes, + perplexity = 30, + learningRate = 10, + iteration = 0, + dim = 2, + N, // length of the nodes. + P, // probability matrix. + Y, // solution. + gains, + ystep; + + /** + * Make a step in t-SNE algorithm and set the velocities for the nodes + * to accumulate the values from solution. + */ + function force() { + // Make a step at each iteration. + step(); + var solution = getSolution(); + + // Set the velocity for each node using the solution. + for (var i = 0; i < nodes.length; i++) { + nodes[i].vx += solution[i][0]; + nodes[i].vy += solution[i][1]; + } + } + + /** + * Calculates the random number from Gaussian distribution. + * @return {number} random number. + */ + function gaussRandom() { + let u = 2 * Math.random() - 1; + let v = 2 * Math.random() - 1; + let r = u * u + v * v; + if (r == 0 || r > 1) { + return gaussRandom(); + } + return u * Math.sqrt(-2 * Math.log(r) / r); + } + + /** + * Return the normalized number. + * @return {number} normalized random number from Gaussian distribution. + */ + function randomN() { + return gaussRandom() * 1e-4; + } + + function sign(x) { + return x > 0 ? 1 : x < 0 ? -1 : 0; + } + + /** + * Create an array of length n filled with zeros. + * @param {number} n - length of array. + * @return {Float64Array} - array of zeros with length n. + */ + function zeros(n) { + if (typeof n === 'undefined' || isNaN(n)) { + return []; + } + return new Float64Array(n); // typed arrays are faster + } + + // Returns a 2d array of random numbers + /** + * Creates a 2d array filled with random numbers. + * @param {number} n - rows. + * @param {number} d - columns. + * @return {array} - 2d array + */ + function random2d(n, d) { + var x = []; + for (var i = 0; i < n; i++) { + var y = []; + for (var j = 0; j < d; j++) { + y.push(randomN()); + } + x.push(y); + } + return x; + } + + /** + * Compute the probability matrix using the provided data. + * @param {array} data - nodes. + * @param {number} perplexity - used to calculate entropy of distribution. + * @param {number} tol - limit for entropy difference. + * @return {2d array} - 2d matrix containing probabilities. + */ + function d2p(data, perplexity, tol) { + N = Math.floor(data.length); + var Htarget = Math.log(perplexity); // target entropy of distribution. + var P1 = zeros(N * N); // temporary probability matrix. + + var prow = zeros(N); // a temporary storage compartment. + for (var i = 0; i < N; i++) { + var betamin = -Infinity; + var betamax = Infinity; + var beta = 1; // initial value of precision. + var done = false; + var maxtries = 50; + + // Perform binary search to find a suitable precision beta + // so that the entropy of the distribution is appropriate. + var num = 0; + while (!done) { + // Compute entropy and kernel row with beta precision. + var psum = 0.0; + for (var j = 0; j < N; j++) { + var pj = Math.exp(-distance(data[i], data[j]) * beta); + // Ignore the diagonals + if (i === j) { + pj = 0; + } + prow[j] = pj; + psum += pj; + } + // Normalize p and compute entropy. + var Hhere = 0.0; + for (j = 0; j < N; j++) { + if (psum == 0) { + pj = 0; + } else { + pj = prow[j] / psum; + } + prow[j] = pj; + if (pj > 1e-7) { + Hhere -= pj * Math.log(pj); + } + } + + // Adjust beta based on result. + if (Hhere > Htarget) { + // Entropy was too high (distribution too diffuse) + // so we need to increase the precision for more peaky distribution. + betamin = beta; // move up the bounds. + if (betamax === Infinity) { + beta = beta * 2; + } else { + beta = (beta + betamax) / 2; + } + + } else { + // Converse case. Make distrubtion less peaky. + betamax = beta; + if (betamin === -Infinity) { + beta = beta / 2; + } else { + beta = (beta + betamin) / 2; + } + } + + // Stopping conditions: too many tries or got a good precision. + num++; + if (Math.abs(Hhere - Htarget) < tol || num >= maxtries) { + done = true; + } + } + + // Copy over the final prow to P1 at row i + for (j = 0; j < N; j++) { + P1[i * N + j] = prow[j]; + } + + } + + // Symmetrize P and normalize it to sum to 1 over all ij + var Pout = zeros(N * N); + var N2 = N * 2; + for (i = 0; i < N; i++) { + for (j = 0; j < N; j++) { + Pout[i * N + j] = Math.max((P1[i * N + j] + P1[j * N + i]) / N2, 1e-100); + } + } + return Pout; + } + + /** + * Initialize a starting (random) solution. + */ + function initSolution() { + Y = random2d(N, dim); + // Step gains to accelerate progress in unchanging directions. + gains = random2d(N, dim, 1.0); + // Momentum accumulator. + ystep = random2d(N, dim, 0.0); + iteration = 0; + } + + /** + * @return {2d array} the solution. + */ + function getSolution() { + return Y; + } + + /** + * Do a single step (iteration) for the layout. + * @return {number} the current cost. + */ + function step() { + iteration += 1; + + var cg = costGrad(Y); // Evaluate gradient. + var cost = cg.cost; + var grad = cg.grad; + + // Perform gradient step. + var ymean = zeros(dim); + for (var i = 0; i < N; i++) { + for (var d = 0; d < dim; d++) { + var gid = grad[i][d]; + var sid = ystep[i][d]; + var gainid = gains[i][d]; + + // Compute gain update. + var newgain = sign(gid) === sign(sid) ? gainid * 0.8 : gainid + 0.2; + if (newgain < 0.01) { + newgain = 0.01; + } + gains[i][d] = newgain; + + // Compute momentum step direction. + var momval = iteration < 250 ? 0.5 : 0.8; + var newsid = momval * sid - learningRate * newgain * grad[i][d]; + ystep[i][d] = newsid; + + // Do the step. + Y[i][d] += newsid; + + // Accumulate mean so that we can center later. + ymean[d] += Y[i][d]; + } + } + + // Reproject Y to have the zero mean. + for (i = 0; i < N; i++) { + for (d = 0; d < dim; d++) { + Y[i][d] -= ymean[d] / N; + } + } + return cost; + } + + /** + * Calculate the cost and the gradient. + * @param {2d array} Y - the current solution to evaluate. + * @return {object} that contains a cost and a gradient. + */ + function costGrad(Y) { + + var pmul = iteration < 100 ? 4 : 1; + + // Compute current Q distribution, unnormalized first. + var Qu = zeros(N * N); + var qsum = 0.0; + for (var i = 0; i < N; i++) { + for (var j = i + 1; j < N; j++) { + var dsum = 0.0; + for (var d = 0; d < dim; d++) { + var dhere = Y[i][d] - Y[j][d]; + dsum += dhere * dhere; + } + var qu = 1.0 / (1.0 + dsum); // Student t-distribution. + Qu[i * N + j] = qu; + Qu[j * N + i] = qu; + qsum += 2 * qu; + } + } + // Normalize Q distribution to sum to 1. + var NN = N * N; + var Q = zeros(NN); + for (var q = 0; q < NN; q++) { + Q[q] = Math.max(Qu[q] / qsum, 1e-100); + } + + var cost = 0.0; + var grad = []; + for (i = 0; i < N; i++) { + var gsum = new Array(dim); // Initialize gradiet for point i. + for (d = 0; d < dim; d++) { + gsum[d] = 0.0; + } + for (j = 0; j < N; j++) { + // Accumulate the cost. + cost += -P[i * N + j] * Math.log(Q[i * N + j]); + var premult = 4 * (pmul * P[i * N + j] - Q[i * N + j]) * Qu[i * N + j]; + for (d = 0; d < dim; d++) { + gsum[d] += premult * (Y[i][d] - Y[j][d]); + } + } + grad.push(gsum); + } + + return { + cost: cost, + grad: grad + }; + } + + /** + * Calculates the stress. Basically, it computes the difference between + * high dimensional distance and real distance. The lower the stress is, + * the better layout. + * @return {number} - stress of the layout. + */ + function getStress() { + var totalDiffSq = 0, + totalHighDistSq = 0; + for (var i = 0, source, target, realDist, highDist; i < nodes.length; i++) { + for (var j = 0; j < nodes.length; j++) { + if (i !== j) { + source = nodes[i], target = nodes[j]; + realDist = Math.hypot(target.x - source.x, target.y - source.y); + highDist = +distance(nodes[i], nodes[j]); + totalDiffSq += Math.pow(realDist - highDist, 2); + totalHighDistSq += highDist * highDist; + } + } + } + return Math.sqrt(totalDiffSq / totalHighDistSq); + } + + // API for initializing the algorithm, setting parameters and querying + // metrics. + force.initialize = function(_) { + nodes = _; + N = nodes.length; + // Initialize the probability matrix. + P = d2p(nodes, perplexity, 1e-4); + initSolution(); + }; + + force.id = function(_) { + return arguments.length ? (id = _, force) : id; + }; + + force.distance = function(_) { + return arguments.length ? (distance = typeof _ === "function" ? _ : constant(+_), force) : distance; + }; + + force.stress = function() { + return getStress(); + }; + + force.learningRate = function(_) { + return arguments.length ? (learningRate = +_, force) : learningRate; + }; + + force.perplexity = function(_) { + return arguments.length ? (perplexity = +_, force) : perplexity; + }; + + return force; +} diff --git a/src/xy.js b/src/xy.js index 137b41d..647cc76 100644 --- a/src/xy.js +++ b/src/xy.js @@ -1,13 +1,13 @@ -/** - * @return x value of a node - */ -export function x(d) { - return d.x; -} - -/** - * @return y value of a node - */ -export function y(d) { - return d.y; -} \ No newline at end of file +/** + * @return x value of a node + */ +export function x(d) { + return d.x; +} + +/** + * @return y value of a node + */ +export function y(d) { + return d.y; +}