ループの代わりにnumpyのベクトル化

Question

私はPythonでいくつかのコードを書きましたが、うまく動作しますが非常に遅いです。私はforループのためだと思います。 numpyコマンドを使用して以下の操作を高速化できることを願っています。私は目標を定義しましょう。

次元数がrow x colの2D numpy配列があるとします。例えば、6 x 11アレイ（下記の図を参照）を考えてみましょう。配列が得られ

1. 私はすべての行の平均値を計算したい、すなわち合計ⱼaᵢⱼ。これはもちろん簡単に行うことができます。 各行ため、私はそれらの合計を計算し、すべての列の数で割って、いくつかの選択された値、特定の閾値以下、すなわち、すべての値の平均を計算する、今

2. （私はこの値CM_tildeを呼び出します） （N）。値がこの定義済みのしきい値を超える場合は、CM_tildeの値（行全体の平均）が加算されます。この値は、その後CM

3. 呼ばれ、CM値は、これに加えて、行

内の各要素から減算され、私はすべてのそれらのCMの値が記載されているnumpyの配列またはリストが欲しいです。

フィギュア：

次のコードが機能していますが、非常に遅い（配列が大きい取得する場合は特に）

CM_tilde = np.mean(data, axis=1) 
N = data.shape[1] 
data_cm = np.zeros((data.shape[0], data.shape[1], data.shape[2])) 
all_CMs = np.zeros((data.shape[0], data.shape[2])) 
for frame in range(data.shape[2]): 
    for row in range(data.shape[0]): 
     CM=0 
     for col in range(data.shape[1]): 
      if data[row, col, frame] < (CM_tilde[row, frame]+threshold): 
       CM += data[row, col, frame] 
      else: 
       CM += CM_tilde[row, frame] 
     CM = CM/N 
     all_CMs[row, frame] = CM 
     # calculate CM corrected value 
     for col in range(data.shape[1]): 
      data_cm[row, col, frame] = data[row, col, frame] - CM 
    print "frame: ", frame 
return data_cm, all_CMs 

任意のアイデア？

Answer 1

それは何をやっているベクトル化するのは非常に簡単です。あなたのオリジナルでこれを比較

import numpy as np 

#generate dummy data 
nrows=6 
ncols=11 
nframes=3 
threshold=0.3 
data=np.random.rand(nrows,ncols,nframes) 

CM_tilde = np.mean(data, axis=1) 
N = data.shape[1] 

all_CMs2 = np.mean(np.where(data < (CM_tilde[:,None,:]+threshold),data,CM_tilde[:,None,:]),axis=1) 
data_cm2 = data - all_CMs2[:,None,:] 

：

In [684]: (data_cm==data_cm2).all() 
Out[684]: True 

In [685]: (all_CMs==all_CMs2).all() 
Out[685]: True 

ロジックは、我々は同時にサイズ[nrows,ncols,nframes]の配列を扱うということです。主なやり方は、サイズ[nrows,nframes]のCM_tildeをサイズ[nrows,1,nframes]のCM_tilde[:,None,:]に変更して、Pythonの放送を利用することです。 Pythonは、この変更されたCM_tildeのシングルトンディメンションであるため、各列に同じ値を使用します。 、我々は再び、我々はdataの対応する値を取得したいかどうか（thresholdに基づく）を選択、またはCM_tildeの放送値をnp.whereを使用することにより

。 np.meanを新たに使用すると、all_CMs2を計算することができます。

最後のステップでは、dataの対応する要素からこの新しいall_CMs2を直接差し引いて放送を利用しました。

一時変数の暗黙的なインデックスを調べることで、このようにコードをベクトル化するのに役立ちます。私の意味は、一時変数CMがループ内で[nrows,nframes]以上に存在し、その値が繰り返しごとにリセットされるということです。これはCMが実質的に数量CM[row,frame]（後で2d配列all_CMsに明示的に割り当てられる）であることを意味し、ここから適切なCMtmp[row,col,frames]量をその列次元に沿って合計することによってそれを構築できることが容易にわかります。役に立つ場合は、np.where(...)の部分をCMtmpとし、それからnp.mean(CMtmp,axis=1)を計算することができます。明らかに同じ結果ですが、おそらくより透明性があります。

Answer 2

ここにあなたの関数のベクトル化があります。私は内部から作業し、私が行ったように以前のバージョンをコメントアウトしました。したがって、私がベクトル化した最初のループには###のコメントがあります。

これはきれいではありません。@Andras's答えがうまくいきますが、うまくいけば、この問題を段階的にどのように解決できるかを教えてください。

def foo2(data, threshold): 
    CM_tilde = np.mean(data, axis=1) 
    N = data.shape[1] 
    #data_cm = np.zeros((data.shape[0], data.shape[1], data.shape[2])) 
    ##all_CMs = np.zeros((data.shape[0], data.shape[2])) 
    bmask = data < (CM_tilde[:,None,:] + threshold) 
    CM = np.zeros_like(data) 
    CM[:] = CM_tilde[:,None,:] 
    CM[bmask] = data[bmask] 
    CM = CM.sum(axis=1) 
    CM = CM/N 
    all_CMs = CM.copy() 
    """ 
    for frame in range(data.shape[2]): 
     for row in range(data.shape[0]): 
      ###print(frame, row) 
      ###mask = data[row, :, frame] < (CM_tilde[row, frame]+threshold) 
      ###print(mask) 
      ##mask = bmask[row,:,frame] 
      ##CM = data[row, mask, frame].sum() 
      ##CM += (CM_tilde[row, frame]*(~mask)).sum() 

      ##CM = CM/N 
      ##all_CMs[row, frame] = CM 
      ## calculate CM corrected value 
      #for col in range(data.shape[1]): 
      # data_cm[row, col, frame] = data[row, col, frame] - CM[row,frame] 
     print "frame: ", frame 
    """ 
    data_cm = data - CM[:,None,:] 
    return data_cm, all_CMs 

出力がこの小さなテストケースに一致しています。これは何よりも正しい次元を得るのに役立ちました。ステップ2では

threshold = .1 
data = np.arange(4*3*2,dtype=float).reshape(4,3,2)