2017-02-12 10 views
0

私はPythonを使って気候モデルの出力を研究して、特定のタイプの嵐を見つける研究者です。私は8つの大きなnumpy配列(次元は109574 x 52 x 57)です。これらの配列は、その日に嵐があったことを示すために1で埋められています(最初の次元は時間です)。嵐がない場合は0です。他の2つの次元は緯度と経度です。重複をチェックする効率の向上 - Python

これらの配列からバックツーバックの日を削除する必要があります。たとえば、1日目と2日目に嵐があった場合、1回の嵐だけをカウントしたいと思います。 1日目、2日目、3日目に嵐があった場合は、合計2回の嵐のために1と3を数えたいと思っています.1〜4日には2回の嵐があります。 np.sumを使用して最後に嵐の数を発見し、時間軸に沿って配列の1を集計しました。

私はこれを達成するために次のコードを実行していますが、非常に遅いという問題に直面しています。他のデータセットに対してこの手順を繰り返す必要があるので、このプロセスを効率化する方法があるかどうか疑問に思っていました。私は下に私のコードを持っています、そして、私は何かを明確にすること以上に満足しています。

# If there is a storm that overlaps two two-day periods, only count it once 
print("Eliminating doubles...") 
for i in range(52): 
    for j in range(57): 
     print(i,j) 
     for k in range(109573): 
      if((storms1[k,i,j]) == 1 and (storms1[k+1,i,j] == 1)): 
       storms1[k,i,j] = 0 
      if((storms2[k,i,j]) == 1 and (storms2[k+1,i,j] == 1)): 
       storms2[k,i,j] = 0 
      if((storms3[k,i,j]) == 1 and (storms3[k+1,i,j] == 1)): 
       storms3[k,i,j] = 0 
      if((storms4[k,i,j]) == 1 and (storms4[k+1,i,j] == 1)): 
       storms4[k,i,j] = 0 
      if((storms5[k,i,j]) == 1 and (storms5[k+1,i,j] == 1)): 
       storms5[k,i,j] = 0 
      if((storms6[k,i,j]) == 1 and (storms6[k+1,i,j] == 1)): 
       storms6[k,i,j] = 0 
      if((storms7[k,i,j]) == 1 and (storms7[k+1,i,j] == 1)): 
       storms7[k,i,j] = 0 
      if((storms8[k,i,j]) == 1 and (storms8[k+1,i,j] == 1)): 
       storms8[k,i,j] = 0 

誰かがループと配列を反復処理を示唆する前に、私はこの質問をする目的のためにそれらを簡素化するために、変数名を変更しました。

ありがとうございます。ここで

答えて

2

はあなたの最も内側のループを置き換えることができますベクトル化機能である:

def do(KK): 
    # find stretches of ones 
    switch_points = np.where(np.diff(np.r_[0, KK, 0]))[0] 
    switch_points.shape = -1, 2 
    # isolate stretches starting on odd days and create mask 
    odd_starters = switch_points[switch_points[:, 0] % 2 == 1, :] 
    odd_mask = np.zeros((KK.shape[0] + 1,), dtype=KK.dtype) 
    odd_mask[odd_starters] = 1, -1 
    odd_mask = np.add.accumulate(odd_mask[:-1]) 
    # apply global 1,0,1,0,1,0,... mask 
    KK[1::2] = 0 
    # invert stretches starting on odd days 
    KK ^= odd_mask 

は、ループ(iとj)の外側のペアの中からそれを呼び出す:

do(storms1[:, i, j]) 
do(storms2[:, i, j]) 
etc. 

これは、配列を変更します所定の位置に。

これは、ループ処理よりもはるかに高速でなければなりません(2つの外側ループは違いがありません)。それはもののブロックの始点と終点を見つけ

:仕組み

。私たちは、そのようなブロックごとに1つおきにゼロがなければならないことを知っています。 グローバル1,0,1,0,1,0、...マスクを使用すると、アルゴリズムは1日おきに0になります。偶数日に開始ブロックで

  • 正しい結果を生成

  • 奇数の日に開始ブロック
  • 及びブロック内の正しいパターンの補数外側変化なし

アルゴリズムの最後のステップは、これらの奇数開始ブロックを反転させることです。

1

最初の軸をシミュレートする1​​Dアレイを使用した例。まず、1のグループがどこから始まるかを見つけます。次に、各グループの長さを見つけます。あなたは、私は彼らがもはや必要とされた後など、変数名を再利用しないことにより、ここに示したものよりも多くのメモリを節約することができ

import numpy 

a = numpy.random.randint(0,2,20) 

# Add an initial 0 
a1 = numpy.r_[0, a] 

# Mark the start of each group of 1's 
d1 = numpy.diff(a1) > 0 

# Indices of the start of groups of 1's 
w1 = numpy.arange(len(d1))[d1] 

# Length of each group 
cs = numpy.cumsum(a) 
c = numpy.diff(numpy.r_[cs[w1], cs[-1]+1]) 

# Apply the counting logic 
storms = c - c//2 

print(a) 
>>> array([0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1]) 
print(c) 
>>> array([1, 2, 4, 1, 3]) 
print(storms) 
>>> array([1, 1, 2, 1, 2]) 

:最後に、あなたのロジックに基づいてイベントの数を計算します

0

だから私は、あなたがしたいと思う:

storms_in[:,i,j] = [0,0,1,1,0,1,1,1,0,1,0,1,1,1,1,0] 
storms_out[:,i,j]= [0,0,1,0,0,1,0,1,0,1,0,1,0,0,1,0] 

これはあなたのコードサンプルが何をしているのかないですが、あなたがあなたの第二段落でやりたいと言うものです。

これを行うには、次の2つのステップ

def storms_disc(storms): # put the whole array here, boolean-safe 
    z = np.zeros((1,) + storms.shape[1:]) # zero-pads for the ends 
    changes = np.r_[storms.astype('int8') ,z] - np.r_[z, storms.astype('int8')] #find where the weather changes 
    changes=((changes[:-1] == 1) | (changes[1:] == -1)).astype('int8') # reduce dimension 
    return ((np.r_[changes, z] - np.r_[z, changes])[:-1] == 1).astype(storms.dtype) #find the first of successive changes 

これは、プロセス全体をベクトル化する必要があるだろう、とあなたはそれを8回コールする必要があると思います。 ( `` a.view(numpy.bool))8あなたが見ることができます

storms=np.random.randint(0,2,90).reshape(10,3,3) 
storms.T 

array([[[1, 0, 0, 1, 1, 1, 1, 1, 1, 0], 
     [0, 0, 1, 1, 0, 1, 1, 0, 0, 1], 
     [1, 0, 1, 0, 1, 0, 1, 0, 1, 0]], 

     [[0, 0, 0, 1, 0, 1, 0, 0, 0, 0], 
     [0, 1, 0, 0, 1, 1, 1, 0, 0, 0], 
     [0, 1, 0, 0, 1, 0, 1, 0, 1, 1]], 

     [[0, 1, 0, 1, 0, 1, 1, 0, 0, 0], 
     [0, 1, 0, 1, 0, 1, 0, 0, 1, 1], 
     [0, 0, 0, 1, 1, 1, 0, 0, 1, 0]]], dtype=int8) 

storms_disc(storms).T 

array([[[1, 0, 0, 1, 0, 0, 0, 0, 1, 0], 
     [0, 0, 1, 0, 0, 1, 0, 0, 0, 1], 
     [1, 0, 1, 0, 1, 0, 1, 0, 1, 0]], 

     [[0, 0, 0, 1, 0, 1, 0, 0, 0, 0], 
     [0, 1, 0, 0, 1, 0, 1, 0, 0, 0], 
     [0, 1, 0, 0, 1, 0, 1, 0, 1, 0]], 

     [[0, 1, 0, 1, 0, 1, 0, 0, 0, 0], 
     [0, 1, 0, 1, 0, 1, 0, 0, 1, 0], 
     [0, 0, 0, 1, 0, 1, 0, 0, 1, 0]]], dtype=int8) 
+0

注:ブール値を減算すると、エラーが発生するためastype呼び出しはその値が1と0

テストであっても、されています-nbityのブール型も8ビットであるため、int配列をブール値として扱います。これは型変換を節約します。 – Benjamin

+0

どこでそれを行うのか分かりません。私はほとんど 'int8'に変換していますので、私は引き算できます。最後に.astype(storms.dtype)を置き換えることができると思います。 –

+0

「ビュー」トリックは逆も同様です...しかし、そうですね、あなたのコードを少し速く読んでいます。 – Benjamin

関連する問題