2017-10-11 60 views
3

我有一個3D numpy數組input_data(q x m x n),我正在使用它來創建直方圖數據以最終繪製,它存儲在plot_data(m x n x 2)中。這一步在我的流程中是一個體面的瓶頸,我想知道是否有更快,更「粗糙」的做法。有沒有辦法使用numpy去除循環?

num_bins = 3 
for i in range(m): 

    for j in range(n): 

     data = input_data[:, i, j] 

     hist, bins = np.histogram(data, bins=num_bins) 

     # Create the (x, y) pairs to plot 
     plot_data[i][j] = np.stack((bins[:-1], hist), axis=1) 
+0

更快你的意思是更簡潔嗎? –

+1

更快的意義運行時間,所以理想地利用numpy的矢量化能力,類似於使用np.sum()計算總和而不是循環並手動計算 – holtc

+0

因此,我查找了一些信息。也許,這個scipy網站會是你想要的?我發現以下內容: https://docs.scipy.org/doc/numpy-1.13.0/reference/arrays.nditer.html – Maderas

回答

0

我認爲你的樣本是相似的,所以直方圖是相似的。在這種情況下,您可以簡化比較,做一個更加量化的方式:

a=np.random.rand(100000,10,10) 


def f(): # roughly your approach. 
    plotdata=np.zeros((10,10,3),np.int32) 
    for i in range(10): 
     for j in range(10): 
      bins,hist=np.histogram(a[:,i,j],3) 
      plotdata[i,j]=bins 
    return plotdata 

def g(): #vectored comparisons 
    u=(a < 1/3).sum(axis=0) 
    w=(a > 2/3).sum(axis=0) 
    v=len(a)-u-w 
    return np.dstack((u,v,w)) 

對於一個8倍的改進:

In [213]: %timeit f() 
548 ms ± 15.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) 

In [214]: %timeit g() 
77.7 ms ± 5.46 ms per loop (mean ± std. dev. of 7 runs, 10 loops each) 
+0

不幸的是,我們不能這樣簡化,我們需要真實的數據 – holtc

1

下面是箱的通用號碼一個量化的方法 -

def vectorized_app(input_data, num_bins): 
    s0 = input_data.min(0) 
    s1 = input_data.max(0) 

    m,n,r = input_data.shape 
    ids = (num_bins*((input_data - s0)/(s1-s0))).astype(int).clip(max=num_bins-1) 
    offset = num_bins*(r*np.arange(n)[:,None] + np.arange(r)) 
    ids3D = ids + offset 
    count3D = np.bincount(ids3D.ravel(), minlength=n*r*num_bins).reshape(n,r,-1) 
    bins3D = create_ranges_nd(s0, s1, num_bins+1)[...,:-1] 

    out = np.empty((n,r,num_bins,2)) 
    out[...,0] = bins3D 
    out[...,1] = count3D 
    return out 

輔助功能(s) -

# https://stackoverflow.com/a/46694364/ @Divakar 
def create_ranges_nd(start, stop, N, endpoint=True): 
    if endpoint==1: 
     divisor = N-1 
    else: 
     divisor = N 
    steps = (1.0/divisor) * (stop - start) 
    return start[...,None] + steps[...,None]*np.arange(N) 

運行測試

原始的方法 -

def org_app(input_data, num_bins): 
    q,m,n = input_data.shape 
    plot_data = np.zeros((m,n,num_bins,2)) 
    for i in range(m): 
     for j in range(n): 
      data = input_data[:, i, j] 
      hist, bins = np.histogram(data, bins=num_bins) 
      plot_data[i][j] = np.stack((bins[:-1], hist), axis=1) 
    return plot_data 

時序和驗證 -

讓我們測試出形狀(100, 100, 100)的大型數據陣列上,並與框的數目爲10

In [967]: # Setup input 
    ...: num_bins = 10 
    ...: m = 100 
    ...: n = 100 
    ...: q = 100 
    ...: input_data = np.random.rand(q,m,n) 
    ...: 
    ...: out1 = org_app(input_data, num_bins) 
    ...: out2 = vectorized_app(input_data, num_bins) 
    ...: print np.allclose(out1, out2) 
    ...: 
True 

In [968]: %timeit org_app(input_data, num_bins) 
1 loop, best of 3: 748 ms per loop 

In [969]: %timeit vectorized_app(input_data, num_bins) 
100 loops, best of 3: 12.7 ms per loop 

In [970]: 748/12.7 # speedup with vectorized one over original 
Out[970]: 58.89763779527559 
+0

這與在大型3d陣列上的for循環相比,性能如何? – holtc

+0

@holtc增加了假設隨機數的形狀(100,100,100)和箱數= 10。 – Divakar

+0

因此,當我在我的真實數據上運行這個時,np.allclose()是錯誤的,因爲我正在碰撞所有我的方法,實際上速度較慢,所以不幸我無法使用此功能。在我的例子中,input_data的數量是(100000,8,360),最多100個bin – holtc

相關問題