Расстояние между точками в массиве с циклическими границами - PullRequest
Новая
       19

Расстояние между точками в массиве с циклическими границами

0 голосов
/ 04 октября 2018

Я знаю, как рассчитать евклидово расстояние между точками в массиве, используя scipy.spatial.distance.cdist 

import scipy.spatial.distance as scidist
import numpy as np

n=4
a=np.zeros([n,n])
i=np.argwhere(a>-1)
dist1=scidist.cdist(i,i,metric='euclidean')

Аналогично ответам на этот вопрос: Рассчитать расстояния между одной точкой в ​​матрицеИз всех других точек

Однако я хотел бы сделать расчет, предполагая циклические граничные условия , например, чтобы точка [0,0] была расстоянием 1 от точки [0,n-1] в данном случае, а не расстояние n-1.

Единственный способ, о котором я могу думать, - это повторить вычисление 4 раза, при этом индексы доменов имеют n, вычтенное из x, y изатем x & y направления, а затем суммируем результаты и находим минимум по 4 срезами: 

distl=[]
for xoff in [0,n]:
    for yoff in [0,n]:
        j=i.copy()
        j[:,0]-=xoff
        j[:,1]-=yoff
        distl.append(scidist.cdist(i,j,metric='euclidean'))
dist=np.amin(np.dstack(distl),2)

Я думаю (!) это работает, но это медленно для моих больших массивов и немного уродливо -там лучше / быстрее?

1 Ответ

0 голосов
/ 06 октября 2018

Вот фиксированная версия вашего кода и другой метод, который немного быстрее.Они дают одинаковые результаты, поэтому я уверен, что они верны: 

import numpy as np
from scipy.spatial.distance import squareform, pdist, cdist
from numpy.linalg import norm

def pb_OP(A, p=1.0):
    distl = []
    for *offs, ct in [(0, 0, 0), (0, p, 1), (p, 0, 1), (p, p, 1), (-p, p, 1)]:
        B = A - offs
        distl.append(cdist(B, A, metric='euclidean'))
        if ct:
            distl.append(distl[-1].T)
    return np.amin(np.dstack(distl), 2)

def pb_pp(A, p=1.0):
    out = np.empty((2, A.shape[0]*(A.shape[0]-1)//2))
    for o, i in zip(out, A.T):
        pdist(i[:, None], 'cityblock', out=o)
    out[out > p/2] -= p
    return squareform(norm(out, axis=0))

test = np.random.random((1000, 2))

assert np.allclose(pb_OP(test), pb_pp(test))

from timeit import timeit
t_OP = timeit(lambda: pb_OP(test), number=10)*100
t_pp = timeit(lambda: pb_pp(test), number=10)*100
print('OP', t_OP)
print('pp', t_pp)

Пробный прогон.1000 баллов: 

OP 210.11001259903423
pp 22.288734700123314

Мы видим, что мой метод в ~ 9 раз быстрее, что, по точному совпадению, является количеством смещенных cponfigurations, которые должна проверить версия OP.Он использует pdist для отдельных координат, чтобы получить абсолютные различия.Там, где они больше половины шага сетки, мы вычитаем один период.Осталось взять евклидову норму и распаковать хранилище.

...