1

現在、Python で大きな for ループを使用して記述したコードをベクトル化しようとしています。ベクトル化されたコードは次のとおりです。

rho[pi,pj] += (rho_coeff*dt)*i_frac*j_frac
rho[pi+1,pj] += (rho_coeff*dt)*ip1_frac*j_frac
rho[pi,pj+1] += (rho_coeff*dt)*i_frac*jp1_frac
rho[pi+1,pj+1] += (rho_coeff*dt)*ip1_frac*jp1_frac

、 、、 、 、のそれぞれは、pi1次元ですべて同じ長さの numpy 配列です。は、2 次元の numpy 配列です。行列のどの要素が変更されたかを示す座標 ( 、 ) のリストを作成します。この変更には、( , ) 要素への用語の追加と、隣接する要素への同様の用語の追加 (( +1, )、( , +1) および ( +1, +1) が含まれます。リスト ( 、) 内の各座標には、固有の 、、、およびが関連付けられています。pjdti_fracj_fracip1_fracjp1_fracrhopipjpipjrho(rho_coeff*dt)*i_frac*j_fracpipjpipjpipjpipjpipjdti_fracj_fracip1_fracjp1_frac

問題は、リストが繰り返し座標を持つことができる (そして常に持つ) ことです。したがってrho、リスト内で同じ座標が検出されるたびに連続して追加する代わりに、最後の繰り返し座標に対応する用語のみを追加します。この問題は、Tentative Numpy Tutorialの仮の Numpy チュートリアルで、インデックスの配列を使用した凝ったインデックス付けの例で簡単に説明されています (ブール値のインデックス付けの前の最後の 3 つの例を参照してください)。残念ながら、彼らはこれに対する解決策を提供しませんでした。

forループに頼らずにこの操作を行う方法はありますか? パフォーマンスを最適化しようとしており、可能であればループをなくしたいと考えています。

参考までに: このコードは、各粒子からの電荷が、体積分率に基づいて粒子の位置を囲むメッシュの 4 つの隣接するノードに追加される 2D 粒子追跡アルゴリズムの一部を形成します。

4

1 に答える 1

1

配列を更新する前に、繰り返されるアイテムを見つけてそれらを追加する必要があります。次のコードは、最初の更新でそれを行う方法を示しています。

rows, cols = 100, 100
items = 1000

rho = np.zeros((rows, cols))
rho_coeff, dt, i_frac, j_frac = np.random.rand(4, items)
pi = np.random.randint(1, rows-1, size=(items,))
pj = np.random.randint(1, cols-1, size=(items,))

# The following code assumes pi and pj have the same dtype
pij = np.column_stack((pi, pj)).view((np.void,
                                      2*pi.dtype.itemsize)).ravel()

unique_coords, indices = np.unique(pij, return_inverse=True)
unique_coords = unique_coords.view(pi.dtype).reshape(-1, 2)
data = rho_coeff*dt*i_frac*j_frac
binned_data = np.bincount(indices, weights=data)
rho[tuple(unique_coords.T)] += binned_data

上記の一意の座標検出をすべて他の更新に再利用できると思うので、次のようにします。

ip1_frac, jp1_frac = np.random.rand(2, items)

unique_coords[:, 0] += 1
data =  rho_coeff*dt*ip1_frac*j_frac
binned_data = np.bincount(indices, weights=data)
rho[tuple(unique_coords.T)] += binned_data

unique_coords[:, 1] += 1
data =  rho_coeff*dt*ip1_frac*jp1_frac
binned_data = np.bincount(indices, weights=data)
rho[tuple(unique_coords.T)] += binned_data

unique_coords[:, 0] -= 1
data =  rho_coeff*dt*i_frac*jp1_frac
binned_data = np.bincount(indices, weights=data)
rho[tuple(unique_coords.T)] += binned_data
于 2013-07-22T17:44:48.273 に答える