1

次のことができる方法を知りたいです。

def GetFlux(self, time):
    bx = self.GetField("bx", time) * self.wpewce
    by = self.GetField("by", time) * self.wpewce
    bz = self.GetField("bz", time) * self.wpewce              

    flux  = np.zeros((self.ncells[0]+1,self.ncells[1]+1),"float32", order='FORTRAN')
    flux2  = np.zeros((self.ncells[0]+1,self.ncells[1]+1),"float32", order='FORTRAN')

    dx = self.dl[0]
    dz = self.dl[1]

    nx = self.ncells[0]
    nz = self.ncells[1]

    j = 0

    for i in np.arange(1, nx):
        flux2[i,0] = flux2[i-1,0] + bz[i-1,0]*dx
    flux[1:,0] = flux[0,0] + np.cumsum(bz[:-1,0]*dx)

    for j in np.arange(1,nz):
        flux2[0,j] = flux2[0,j-1] - bx[0,j-1]*dz
    flux[0,1:] = flux[0,0] - np.cumsum(bx[0,:-1]*dz)

    for i in np.arange(1,nx):
        for j in np.arange(1,nz):
            flux2[i,j] = 0.5*(flux2[i-1,j] + bz[i-1,j]*dx) + 0.5*(flux2[i,j-1] - bx[i,j-1]*dz)

    return flux2

ただし、非常に長い時間がかかる 2 つのネストされたループがありません。Bx、Bzおよびflux同じサイズの配列です。

最初の 2 つの単一ループを配列インデックスと cumsum で置き換えることができましたが、ネストされたループを置き換える方法がわかりません。

何か案が ？

ありがとう

4

3 に答える 3

1

内側のループは、ベクトル化するのがかなり簡単です。次のような基本的な方程式があります。

X[n] = a * X[n-1] + b[n]

この方程式は、X[n-1]に依存することなく拡張および書き換えることができます。

X[n] = a^n * X[0] + a^(n-1) * b[0] + a^(n-2) * b[1] + ... + a^0 * b[n]

したがって、元のコードが次のようになっている場合：

for i in np.arange(1,nx+1):
    for j in np.arange(1,nz+1):
        flux2[i,j] = 0.5*(flux2[i-1,j] + bz[i-1,j]*dx) \
                   + 0.5*(flux2[i,j-1] - bx[i,j-1]*dz)

あなたはこのように内側のループを取り除くことができます：

a = 0.5
aexp = np.arange(nz).reshape(nz, 1) - np.arange(nz).reshape(1, nz)
abcoeff = a**aexp
abcoeff[aexp<0] = 0
for i in np.arange(1,nx+1):
    b = 0.5*flux2[i-1, 1:] + 0.5*bz[i-1, 1:]*dx - 0.5*bx[i,:-1]*dz
    bvals = (abcoeff * b.reshape(1, nz)).sum(axis=1)
    n = np.arange(1, nz+1)
    x0 = flux2[i, 0]
    flux2[i, 1:] = a**n * x0 + bvals

浮動小数点エラーのため、値はまったく同じにはなりませんが、十分に近い値になります。理論的には、同じ手順を適用して両方のループを取り除くことができると思いますが、それは非常に複雑になり、アレイの形状によっては、パフォーマンス上のメリットがあまりない場合があります。

于 2012-08-08T18:17:41.980 に答える
1

この種の問題に対して scipy.ndimage.convolve を (ab) 使用する可能性があります。scipy でいくつかのフィルター メソッドを使用することも機能する可能性があり、適切に機能する scipy.ndimage.convolve に依存しないため、より良いものになる可能性があります (そして、これは遠い将来に変更されると想像できます)。（編集：numpy.convolveのように、これを行うことができないscipy.signal.convolveを最初に書きました）

秘訣は、この畳み込み関数をその場で使用できるため、二重の for ループです。

for i in xrange(1, flux.shape[0]):
    for j in xrange(1, flux.shape[1]):
        flux[i,j] = 0.5*(flux[i-1,j] + bz[i-1,j]*dx) + 0.5*(flux[i,j-1] - bx[i,j-1]*dz)

次のように置き換えることができます (非常に多くの一時配列が必要なため、申し訳ありません...):

from scipy.ndimage import convolve
_flux = np.zeros((flux.shape[0]+1, flux.shape[1]+1), dtype=flux.dtype)
temp_bx = np.zeros((bx.shape[0]+1, bx.shape[1]+1), dtype=bx.dtype)
temp_bz = np.zeros((bz.shape[0]+1, bz.shape[1]+1), dtype=bz.dtype)

_flux[:-1,:-1] = flux
convolve(_flux[:-1,:-1], [[0, 0.5], [0.5, 0]], _flux[1:,1:])

temp_bz[1:,1:-1] = bz[:,1:]*dx
temp_bx[1:-1,1:] = bx[1:,:]*dz

conv_b = np.array([[0.0, 0.5], [0.5, 0.5]])
convolve(temp_bz[:-1,:-1], [[0.5, 0.5], [0.5, 0.]], temp_bz[1:,1:])
convolve(temp_bx[:-1,:-1], [[-0.5, 0.5], [0.5, 0.]], temp_bx[1:,1:])

flux = _flux[:-1,:-1] + temp_by[:-1,:-1] + temp_bx[:-1,:-1]

残念ながら、bx、bz が最終結果にどのように入るかをいじる必要があることを意味しますが、このアプローチは 2 の大きな累乗を作成することを回避し、前の回答よりも著しく高速になるはずです。

(numpy convolve 関数では、このインプレースでの使用が許可されていないことに注意してください。)

于 2012-08-08T19:26:53.777 に答える
0

convolve を使用する方法は非常に優れていますが、ステンシルの実行方法は明らかではありません...この行の場合

flux[i,j] = 0.5*(flux[i-1,j] + bz[i-1,j]*dx) + 0.5*(flux[i,j-1] - bx[i,j-1]*dz)

に置き換えられました

flux[i,j] = a*flux[i-1,j] + b*bz[i-1,j] + c*flux[i,j-1] - d*bx[i,j-1]

最初の畳み込み (_flux で) はステンシル [[0, a], [b, 0]] を使用すると思います。しかし、a、b、c、および d がスカラーであると仮定すると、bz および bx の他の 2 つのステンシルは何でしょうか?

于 2016-01-20T16:26:17.430 に答える