在考慮週期邊界條件的同時優化Python距離計算

我已經編寫了一個Python腳本來計算三維空間中兩個點之間的距離，同時考慮了週期邊界條件。問題是我需要對許多點進行計算，計算速度很慢。這是我的功能。在考慮週期邊界條件的同時優化Python距離計算

def PBCdist(coord1,coord2,UC): 
    dx = coord1[0] - coord2[0] 
    if (abs(dx) > UC[0]*0.5): 
     dx = UC[0] - dx 
    dy = coord1[1] - coord2[1] 
    if (abs(dy) > UC[1]*0.5): 
     dy = UC[1] - dy 
    dz = coord1[2] - coord2[2] 
    if (abs(dz) > UC[2]*0.5): 
     dz = UC[2] - dz 
    dist = np.sqrt(dx**2 + dy**2 + dz**2) 
    return dist

那麼我所說的功能，從而

for i, coord2 in enumerate(coordlist): 
    if (PBCdist(coord1,coord2,UC) < radius): 
     do something with i

最近我讀，我可以極大地使用列表理解提高性能。對於非PBC情況下工作，而不是爲PBC情況

coord_indices = [i for i, y in enumerate([np.sqrt(np.sum((coord2-coord1)**2)) for coord2 in coordlist]) if y < radius] 
for i in coord_indices: 
    do something

是否有某種方式來做到這對中國人民銀行情況相同呢？有沒有更好的替代方案？

來源

2012-06-19 johnjax

您正在使用NumPy，因此您應該引導循環以提高性能。「coordlist」的結構究竟是什麼？它應該是一個二維NumPy數組，以便能夠使用NumPy ufuncs優化循環。 –

coordlist是一個形狀約爲（5711,3）的numpy數組。 coordlist本身來自一個更大的列表，所以我有效地循環了coordlist 20,000次，並且coordlist的列表循環了大約50次......你得到了圖片。 – johnjax

我查閱了NumPy中的矢量化函數。該文檔說：[「矢量化函數主要是爲了方便，而不是爲了性能。該實現本質上是一個for循環。」]（http://docs.scipy.org/doc/numpy/reference/generated/numpy。 vectorize.html） – johnjax

你應該寫你的distance()功能的方式可以通過5711點矢量化循環。以下實現接受點作爲任一x0或x1參數的數組：

def distance(x0, x1, dimensions): 
    delta = numpy.abs(x0 - x1) 
    delta = numpy.where(delta > 0.5 * dimensions, delta - dimensions, delta) 
    return numpy.sqrt((delta ** 2).sum(axis=-1))

實施例：

>>> dimensions = numpy.array([3.0, 4.0, 5.0]) 
>>> points = numpy.array([[2.7, 1.5, 4.3], [1.2, 0.3, 4.2]]) 
>>> distance(points, [1.5, 2.0, 2.5], dimensions) 
array([ 2.22036033, 2.42280829])

的結果作爲第二個參數傳遞點之間的距離的陣列distance()和各指向points。

來源

2012-06-19 20:58:26

這導致我的代碼加快了5倍。謝謝！對於那些尋找替代品的人來說，lazyr的答案也表現得很好。 – johnjax

採取規範後沒有什麼區別，但我喜歡通過用delta = numpy.where（「，delta - dimensions，」）替換第三行來得到delta上的正確符號。還要注意''np.hypot'在避免溢出方面比sqrt更好（sum（x ** 2）） – Patrick

@Patrick'numpy.hypot（）'只適用於兩維，而OP需要適用於三維三維點。關於'delta'的標誌，好吧，如果你在意，可以隨意編輯帖子。 :) –

看看Ian Ozsvalds high performance python tutorial。它包含很多關於下一步你可以去哪裏的建議。

包括：

矢量
用Cython
pypy
numexpr
pycuda
multiprocesing
平行蟒

來源

2012-06-19 20:33:20

import numpy as np 

bounds = np.array([10, 10, 10]) 
a = np.array([[0, 3, 9], [1, 1, 1]]) 
b = np.array([[2, 9, 1], [5, 6, 7]]) 

min_dists = np.min(np.dstack(((a - b) % bounds, (b - a) % bounds)), axis = 2) 
dists = np.sqrt(np.sum(min_dists ** 2, axis = 1))

這裏a和b是向量的名單要計算距離之間bounds是空間的界限（所以這裏所有三個維度從0到10，然後換行）。它計算了a[0]和b[0],a[1]和b[1]之間的距離，依此類推。

我相信numpy的專家可以做的更好，但是這可能會是一個數量級比你在做什麼更快，因爲大部分的工作在C

現在做

來源

2012-06-19 20:52:15

我也嘗試了這種方法。它還導致代碼增加了大約5倍。不幸的是，我只能選擇1個答案作爲正確答案：/ – johnjax

@johnjax如果值得的話，我會接受Sven Marnach的回答，如果我一直在你的鞋子裏。它比我的更直接適用。 –

我發現meshgrid對於生成距離非常有用。例如：

import numpy as np 
row_diff, col_diff = np.meshgrid(range(7), range(8)) 
radius_squared = (row_diff - x_coord)**2 + (col_diff - y_coord)**2

我現在有一個數組（radius_squared）其中每一條目指定從陣列位置[x_coord, y_coord]的距離的平方。

以環陣，我可以做到以下幾點：

row_diff, col_diff = np.meshgrid(range(7), range(8)) 
row_diff = np.abs(row_diff - x_coord) 
row_circ_idx = np.where(row_diff > row_diff.shape[1]/2) 
row_diff[row_circ_idx] = (row_diff[row_circ_idx] - 
         2 * (row_circ_idx + x_coord) + 
         row_diff.shape[1]) 
row_diff = np.abs(row_diff) 
col_diff = np.abs(col_diff - y_coord) 
col_circ_idx = np.where(col_diff > col_diff.shape[0]/2) 
col_diff[row_circ_idx] = (row_diff[col_circ_idx] - 
         2 * (col_circ_idx + y_coord) + 
         col_diff.shape[0]) 
col_diff = np.abs(row_diff) 
circular_radius_squared = (row_diff - x_coord)**2 + (col_diff - y_coord)**2

我現在都用矢量數學覽的排列距離。

來源

2014-11-21 17:35:34 carsonc

在考慮週期邊界條件的同時優化Python距離計算

回答

相關問題