2

ベクトル化により、以下のコードを高速化しようとしています。

[rows,cols] = flow_direction_np.shape
elevation_gain = np.zeros((rows,cols), np.float)

for [i, j], flow in np.ndenumerate(flow_direction_np):
    try:
        if flow == 32:
            elevation_gain[i - 1, j - 1]  = elevation_gain[i - 1, j - 1] + sediment_transport_np[i, j]
        elif flow == 64:
            elevation_gain[i - 1, j]  = elevation_gain[i - 1, j] + sediment_transport_np[i, j]
        elif flow == 128:
            elevation_gain[i - 1, j + 1]  = elevation_gain[i - 1, j + 1] + sediment_transport_np[i, j]
        elif flow == 16:
            elevation_gain[i, j - 1]  = elevation_gain[i, j - 1] + sediment_transport_np[i, j]
        elif flow == 1:
            elevation_gain[i, j + 1]  = elevation_gain[i, j + 1] + sediment_transport_np[i, j]
        elif flow == 2:
            elevation_gain[i + 1, j + 1]  = elevation_gain[i + 1, j + 1] + sediment_transport_np[i, j]
        elif flow == 4:
            elevation_gain[i + 1, j]  = elevation_gain[i + 1, j] + sediment_transport_np[i, j]
        elif flow == 8:
            elevation_gain[i + 1, j - 1]  = elevation_gain[i + 1, j - 1] + sediment_transport_np[i, j]
    except IndexError:
            elevation_gain[i, j] = 0

これは私のコードが現時点でどのように見えるかです:

elevation_gain = np.zeros_like(sediment_transport_np)
nrows, ncols = flow_direction_np.shape
lookup = {32: (-1, -1),
            16:  (0, -1), 
            8:   (+1, -1),
            4:   (+1,  0),
            64: (-1,  0),
            128:(-1,  +1),
            1:   (0,  +1),
            2:   (+1,  +1)}

# Initialize an array for the "shifted" mask
shifted = np.zeros((nrows+2, ncols+2), dtype=bool)

# Pad elevation gain with zeros
tmp = np.zeros((nrows+2, ncols+2), elevation_gain.dtype)
tmp[1:-1, 1:-1] = elevation_gain
elevation_gain = tmp

for value, (row, col) in lookup.iteritems():
    mask = flow_direction_np == value

    # Reset the "shifted" mask
    shifted.fill(False)
    shifted[1:-1, 1:-1] = mask

    # Shift the mask by the right amount for the given value
    shifted = np.roll(shifted, row, 0)
    shifted = np.roll(shifted, col, 1)

    # Set the values in elevation change to the offset value in sed_trans
    elevation_gain[shifted] = elevation_gain[shifted] + sediment_transport_np[mask]

私が抱えている問題は、彼らが最後に同じ結果を私に与えていないということです。

4

2 に答える 2

0

異なる結果が得られる理由は、python が負のインデックスを処理する方法によるものです。

読んでいる他の人のために、この質問（および回答）はここからのフォローアップです：numpy array を反復し、別の配列の値にインデックスを付けます

まず、「ベクトル化された」コードが非常に密集していることをお詫びします。以前の回答に詳しい説明があるので、ここでは繰り返しません。

元のコード (元の質問) は、実際には、ここに投稿したバージョンとは微妙に異なります。

基本的に、あなたが持っている前に

for [i, j], flow in np.ndenumerate(flow_direction_np):
    try:
        if flow == 32:
            ...
        elif ...
            ...

i+1またはj+1がグリッドのサイズより大きい場合、インデックスエラーが発生していました。

ただやっている：

for [i, j], flow in np.ndenumerate(flow_direction_np):
    try:
        if flow == 32:
            ...
        elif ...
            ...
    except IndexError:
        elevation_change[i, j] = 0

グリッドの異なる側で異なる境界条件を定義するため、実際には正しくありません。

2 番目のケースでは、j-1ori-1が負の場合、グリッドの反対側の値が返されます。ただし、j+1またはi+1がグリッドのサイズより大きい場合は、0が返されます。(したがって、「異なる境界条件」。)

コードのベクトル化されたバージョンでは、インデックスが負の場合とグリッドの境界を超えている場合の両方0でが返されます。

簡単な例として、次の場合に何が起こるかに注意してください。

In [1]: x = [1, 2, 3]

In [2]: x[0]
Out[2]: 1

In [3]: x[1]
Out[3]: 2

In [4]: x[2]
Out[4]: 3

In [5]: x[3]
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-5-ed224ad0520d> in <module>()
----> 1 x[3]

IndexError: list index out of range

In [6]: x[-1]
Out[6]: 3

In [7]: x[-2]
Out[7]: 2

In [8]: x[-3]
Out[8]: 1

In [9]: x[-4]
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-9-f9c639f21256> in <module>()
----> 1 x[-4]

IndexError: list index out of range

In [10]:

シーケンスのサイズまでの負のインデックスは有効であり、シーケンスの「反対側」を返すことに注意してください。したがって、x[3]エラーが発生しx[-1]ますが、もう一方の端が返されます。

うまくいけば、それは少し明確です。

于 2013-07-19T14:26:23.013 に答える
0

np.where条件が発生するインデックスを取得するために使用すると、パフォーマンスを大幅に向上させることができます。

ind = np.where( flow_direction_np==32 )

2 つの要素を持つタプルであることがわかります。1 つ目は配列indの最初の軸のインデックスで、2 つ目は 2 番目の軸のインデックスです。flow_direction_np

このインデックスを使用してシフトを適用できます:i-1などj-1...

ind_32 = (ind[0]-1, ind[1]-1)

次に、ファンシー インデックスを使用して配列を更新します。

elevation_gain[ ind_32 ] += sediment_transport_np[ ind ]

編集:この概念をあなたのケースに適用すると、次のようになります:

lookup = {32: (-1, -1),
          16: ( 0, -1),
           8: (+1, -1),
           4: (+1,  0),
          64: (-1,  0),
         128: (-1, +1),
           1: ( 0, +1),
           2: (+1, +1)}

for num, shift in lookup.iteritems():
    ind = np.where( flow_direction_np==num )
    ind_num = ind[0] + shift[0], ind[1] + shift[1]
    elevation_gain[ ind_num] += sediment_transport_np[ ind ]
于 2013-07-19T12:16:21.373 に答える