The neighbours within a radius¶

Every pair of a list's items no farther apart than r, found through the bins.

Requires

each item's Position p (m, a vector of the space's dimension), the Radius r (m, positive).

Provides

the neighbour relation, the pairs (i, j) with i < j and |p_i - p_j| <= r, each once, ordered by i then j.

Law

the space is binned in cubes of side r, an item in the bin floor(p / r); any two items within r of each other lie in the same bin or in adjacent ones, so an item's candidates are the items of its bin and the 3^d - 1 around it, and a candidate is a neighbour where the squared distance is at most r^2.

Source

Hockney & Eastwood, Computer Simulation Using Particles, 1988, 8-4 (the chaining mesh); Allen & Tildesley, Computer Simulation of Liquids, 2nd ed., 5.3 (cell lists).

Proof

invariant, limit.

No Trait form: waits on the middle: a key kernel building the relation of a list's near pairs that a law folds through, which draft 5 does not write yet; until then a reference.

space.neighbours_within(position, radius)

Every pair of a list's items no farther apart than r, found through the bins: the relation a law over near pairs folds through, contact's broad phase among them.

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import library
from lusine.kernels.space.neighbours_within import neighbours_within

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) import numpy as np import matplotlib.pyplot as plt from lusine import library from lusine.kernels.space.neighbours_within import neighbours_within

The reference¶

The kernel's reference, in numpy; its Trait form is not written yet.

In [2]:

library.reference(neighbours_within)

library.reference(neighbours_within)

Out[2]:

the reference, neighbours_within.py

def neighbours_within(position, radius):
    """THE NEIGHBOURS WITHIN A RADIUS: every pair of a list's items no
    farther apart than r, found through the bins.

    Requires: each item's `Position` p (m, a vector of the space's
    dimension), the `Radius` r (m, positive).
    Provides: the neighbour relation, the pairs (i, j) with i < j and
    |p_i - p_j| <= r, each once, ordered by i then j.
    Law: the space is binned in cubes of side r, an item in the bin
    floor(p / r); any two items within r of each other lie in the same bin
    or in adjacent ones, so an item's candidates are the items of its bin and
    the 3^d - 1 around it, and a candidate is a neighbour where the squared
    distance is at most r^2.
    Source: Hockney & Eastwood, Computer Simulation Using Particles, 1988,
    8-4 (the chaining mesh); Allen & Tildesley, Computer Simulation of
    Liquids, 2nd ed., 5.3 (cell lists).
    Proof: invariant, limit."""
    p = np.asarray(position, dtype=np.float64)
    if p.ndim != 2:
        raise ValueError(f"positions are a list of vectors, shape (items, dimension), not {p.shape}")
    r = float(radius)
    if not r > 0.0:
        raise ValueError(f"the radius is positive, not {radius!r}")
    keys = np.floor(p / r).astype(np.int64)
    bins = {}
    for i, key in enumerate(map(tuple, keys.tolist())):
        bins.setdefault(key, []).append(i)
    around = list(itertools.product((-1, 0, 1), repeat=p.shape[1]))
    rr = r * r
    pairs = []
    for i, key in enumerate(map(tuple, keys.tolist())):
        for d in around:
            for j in bins.get(tuple(a + b for a, b in zip(key, d)), ()):
                if j > i:
                    v = p[j] - p[i]
                    if float(np.dot(v, v)) <= rr:
                        pairs.append((i, j))
    pairs.sort()
    out = np.asarray(pairs, dtype=np.int64).reshape(-1, 2)
    return out[:, 0], out[:, 1]

Invariant: the bins find what every pair's distance finds¶

The bins are only a search: the pairs they find are exactly those a test of every pair finds, each once.

In [3]:

rng = np.random.default_rng(1)
p = rng.uniform(0, 10, (400, 2))
i, j = neighbours_within(p, 0.6)
d2 = ((p[:, None, :] - p[None, :, :]) ** 2).sum(-1)
want = [(a, b) for a in range(len(p)) for b in range(a + 1, len(p)) if d2[a, b] <= 0.36]
assert list(zip(i.tolist(), j.tolist())) == want
fig, ax = plt.subplots(figsize=(4.5, 4.5))
ax.plot(*p.T, '.', ms=3)
for a, b in zip(i, j):
    ax.plot(*p[[a, b]].T, lw=0.6, color='tab:orange')
ax.set_aspect('equal'); ax.set_title(f'{len(i)} pairs within 0.6')
plt.show()

rng = np.random.default_rng(1) p = rng.uniform(0, 10, (400, 2)) i, j = neighbours_within(p, 0.6) d2 = ((p[:, None, :] - p[None, :, :]) ** 2).sum(-1) want = [(a, b) for a in range(len(p)) for b in range(a + 1, len(p)) if d2[a, b] <= 0.36] assert list(zip(i.tolist(), j.tolist())) == want fig, ax = plt.subplots(figsize=(4.5, 4.5)) ax.plot(*p.T, '.', ms=3) for a, b in zip(i, j): ax.plot(*p[[a, b]].T, lw=0.6, color='tab:orange') ax.set_aspect('equal'); ax.set_title(f'{len(i)} pairs within 0.6') plt.show()

Invariant: the storage order changes nothing¶

The items stored in another order are the same pairs, named by the same items.

In [4]:

perm = rng.permutation(len(p))
pi, pj = neighbours_within(p[perm], 0.6)
assert sorted(tuple(sorted((int(perm[a]), int(perm[b])))) for a, b in zip(pi, pj)) == want
print('the same relation')

perm = rng.permutation(len(p)) pi, pj = neighbours_within(p[perm], 0.6) assert sorted(tuple(sorted((int(perm[a]), int(perm[b])))) for a, b in zip(pi, pj)) == want print('the same relation')

the same relation

Limit: exactly r, and everything¶

A pair exactly r apart is within it; a radius past the whole set pairs every item with every other.

In [5]:

i, j = neighbours_within(np.array([[0.0, 0], [1, 0], [2, 0], [3.5, 0]]), 1.0)
assert list(zip(i.tolist(), j.tolist())) == [(0, 1), (1, 2)]
i, _ = neighbours_within(rng.uniform(0, 1, (10, 3)), 5.0)
assert len(i) == 45
print('the edge case and the whole set')

i, j = neighbours_within(np.array([[0.0, 0], [1, 0], [2, 0], [3.5, 0]]), 1.0) assert list(zip(i.tolist(), j.tolist())) == [(0, 1), (1, 2)] i, _ = neighbours_within(rng.uniform(0, 1, (10, 3)), 5.0) assert len(i) == 45 print('the edge case and the whole set')

the edge case and the whole set