计算周期性二维晶格中的成对空间距离

Question

我一直在徒劳地寻找一种解决方案来计算 periodic/wrapped 边格子上离散的等距正方形站点之间的空间距离。例如在9x9的格子上设置相邻站点为2个单位：

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')

其中站点用索引注释并且红色箭头表示与站点 0（2 个单位）的等距离，我如何创建一个数组 sdist 给我站点之间的最短成对距离，即 sdist[7,1]=2，sdist[2,6] = 2.8284。常用工具 scipy.spatial.distance.cdist 似乎不适用于此处。

Answer 1

由于 distance.cdist 接受任意指标，作为可调用提供，问题在于为您的指标编写函数。

如果它是点 p 和 q 之间的包裹距离，那就是

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

其中m是晶格的大小（遗憾的是cdist不支持向距离函数传递额外的参数，所以m必须取更大的范围。 )

但在您的情况下，您希望边长为 1 的正方形之间的距离最小。这意味着在计算包裹差异向量 np.minimum(diff, m - diff) 之后，我们可以将其每个分量减少 1，因为，比如说，两个正方形点之间的最小 x 差比这些正方形中心之间的 x 差小 1。当然这个减法不应该造成差异负。所以函数变成了

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))

其中 clip 处理两个中心具有相同 x 坐标或相同 y 坐标的情况。

剩下的就是两行（不算from scipy.spatial import distance）：

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

示例的输出：

[[ 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.        ]]

Answer 2

我试图理解请求，我的解释在这里说明

我没有为格点编号，而是从 1 开始编号

然后我 'augmented' lattice 中给出的单元格通过切片并使用 np.concatenate 在起始格子的每个边界上添加半个单元格图中的黄线显示单元格寄宿生

给定起始格点标签编号，我找到我的 'fullLaug' 数组中的所有实例并计算平方欧氏距离，并按它排序

根据示例，我取最短的一条线，从起点 1 到 'fullLaug' 中每个其余点的最短距离实例，并打印实际距离

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.    ]]

计算周期性二维晶格中的成对空间距离

Calculating pairwise spatial distances in periodic 2D lattice

python

arrays

numpy

scipy

euclidean-distance