Calculating pairwise spatial distances in periodic 2D lattice

Question

I've been searching in vain for a solution to compute spatial distances between discrete, equidistant square sites on a lattice of periodic/wrapped edges. For example, setting adjacent sites to be 2 units apart on 9x9 lattice:

m = 9
lattice = np.zeros((m,m))
for i in arange(0,6+1,3):
    for j in arange(0,6+1,3):
        lattice[i,j]=1

In  []: lattice
Out []: array([[ 1.,  0.,  0.,  1.,  0.,  0.,  1.,  0.,  0.],
               [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
               [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
               [ 1.,  0.,  0.,  1.,  0.,  0.,  1.,  0.,  0.],
               [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
               [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
               [ 1.,  0.,  0.,  1.,  0.,  0.,  1.,  0.,  0.],
               [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
               [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])

In  []: plt.imshow(lattice, origin='origin', cmap='binary', interpolation='none')

where the sites are annotated with indices and the red arrows indicate equidistances from site 0 (2 units), how can I create an array sdist that gives me the shortest pairwise distances between sites, i.e. sdist[7,1]=2, sdist[2,6] = 2.8284. The usual tool scipy.spatial.distance.cdist doesn't seem to apply here.

Answer 1

I tried to understand the request, my interpertaion is illustrated here

Instead of ones for the lattice points, I numbered them, starting from 1

I then 'augmented' the cell given in lattice by slicing and using np.concatenate to add half cells on each boarder of starting lattice Yellow lines on the plot show cell boarders

Given starting lattice point label number, I find all instances in my 'fullLaug' array and calc squared Euclidean disatnce, and sort by it

Per the example I take the shortest and plot a line from the starting point 1 to the shortest distance instance of the each of the rest of the points in 'fullLaug' as well as printing the actual distance

import numpy as np
import matplotlib.pyplot as plt

" create lattice, lattice points nonzero, incremtning count "

m = 9
pcnt = 0
lattice = np.zeros((m, m))
for i in range(0, 6+1, 3):
    for j in range(0, 6+1, 3):
        pcnt += 1
        lattice[i, j] = pcnt # lable lattice points with count

#print(*lattice, sep='\n')

"augment lattice with duplicated halves, up/down, left/right"

def halves(a):
    n = len(a)
    return a[:(n+n%2)//2-n%2], a[(n-n%2)//2+n%2:]


rightL, leftL = halves(lattice.T)
n = len(rightL)
cornerL = np.zeros((n, n))
rightL, leftL = rightL.T, leftL.T

rightLaug = np.concatenate((cornerL, rightL, cornerL), axis=0)
leftLaug = np.concatenate((cornerL, leftL, cornerL), axis=0)

lowL, upL = halves(lattice)

centerL = np.concatenate((upL, lattice, lowL), axis=0)

fullLaug = np.concatenate((leftLaug, centerL, rightLaug), axis=1)

"plot fully agumented lattice, yellow lines are boarders of cells"

plt.imshow(fullLaug, origin='origin', cmap='jet', interpolation='none')

plt.plot((n, n), (0, 2*n+m-1), 'y', (n+m-1, n+m-1), (2*n+m-1, 0), 'y')
plt.plot((0, 2*n+m-1), (n, n), 'y', (2*n+m-1, 0), (n+m-1, n+m-1), 'y')

#print(*zip(*np.where(fullLaug == 3)))


def distsq(a, b, n):
    """takes lattice point pcnt labels a, b
    finds a indices in lattice, offests to indices in fullaug
    finds all instances of b indices in fullaug
    calcs all a-b distances
    returns sorted list of distances, indices, a-a appended to start
    """
    a = list(*zip(*np.where(lattice == a)))
    a[0] += n; a[1] += n

    bz = zip(*np.where(fullLaug == b))

    ds = [[(a[0] - b[0])**2 + (a[1] - b[1])**2, [*b]] for b in bz]

    return [[0, a]] + sorted(ds, key=lambda x: x[0])

""" plot shortest distance lines from point a to each point in lattice
    fill dist_lattice"""

dist_lattice = np.zeros((m, m))
for i in range(1, pcnt + 1):
    dps = distsq(1, i, n)[0:2]
    dist_lattice[np.where(lattice == i)] = np.sqrt(dps[1][0])
    plt.plot(*zip(dps[0][1], dps[1][1]))

np.set_printoptions(precision=4, threshold=1000)
print(dist_lattice)

[[ 0.      0.      0.      3.      0.      0.      3.      0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.      0.      0.      0.    ]
 [ 3.      0.      0.      4.2426  0.      0.      4.2426  0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.      0.      0.      0.    ]
 [ 3.      0.      0.      4.2426  0.      0.      6.7082  0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.      0.      0.      0.    ]]

Answer 2

Since distance.cdist accepts an arbitrary metric, provided as a callable, the problem is just in writing a function for your metric.

If it was wrapped distance between points p and q, that would be

def wrapped_euclidean_points(p, q):
    diff = np.abs(p - q)
    return np.linalg.norm(np.minimum(diff, m - diff))

where m is the size of the lattice (unfortunately cdist does not support passing extra parameters to the distance function, so m has to be taken from the larger scope.)

But in your case, you want minimal distances between squares of sidelength 1. What this means is that after computing the wrapped-difference vector np.minimum(diff, m - diff) we can reduce each of its components by 1, because, say, the smallest x-difference between points of two squares is 1 less than the x-difference between the centers of those squares. Of course this subtraction should not make the difference negative. So the function becomes

def wrapped_euclidean_squares(p, q):
    diff = np.abs(p - q)
    return np.linalg.norm(np.clip(np.minimum(diff, m - diff) - 1, 0, None))

where clip takes care of the case when two centers had the same x-coordinate or the same y-coordinate.

The rest is just two lines (not counting from scipy.spatial import distance):

coords = np.vstack(np.nonzero(lattice)).T
dist = distance.cdist(coords, coords, metric=wrapped_euclidean_squares)

Output for your example:

[[ 0.          2.          2.          2.          2.82842712  2.82842712       2.          2.82842712  2.82842712]
 [ 2.          0.          2.          2.82842712  2.          2.82842712       2.82842712  2.          2.82842712]
 [ 2.          2.          0.          2.82842712  2.82842712  2.       2.82842712  2.82842712  2.        ]
 [ 2.          2.82842712  2.82842712  0.          2.          2.          2.       2.82842712  2.82842712]
 [ 2.82842712  2.          2.82842712  2.          0.          2.       2.82842712  2.          2.82842712]
 [ 2.82842712  2.82842712  2.          2.          2.          0.       2.82842712  2.82842712  2.        ]
 [ 2.          2.82842712  2.82842712  2.          2.82842712  2.82842712       0.          2.          2.        ]
 [ 2.82842712  2.          2.82842712  2.82842712  2.          2.82842712       2.          0.          2.        ]
 [ 2.82842712  2.82842712  2.          2.82842712  2.82842712  2.          2.       2.          0.        ]]