我正試圖模擬一些二維粒子。每個粒子都是一個有一個方向的圓。該方向由一個二維單位矢量指定。
在我的模擬的一個部分中，我想計算一個粒子對之間的角度和每個粒子的方向的函式。這應該是針對每個粒子對進行的。在視覺上，我想計算一個$ heta_i$和$ heta_j$的函式（見圖片鏈接）。
我已經為每個粒子對計算了成對的位移單位向量。這是一個名為r的numpy陣列，其形狀為(N, N, 2)，其中N是粒子的總數。我還計算了每個粒子在直角坐標中的方向。這是一個名為orientation of shape (N, 2)的numpy陣列，
。我已經能夠把我需要的代碼寫成一個雙for-loop。
output = np.zeros(r.shape)
for i in range（len（r））。
for j in range(len(r[i])) 。
if（i < j）。
theta_i = angle_between(orientation[i], r[i,j])
theta_j = angle_between(orientation[j], r[i,j])
output[i][j] = (1 / (1 np.exp(-w*(np.cos(theta_i) - np.cos(alpha)))))) * (1 / (1 np.exp(-w*(np.cos(theta_j) - np.cos(alpha))))))
不過，運行這段代碼需要的時間太長了，尤其是對于大型系統。是否有辦法將雙倍for回圈矢量化，使其運行更快？
uj5u.com熱心網友回復：
下面是我運行你的代碼的一個例子：
2)
moment_cartesian = np.random.rand(N, 2)
output = np.zeros(r_hat.shape)
輸出 = np.zeros((N, N))
w = 1
alpha = 1for i in range（len（r_hat））。
for j in range(len(r_hat[i]) ):
if（i < j）。
theta_i = angle_between(moment_cartesian[i], r_hat[i,j])
theta_j = angle_between(moment_cartesian[j], r_hat[i,j])
theta_i_ij[i, j] = theta_i
theta_j_ij[i, j] = theta_j
output[i][j] = (1 / (1 np.exp(-w*(np.cos(theta_i) - np.cos(alpha)))))) * (1 / (1 np.exp(-w*(np.cos(theta_j) - np.cos(alpha))))))
print（f'theta_i_ij =
{theta_i_ij}'/span>)
print(f'theta_j_ij =
{theta_j_ij}')
print（f'output =
{output}')
import numpy as np
np.random.seed(0)
N = 3
輸出：
theta_i_ij =
[[0. 0.99438035 0.98889055]
[0. 0. 0.9954101 ]
[0. 0. 0. ]]
theta_j_ij =
[[0. 0.99873953 0.99891568]
[0. 0. 0.90135909]
[0. 0. 0. ] ]
輸出=
[[0. 0.25072287 0.25127888]
[0. 0. 0.26052633]
[0. 0. 0. ]]
下面是一個代碼矢量化的例子：
innerp = np.einsum('jk,ijk->ij', moment_cartesian, r_hat)
norm = np.linalg.norm(moment_cartesian, axis = -1) [None, :]*np.linalg.norm(r_hat, axis = -1)
theta_j_v = innerp/norm
innerp = np.einsum('ik,ijk->ij', moment_cartesian, r_hat)
norm = np.linalg.norm(moment_cartesian, axis = -1)[:, None]*np.linalg.norm(r_hat, axis = -1)
theta_i_v = innerp/norm
output_v = (1 / (1 np.exp(-w*(np.cos(theta_i_v) - np.cos(alpha)))))) * (1 / (1 np.exp(-w*(np.cos(theta_j_v) - np.cos(alpha))))))
print（f'theta_i_v =
{theta_i_v}')
print（f'theta_j_v =
{theta_j_v}')
print（f'output_v =
{output_v}')
輸出：
theta_i_v =
[[0.99717384 0.99438035 0.98889055]
[0.9090371 0.95351666 0.9954101 ]
[0.99986414 0.98876437 0.87290329]]
theta_j_v =
[[0.99717384 0.99873953 0.99891568]
[0.96281813 0.95351666 0.90135909]
[0.98397106 0.9796661 0.87290329] ]
output_v =
[[0.25059435 0.25072287 0.25127888]
[0.26327747 0.25972068 0.26052633]
[0.2516916 0.25331216 0.27620631 ]]
但我必須說，我沒有想到有什么方法可以將你的代碼矢量化而不需要不必要地計算矩陣的下三角。 實際上，你根本不需要對你的代碼進行矢量化。 如果你使用 Numba 來編譯一個包含你的代碼的函式，它將執行同樣的操作，如果不是比矢量化的解決方案更快的話。
下面是我的 Numba 例子：
from numba import jit
)
輸出：
output =
[[0. 0.25072287 0.25127888]
[0. 0. 0.26052633]
[0. 0. 0. ]]
uj5u.com熱心網友回復：
如果你正在尋找在JAX中這樣做的方法，這個notebook提供了一個很好的視覺教程。
轉載請註明出處，本文鏈接：https://www.uj5u.com/shujuku/326860.html
標籤：
下一篇：尋找影像中一條線上的像素坐標