2017-09-05 53 views
1

有沒有辦法擺脫下面的代碼中的循環,並用矢量化操作代替它?向量化每個矩陣行的變量範圍的numpy.argmin搜索

給定一個數據矩陣,對於每一行,我想找到適合在單獨數組中定義的範圍內(每行)的最小值的索引。

下面是一個例子:

import numpy as np 
np.random.seed(10) 

# Values of interest, for this example a random 6 x 100 matrix 
data = np.random.random((6,100)) 

# For each row, define an inclusive min/max range 
ranges = np.array([[0.3, 0.4], 
        [0.35, 0.5], 
        [0.45, 0.6], 
        [0.52, 0.65], 
        [0.6, 0.8], 
        [0.75, 0.92]]) 


# For each row, find the index of the minimum value that fits inside the given range 
result = np.zeros(6).astype(np.int) 
for i in xrange(6): 
    ind = np.where((ranges[i][0] <= data[i]) & (data[i] <= ranges[i][1]))[0] 
    result[i] = ind[np.argmin(data[i,ind])] 

print result 
# Result: [35 8 22 8 34 78] 

print data[np.arange(6),result] 
# Result: [ 0.30070006 0.35065639 0.45784951 0.52885388 0.61393513 0.75449247] 
+0

做什麼,如果所有的'data'超出'range'對於一個給定的行? –

回答

1

方法1:使用broadcastingnp.minimum.reduceat -

mask = (ranges[:,None,0] <= data) & (data <= ranges[:,None,1]) 
r,c = np.nonzero(mask) 
cut_idx = np.unique(r, return_index=1)[1] 
out = np.minimum.reduceat(data[mask], cut_idx) 

改進,以避免np.nonzero,並直接從mask計算cut_idx

cut_idx = np.concatenate(([0], np.count_nonzero(mask[:-1],1).cumsum())) 

方法2:使用broadcasting和填充NaNs無效的地方,然後用np.nanargmin -

mask = (ranges[:,None,0] <= data) & (data <= ranges[:,None,1]) 
result = np.nanargmin(np.where(mask, data, np.nan), axis=1) 
out = data[np.arange(6),result] 

方法3:如果不迭代就夠了(就像你有一個6環在樣本中迭代),您可能想要堅持循環以提高內存效率,但要使用更高效的masking而不是布爾數組 -

out = np.zeros(6) 
for i in xrange(6): 
    mask_i = (ranges[i,0] <= data[i]) & (data[i] <= ranges[i,1]) 
    out[i] = np.min(data[i,mask_i]) 

方法#4:這裏還有一個更多的loopy解決方案。這個想法將是對每一行數據進行排序。然後,在np.searchsorted的幫助下,使用每行的兩個範圍限制來決定啓動和停止索引。此外,我們將使用這些指數來切片,然後得到最小值。切片的好處是,我們將與views合作,因此在內存和性能方面都非常高效。

的實施將是這個樣子 -

out = np.zeros(6) 
sdata = np.sort(data, axis=1) 
for i in xrange(6): 
    start = np.searchsorted(sdata[i], ranges[i,0]) 
    stop = np.searchsorted(sdata[i], ranges[i,1], 'right')  
    out[i] = np.min(sdata[i,start:stop]) 

此外,我們可以在一個量化的方式讓那些啓動,停止索引以下的vectorized searchsorted的實現。

基於由@Daniel F因爲當我們面對的是給data限度內的範圍內時建議,我們可以簡單地使用開始指數 -

out[i] = sdata[i, start] 
+0

對於#4:不是'sdata [i,np.searchsorted(sdata [i],ranges [i,0])]'已經是最小值,除非它超出邊界?只要'在哪裏'測試該值並返回'nan',如果它超出範圍。 –

+0

@DanielF不知道我有你。爲什麼'sdata [i,np.searchsorted(sdata [i],ranges [i,0])]'最小?在有序數組'sdata'中,我們正在尋找第一個索引,其中'ranges [i,0]'在它的左邊。這個'stop'需要一個編輯:'np.searchsorted(sdata [i],ranges [i,1],'right')'來覆蓋這樣一個排序數據右邊的第一個索引。 – Divakar

+0

因爲如果'sdata'被排序,'min(sdata [start:stop])'總是'sdata [start]'。實際上,如果所有值都低於下限,那麼現在你會得到一個錯誤,因爲'start'和'stop'將會是'sdata.shape [1]' –

1

在範圍假設至少一個值,你甚至不用與上限打擾:

result = np.empty(6) 
for i in xrange(6): 
    lt = (ranges[i,0] >= data[i]).sum() 
    result[i] = np.argpartition(data[i], lt)[lt] 

事實上,你甚至可以矢量化整個事情用argpartition

lt = (ranges[:,None,0] >= data).sum(1) 
result = np.argpartition(data, lt)[np.arange(data.shape[0]), lt] 

當然,這僅僅是有效的,如果data.shape[0] < < data.shape[1],否則你基本上排序

+0

我沒有機會去尋找問題,但是當運行你的第一個例子時,我在這一行得到一個索引超出界限的錯誤:'result [i] = np.argpartition(data,lt)[lt]'。 ..第二個例子工作。 – Fnord

+0

是的,修正了這個問題。 –