在Python/Numpy中使用不同數組形狀的向量化三元組

我是Python/Numpy中的新成員，並且正在嘗試將我的三重循環改進爲更高效的計算，但無法安靜地弄清楚如何執行此操作。在Python/Numpy中使用不同數組形狀的向量化三元組

上的大小（25,35）和陣列的形狀的網格進行的計算是：

A = (8760,25,35) 
B = (12,25,35)

第一尺寸A對應於在一年（數小時〜8760 ），B中的第一維是月數（12）。我想使用B [0，：，]中的值作爲第一個月，而B [1，：，]代表第二個等等。

到目前爲止，我以非常不確定的方式創建了一個指數陣列充滿1,1,1 ...... 2,2,2 ......，12從B.我的代碼一些隨機數

N,M = (25, 35) 
    A = np.random.rand(8760,N,M) 
    B = np.random.rand(12,N,M) 
    q = len(A)/12 
    index = np.hstack((np.full((1,q),1),np.full((1,q),2),np.full((1,q),3),np.full((1,q),4),np.full((1,q),5),np.full((1,q),6),np.full((1,q),7),np.full((1,q),8),np.full((1,q),9),np.full((1,q),10),np.full((1,q),11),np.full((1,q),12)))-1 
    results = np.zeros((len(A),N,M)) 

    for t in xrange(len(A)): 
     for i in xrange(N): 
      for j in xrange(M): 
       results[t][i][j] = some_function(A[t][i][j], B[index[(0,t)]][i][j],H = 80.) 


    def some_function(A,B,H = 80.0): 
     results = A*np.log(H/B)/np.log(10./B) 
     return results

提取值怎樣才能提高這個速度如何計算？

來源

2015-09-24 HolmN

您能否分享一下'some_function'細節？因爲對於真正的矢量化，嵌套循環內的*深*嵌入的函數必須*出來循環，這取決於它的實現。 – Divakar

我編輯了包含該功能的問題！我感謝您的幫助！ – HolmN

由於index [（0，t）]是一個標量，所以'index'是一維數組會在'index [（0，t）] [i] [j]'處出錯。 – Divakar

NumPy支持broadcasting，它允許以高度優化的方式在不同形狀的陣列上執行元素操作。在你的情況下，你的行數和列數在A和B是一樣的。但是，在第一維中，這兩個數組的元素數量是不同的。看看實現，看起來B的第一維元素按q編號重複，直到我們轉到第一維中的下一個元素。這與B的第一維中的元素的數量是第一維中的元素的數量q乘以A的事實相符。現在

，回到broadcasting，該解決方案將是分裂的A第一尺寸具有四維陣列，例如，我們在此第一維的元素的數量重塑四維陣列匹配了B的第一維中元素的數量。接下來，reshapeB也可以通過在B[:,None,:,:]的第二維上創建singleton維度（無元素的維度）。然後，NumPy將使用廣播魔術並執行廣播元素乘法，因爲這正是我們在我們的some_function中所做的。

下面是一個使用NumPy's broadcasting能力的量化執行 -

H = 80.0 
M,N,R = B.shape 
B4D = B[:,None,:,:] 
out = ((A.reshape(M,-1,N,R)*np.log(H/B4D))/np.log(10./B4D)).reshape(-1,N,R)

運行測試和驗證輸出 -

In [50]: N,M = (25, 35) 
    ...: A = np.random.rand(8760,N,M) 
    ...: B = np.random.rand(12,N,M) 
    ...: H = 80.0 
    ...: 

In [51]: def some_function(A,B,H = 80.0): 
    ...:  return A*np.log(H/B)/np.log(10./B) 
    ...: 

In [52]: def org_app(A,B,H): 
    ...: q = len(A)/len(B) 
    ...: index = np.repeat(np.arange(len(B))[:,None],q,axis=1).ravel()[None,:] # Simpler 
    ...: results = np.zeros((len(A),N,M)) 
    ...: for t in xrange(len(A)): 
    ...:  for i in xrange(N): 
    ...:   for j in xrange(M): 
    ...:    results[t][i][j] = some_function(A[t][i][j], B[index[(0,t)]][i][j]) 
    ...: return results 
    ...: 

In [53]: def vectorized_app(A,B,H): 
    ...: M,N,R = B.shape 
    ...: B4D = B[:,None,:,:] 
    ...: return ((A.reshape(M,-1,N,R)*np.log(H/B4D))/np.log(10./B4D)).reshape(-1,N,R) 
    ...: 

In [54]: np.allclose(org_app(A,B,H), vectorized_app(A,B,H)) 
Out[54]: True 

In [55]: %timeit org_app(A,B,H) 
1 loops, best of 3: 1min 32s per loop 

In [56]: %timeit vectorized_app(A,B,H) 
10 loops, best of 3: 217 ms per loop

來源

2015-09-24 09:30:17 Divakar

在Python/Numpy中使用不同數組形狀的向量化三元組

回答

相關問題