6 回答

TA貢獻1155條經驗 獲得超0個贊
無循環解決方案是max在由以下命令創建的窗口上使用skimage.util.view_as_windows:
list(map(max, view_as_windows(A, (2,))))
[8, 33, 33, 4, 6]
復制/粘貼示例:
import numpy as np
from skimage.util import view_as_windows
A = np.array([8, 2, 33, 4, 3, 6])
list(map(max, view_as_windows(A, (2,))))

TA貢獻1796條經驗 獲得超7個贊
在本次問答中,我們基本上要求滑動最大值。之前已經探討過這一點 -?NumPy 數組中滑動窗口中的 Max。由于我們希望提高效率,因此我們可以看得更遠。其中之一是numba
,這里有兩個最終變體,我最終得到了杠桿parallel
指令,它比沒有版本提高了性能:
import numpy as np
from numba import njit, prange
@njit(parallel=True)
def numba1(a, W):
? ? L = len(a)-W+1
? ? out = np.empty(L, dtype=a.dtype)
? ? v = np.iinfo(a.dtype).min
? ? for i in prange(L):
? ? ? ? max1 = v
? ? ? ? for j in range(W):
? ? ? ? ? ? cur = a[i + j]
? ? ? ? ? ? if cur>max1:
? ? ? ? ? ? ? ? max1 = cur? ? ? ? ? ? ? ??
? ? ? ? out[i] = max1
? ? return out?
@njit(parallel=True)
def numba2(a, W):
? ? L = len(a)-W+1
? ? out = np.empty(L, dtype=a.dtype)
? ? for i in prange(L):
? ? ? ? for j in range(W):
? ? ? ? ? ? cur = a[i + j]
? ? ? ? ? ? if cur>out[i]:
? ? ? ? ? ? ? ? out[i] = cur? ? ? ? ? ? ? ??
? ? return out?
從之前鏈接的問答中,等效的 SciPy 版本將是 -
from scipy.ndimage.filters import maximum_filter1d
def scipy_max_filter1d(a, W):
? ? L = len(a)-W+1
? ? hW = W//2 # Half window size
? ? return maximum_filter1d(a,size=W)[hW:hW+L]
標桿管理
其他發布的通用窗口 arg 的工作方法:
from skimage.util import view_as_windows
def rolling(a, window):
? ? shape = (a.size - window + 1, window)
? ? strides = (a.itemsize, a.itemsize)
? ? return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)
# @mathfux's soln
def npmax_strided(a,n):
? ? return np.max(rolling(a, n), axis=1)
# @Nicolas Gervais's soln
def mapmax_strided(a, W):
? ? return list(map(max, view_as_windows(a,W)))
cummax = np.maximum.accumulate
def pp(a,w):
? ? N = a.size//w
? ? if a.size-w+1 > N*w:
? ? ? ? out = np.empty(a.size-w+1,a.dtype)
? ? ? ? out[:-1] = cummax(a[w*N-1::-1].reshape(N,w),axis=1).ravel()[:w-a.size-1:-1]
? ? ? ? out[-1] = a[w*N:].max()
? ? else:
? ? ? ? out = cummax(a[w*N-1::-1].reshape(N,w),axis=1).ravel()[:w-a.size-2:-1]
? ? out[1:N*w-w+1] = np.maximum(out[1:N*w-w+1],
? ? ? ? ? ? ? ? ? ? ? ? ? ? cummax(a[w:w*N].reshape(N-1,w),axis=1).ravel())
? ? out[N*w-w+1:] = np.maximum(out[N*w-w+1:],cummax(a[N*w:]))
? ? return out
使用benchit包(很少有基準測試工具打包在一起;免責聲明:我是它的作者)對建議的解決方案進行基準測試。
import benchit
funcs = [mapmax_strided, npmax_strided, numba1, numba2, scipy_max_filter1d, pp]
in_ = {(n,W):(np.random.randint(0,100,n),W) for n in 10**np.arange(2,6) for W in [2, 10, 20, 50, 100]}
t = benchit.timings(funcs, in_, multivar=True, input_name=['Array-length', 'Window-length'])
t.plot(logx=True, sp_ncols=1, save='timings.png')
因此,numba 非常適合小于 的窗口大小10
,此時沒有明顯的贏家,而在較大的窗口大小上,pp
SciPy 則排名第二。

TA貢獻1839條經驗 獲得超15個贊
這是專門為較大窗戶設計的方法。窗口大小為 O(1),數據大小為 O(n)。
我已經完成了純 numpy 和 pythran 實現。
我們如何在窗口大小上實現 O(1) ?我們使用“鋸齒”技巧:如果 w 是窗口寬度,我們將數據分為很多 w ,對于每個組,我們從左到右和從右到左求累積最大值。任何中間窗口的元素分布在兩組中,并且交集的最大值位于我們之前計算的累積最大值之中。因此,每個數據點總共需要 3 次比較。
benchit w=100;我的函數是 pp (numpy) 和 winmax (pythran):
對于小窗口尺寸 w=5,圖像更加均勻。有趣的是,即使對于非常小的尺寸,pythran 仍然具有巨大的優勢。他們必須采取正確的措施來最大限度地減少呼叫開銷。
蟒蛇代碼:
cummax = np.maximum.accumulate
def pp(a,w):
? ? N = a.size//w
? ? if a.size-w+1 > N*w:
? ? ? ? out = np.empty(a.size-w+1,a.dtype)
? ? ? ? out[:-1] = cummax(a[w*N-1::-1].reshape(N,w),axis=1).ravel()[:w-a.size-1:-1]
? ? ? ? out[-1] = a[w*N:].max()
? ? else:
? ? ? ? out = cummax(a[w*N-1::-1].reshape(N,w),axis=1).ravel()[:w-a.size-2:-1]
? ? out[1:N*w-w+1] = np.maximum(out[1:N*w-w+1],
? ? ? ? ? ? ? ? ? ? ? ? ? ? cummax(a[w:w*N].reshape(N-1,w),axis=1).ravel())
? ? out[N*w-w+1:] = np.maximum(out[N*w-w+1:],cummax(a[N*w:]))
? ? return out
pythran 版本;編譯用pythran -O3 <filename.py>; 這將創建一個可以導入的已編譯模塊:
import numpy as np
# pythran export winmax(float[:],int)
# pythran export winmax(int[:],int)
def winmax(data,winsz):
? ? N = data.size//winsz
? ? if N < 1:
? ? ? ? raise ValueError
? ? out = np.empty(data.size-winsz+1,data.dtype)
? ? nxt = winsz
? ? for j in range(winsz,data.size):
? ? ? ? if j == nxt:
? ? ? ? ? ? nxt += winsz
? ? ? ? ? ? out[j+1-winsz] = data[j]
? ? ? ? else:
? ? ? ? ? ? out[j+1-winsz] = out[j-winsz] if out[j-winsz]>data[j] else data[j]
? ? running = data[-winsz:N*winsz].max()
? ? nxt -= winsz << (nxt > data.size)
? ? for j in range(data.size-winsz,0,-1):
? ? ? ? if j == nxt:
? ? ? ? ? ? nxt -= winsz
? ? ? ? ? ? running = data[j-1]
? ? ? ? else:
? ? ? ? ? ? running = data[j] if data[j] > running else running
? ? ? ? ? ? out[j] = out[j] if out[j] > running else running
? ? out[0] = data[0] if data[0] > running else running
? ? return out

TA貢獻1815條經驗 獲得超6個贊
如果存在連續的n項目,則擴展解決方案需要循環:
np.maximum(*[A[i:len(A)-n+i+1] for i in range(n)])
為了避免這種情況,您可以使用跨步技巧并將其轉換A為n-length 塊的數組:
def rolling(a, window):
shape = (a.size - window + 1, window)
strides = (a.itemsize, a.itemsize)
return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)
例如:
>>> rolling(A, 3)
array([[ 8, 2, 8],
[ 2, 8, 33],
[ 8, 33, 33],
[33, 33, 4]])
完成后你可以用 殺死它np.max(rolling(A, n), axis=1)。
盡管它很優雅,但這個解決方案和第一個解決方案都不是高效的,因為我們在僅相差兩項的相鄰塊上重復應用最大值。

TA貢獻1830條經驗 獲得超3個贊
對于所有 n 的遞歸解決方案
import numpy as np
import sys
def recursive(a: np.ndarray, n: int, b=None, level=2):
if n <= 0 or n > len(a):
raise ValueError(f'len(a):{len(a)} n:{n}')
if n == 1:
return a
if len(a) == n:
return np.max(a)
b = np.maximum(a[:-1], a[1:]) if b is None else np.maximum(a[level - 1:], b)
if n == level:
return b
return recursive(a, n, b[:-1], level + 1)
test_data = np.array([8, 2, 33, 4, 3, 6])
for test_n in range(1, len(test_data) + 2):
try:
print(recursive(test_data, n=test_n))
except ValueError as e:
sys.stderr.write(str(e))
輸出
[ 8 2 33 4 3 6]
[ 8 33 33 4 6]
[33 33 33 6]
[33 33 33]
[33 33]
33
len(a):6 n:7
關于遞歸函數
你可以觀察下面的數據,然后你就會知道如何寫遞歸函數了。
"""
np.array([8, 2, 33, 4, 3, 6])
n=2: (8, 2), (2, 33), (33, 4), (4, 3), (3, 6) => [8, 33, 33, 4, 6] => B' = [8, 33, 33, 4]
n=3: (8, 2, 33), (2, 33, 4), (33, 4, 3), (4, 3, 6) => B' [33, 4, 3, 6] => np.maximum([8, 33, 33, 4], [33, 4, 3, 6]) => 33, 33, 33, 6
...
"""
添加回答
舉報