31:T6cb,md`# Point Neighborhood Lookup Grid ## Implementations for ${tex`\R^2`} and ${tex`\R^n`} **${compileText([ "For a point a set of points ", tex`P \subset \R^2`, ", find the subset ", tex`S \subset P`, " of points within a radius of influence ", tex`r`, " of any given point ", tex`p \in \R^2.` ])}** An example situation where this is useful, and the raison d'être for this notebook, is to compute local pairwise forces between point particles in a physics simulation. Do this efficiently by precomputing a two level lookup table where the points are grouped by their coordinates rounded to the nearest multiple of ${tex`r`}. This algorithm probably has an established name, but I didn't try to find out. Let me know about it! Another approach is to use a quadtree, as demoed in this other notebook by [@fil](https://observablehq.com/@fil): [Quadtree.findInCircle](https://observablehq.com/@d3/quadtree-findincircle?collection=@d3/d3-quadtree). My _intuition_ is that a quadtree probably gives faster queries, but that a lookup grid is faster to initialize. In my target use case, the physics simulation, I perform one search per particle per iteration. Particles generally move on each iteration, so that any construct of the particle positions used to facilitate the search must be recreated, or at least updated, on every iteration. The simulation is supposed to converge to an equilibrium, so there might be _some_ opportunity for reuse as the particles settle, but that also increases implementation complexity. Future work: - Benchmarking for different datasets. - Quantitative comparison to quadtree approach. - Investigate how the relation between point density and radius of influence affects performance. `32:Ta7d,{ let elapsedMillisMovingAverage = null; const movingAverageAlpha = 0.01; const lookupGrid = new LookupGrid(points, influenceRadius); const cursor = svg``; const candidates = svg``; const neighbors = svg``; const gridLinesV = d3 .range(0, width, influenceRadius * 2) .map(x => svg``); const gridLinesH = d3 .range(0, height, influenceRadius * 2) .map(y => svg``); const theSvg = html``; const log = html`

`; update({ x: width / 2, y: height / 2 }); function update(cursorPosition) { cursor.setAttribute("cx", cursorPosition.x); cursor.setAttribute("cy", cursorPosition.y); candidates.setAttribute( "d", computePointsPathAttrD( lookupGrid.getPointsInAdjacentCells(cursorPosition) ) ); const { result: neighborPoints, elapsedMillis } = time(() => lookupGrid.getNeighbors(cursorPosition) ); elapsedMillisMovingAverage = (1 - movingAverageAlpha) * (elapsedMillisMovingAverage ?? elapsedMillis) + elapsedMillis * movingAverageAlpha; log.innerHTML = ""; log.appendChild( md`Last search duration: ${elapsedMillis.toFixed( 4 )}ms.
Exponential moving average (${tex`\alpha`}=${movingAverageAlpha.toFixed( 3 )}): ${elapsedMillisMovingAverage.toFixed( 4 )}ms.
(Render duration ≠ Search duration)` ); neighbors.setAttribute("d", computePointsPathAttrD(neighborPoints)); } theSvg.addEventListener("pointermove", function(event) { const origin = theSvg.getBoundingClientRect(); const x = event.clientX - origin.x; const y = event.clientY - origin.y; update({ x, y }); }); return html`

${theSvg} ${log} The duration is measured for the lookup call exclusively. Keep in mind that rendering may be slow even though the search is fast. This demo also does a secondary search to find and highlight all points in the closest four grid cells, to illustrate the algorithm.

`; }33:T695,class LookupGrid { constructor(points, neighborhoodRadius) { this.r = neighborhoodRadius; this.rDouble = 2 * neighborhoodRadius; this.tree = d3.group( points, d => this.computeKey(d.x - this.r), d => this.computeKey(d.y - this.r), ) } getPointsInAdjacentCells(point) { const xa = this.computeKey(point.x); const ya = this.computeKey(point.y); const xb = xa - 1; const yb = ya - 1; // Get the points in all four grid cells around the // grid line crossing closest to the point. return [ this.getPointsInCell(xa, ya), this.getPointsInCell(xa, yb), this.getPointsInCell(xb, yb), this.getPointsInCell(xb, ya), ].flat() } getNeighbors(point, measureDistance = (dx, dy) => dx*dx + dy*dy) { const maxDistance = measureDistance(this.r, 0); return this.getPointsInAdjacentCells(point) // I considered using an intermediate filtering step here that removes // all points outside the bounding box with radius r around point, // but opted to just do the distance calculation right away. Two multiplications // and one addition for the squared distance shouldn't be too bad... .filter(({x, y}) => measureDistance(x - point.x, y - point.y) <= maxDistance) } // Using another distance measure is easy, but bear in mind the extents of the grid cells getManhattanNeighbors(point) { return this.getNeighbors(point, (dx, dy) => Math.abs(dx) + Math.abs(dy)); } getPointsInCell(xKey, yKey) { return this.tree.get(xKey)?.get(yKey) ?? []; } computeKey(coordinate) { return Math.round(coordinate / this.rDouble); } }34:T884,class LookupGridRn { // With points expressed as an array of their coordinates, we can extend to work in any dimensional space constructor(points, neighborhoodRadius) { this.r = neighborhoodRadius; this.rDouble = 2 * neighborhoodRadius; this.nDims = points[0]?.length ?? 0; const keyExtractors = Array.from( { length: this.nDims }, (_,iDim) => point => this.computeKey(point[iDim] - this.r) ); this.tree = d3.group( points, ...keyExtractors ) } getPointsInAdjacentCells(point) { const keys = this.computeKeys(point); const powers = keys.map((x, i) => Math.pow(2, i)).reverse(); // The keys for adjacent cells are produced by subtracting 1 from some coordinates. // The specific coordinates to modify to get cell n are those with ones // in row n of a boolean truth table with nDims variables. const adjacentKeySets = Array.from( { length: Math.pow(2, this.nDims) }, () => keys.slice() ).map((arr, i) => arr.map((x, j) => x - (Math.floor(i / powers[j]) % 2))) return adjacentKeySets.map(this.getPointsInCell, this).flat(); } getNeighbors(point, norm = point => Math.hypot(...point)) { const pointOnRadius = Array.from( {length: this.nDims}, (_, i) => i === 0 ? this.r : 0 ); const maxNorm = norm(pointOnRadius); return this.getPointsInAdjacentCells(point) .filter(candidate => norm(candidate.map((coordinate, i) => coordinate - point[i])) <= maxNorm) } // Using another distance measure is easy, but bear in mind the extents of the grid cells getManhattanNeighbors(point) { return this.getNeighbors(point, point => point.reduce((s, x) => s + Math.abs(x), 0)); } getPointsInCell(keys) { let subtree = this.tree; let i = 0; while (subtree && i < keys.length) { const key = keys[i++]; subtree = subtree.get(key); } return subtree instanceof Array ? subtree : []; } computeKeys(point) { return point.map(coordinate => this.computeKey(coordinate)); } computeKey(coordinate) { return Math.round(coordinate / this.rDouble); } }35:[] 2f:{"title":"Point Neighborhood Lookup Grid","cells":[{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":71,"value":"viewof pointCount = Range([0, 50000], {\n label: \"Number of points\",\n value: 2500,\n step: 1,\n transform: Math.sqrt\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":77,"value":"viewof influenceRadius = Range([0, 400], {\n label: \"Radius of influence, as percentage of maximum\",\n step: 1,\n value: Math.round(width / 15)\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":592,"value":"md`The orange dots are those found within the radius of influence of the mouse cursor, the so called _neighbors_. Black dots are _candidate neighbors_ that were returned from the lookup and then filtered by distance to get the actual neighbors.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":96,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":375,"value":"md`Expand the example instance below to explore its structure. The keys in the tree are multiples of the influence radius.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":135,"value":"lookupGridExampleInstance = new LookupGrid(points, influenceRadius)","pinned":false,"hidden":false,"mode":"ojs"},{"id":387,"value":"md`The next cell contains the implementation; a class called \\`${LookupGrid.name}\\`.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":175,"value":"$33","pinned":true,"hidden":false,"mode":"ojs"},{"id":514,"value":"md`\n---\n## Extension to ${tex`\\R^n`}\nJust for the challenge, I wrote \nanother version, \\`${\n LookupGridRn.name\n}\\`, that extends the concept to ${tex`\\R^n`}, with points represented as arrays of coordinates.\n\nOf course, the implementation is way more convoluted than the version specific for ${tex`\\R^2`}. First, lets look at some examples of instantiating \\`${\n LookupGridRn.name\n}\\` and performing some searches with it, then read on for the impementation.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":440,"value":"// Generate some high dimensional points\nhighDimPoints = createPoints(10, 4, () => Math.floor(Math.random() * 10))","pinned":true,"hidden":false,"mode":"ojs"},{"id":438,"value":"// Explore this example instance!\nhighDimLookupGrid = new LookupGridRn(highDimPoints, 10)","pinned":true,"hidden":false,"mode":"ojs"},{"id":469,"value":"// Using euclidean distance\nhighDimLookupGrid.getNeighbors(highDimPoints[0])","pinned":true,"hidden":false,"mode":"ojs"},{"id":530,"value":"// Using the included other example distance measure, Manhattan distance\nhighDimLookupGrid.getManhattanNeighbors(highDimPoints[0])","pinned":true,"hidden":false,"mode":"ojs"},{"id":535,"value":"// Another exotic distance measure\nnew LookupGridRn(highDimPoints, 4).getNeighbors(highDimPoints[0], point =>\n Math.max(...point)\n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":584,"value":"md`Now, here's the implementation:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":429,"value":"$34","pinned":true,"hidden":false,"mode":"ojs"},{"id":52,"value":"md`## Misc`","pinned":false,"hidden":false,"mode":"ojs"},{"id":567,"value":"points = Array.from({ length: pointCount }, (value, index) => {\n return { x: Math.random() * width, y: Math.random() * height };\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":432,"value":"function createPoints(\n nPoints,\n nDimensions = 2,\n coordinateInitializer = Math.random\n) {\n const pointOptions = { length: nDimensions };\n return Array.from({ length: nPoints }, () =>\n Array.from(pointOptions, coordinateInitializer)\n );\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":92,"value":"height = width / 2","pinned":false,"hidden":false,"mode":"ojs"},{"id":194,"value":"computePointsPathAttrD = points =>\n points.map(({ x, y }) => `M${x},${y}h0`).join(\" \")","pinned":false,"hidden":false,"mode":"ojs"},{"id":124,"value":"html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":329,"value":"function time(f, ...args) {\n const t1 = performance.now();\n const result = f(...args);\n const t2 = performance.now();\n return { result, elapsedMillis: t2 - t1 };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":66,"value":"import { Range } from \"@observablehq/inputs\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":57,"value":"d3 = require(\"d3@6\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":28,"value":"compileText = pieces => md`${pieces.map(piece => html`${piece}`)}`","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W35","stdlib":"1","titling":"auto"}