在一维NumPy数组中寻找局部极大值/极小值的奇点/集合(再次)

完整的矢量化解决方案:

test03 = np.array([2,2,10,4,4,4,5,6,7,2,6,5,5,7,7,1,1])  # Size 17
extended = np.empty(len(test03)+2)  # Rooms to manage edges, size 19
extended[1:-1] = test03
extended[0] = extended[-1] = np.inf

flag_left = extended[:-1] <= extended[1:]  # Less than successor, size 18
flag_right = extended[1:] <= extended[:-1]  # Less than predecessor, size 18

flagmini = flag_left[1:] & flag_right[:-1]  # Local minimum, size 17
mini = np.where(flagmini)[0]  # Indices of minimums
spl = np.where(np.diff(mini)>1)[0]+1  # Places to split
result = np.split(mini, spl)

result :

[0, 1] [3, 4, 5] [9] [11, 12] [15, 16]

编辑

不幸的是,当它们至少有3个项目大时,这也会检测到最大值,因为它们被视为平坦的局部最小值。这样的话,一块裸体补丁会很难看。

为了解决这个问题,我提出了另外两种解决方案,先是numpy,然后是numba。

numpy使用的是什么 np.diff :

import numpy as np
test03=np.array([12,13,12,4,4,4,5,6,7,2,6,5,5,7,7,17,17])
extended=np.full(len(test03)+2,np.inf)
extended[1:-1]=test03

slope = np.sign(np.diff(extended))  # 1 if ascending,0 if flat, -1 if descending
not_flat,= slope.nonzero() # Indices where data is not flat.   
local_min_inds, = np.where(np.diff(slope[not_flat])==2) 

#local_min_inds contains indices in not_flat of beginning of local mins. 
#Indices of End of local mins are shift by +1:   
start = not_flat[local_min_inds]
stop =  not_flat[local_min_inds+1]-1

print(*zip(start,stop))
#(0, 1) (3, 5) (9, 9) (11, 12) (15, 16)    

与numba acceleration兼容的直接解决方案:

#@numba.njit
def localmins(a):
    begin= np.empty(a.size//2+1,np.int32)
    end  = np.empty(a.size//2+1,np.int32)
    i=k=0
    begin[k]=0
    search_end=True
    while i<a.size-1:
         if a[i]>a[i+1]:
                begin[k]=i+1
                search_end=True
         if search_end and a[i]<a[i+1]:   
                end[k]=i
                k+=1
                search_end=False
        i+=1
    if search_end and i>0  : # Final plate if exists 
        end[k]=i
        k+=1 
    return begin[:k],end[:k]

    print(*zip(*localmins(test03)))
    #(0, 1) (3, 5) (9, 9) (11, 12) (15, 16)  

我认为另一个功能来自 scipy.signal 会很有趣的。

from scipy.signal import find_peaks

test03 = np.array([2,2,10,4,4,4,5,6,7,2,6,5,5,7,7,1,1])
find_peaks(test03)

Out[]: (array([ 2,  8, 10, 13], dtype=int64), {})

find_peaks 有很多选择,可能非常有用,尤其是对于嘈杂的信号。

使现代化

这个功能非常强大,用途广泛。可以为峰值最小宽度、高度、彼此之间的距离等设置多个参数。例如:

test04 = np.array([1,1,5,5,5,5,5,5,5,5,1,1,1,1,1,5,5,5,1,5,1,5,1])
find_peaks(test04, width=1)

Out[]: 
(array([ 5, 16, 19, 21], dtype=int64),
 {'prominences': array([4., 4., 4., 4.]),
  'left_bases': array([ 1, 14, 18, 20], dtype=int64),
  'right_bases': array([10, 18, 20, 22], dtype=int64),
  'widths': array([8., 3., 1., 1.]),
  'width_heights': array([3., 3., 3., 3.]),
  'left_ips': array([ 1.5, 14.5, 18.5, 20.5]),
  'right_ips': array([ 9.5, 17.5, 19.5, 21.5])})

看见 documentation 更多例子。

有多种方法可以解决这个问题。这里列出了一种方法。 您可以创建一个自定义函数,并在查找mimima时使用最大值来处理边缘情况。

import numpy as np
a = np.array([2,2,10,4,4,4,5,6,7,2,6,5,5,7,7,1,1])

def local_min(a):
    temp_list = list(a)
    maxval = max(a) #use max while finding minima
    temp_list = temp_list + [maxval] #handles last value edge case.

    prev = maxval #prev stores last value seen
    loc = 0 #used to store starting index of minima
    count = 0 #use to count repeated values
    #match_start = False
    matches = []
    for i in range(0, len(temp_list)): #need to check all values including the padded value
        if prev == temp_list[i]:
            if count > 0: #only increment for minima candidates
                count += 1
        elif prev > temp_list[i]:
            count = 1
            loc = i
    #        match_start = True
        else: #prev < temp_list[i]
            if count > 0:
                matches.append((loc, count))
            count = 0
            loc = i
        prev = temp_list[i]
    return matches

result = local_min(a)

for match in result:
    print ("{} minima found starting at location {} and ending at location {}".format(
            match[1], 
            match[0],
            match[0] + match[1] -1))

如果这对你有用,请告诉我。这个想法很简单,你想在列表中迭代一次,并在看到它们时保持存储极小值。通过在每一端填充最大值来处理边。(或填充最后一端,并使用最大值进行初始比较)

下面是一个基于将数组限制为一系列窗口的答案:

import numpy as np
from numpy.lib.stride_tricks import as_strided

def windowstride(a, window):
    return as_strided(a, shape=(a.size - window + 1, window), strides=2*a.strides)

def local_min(a, maxwindow=None, doends=True):
    if doends: a = np.pad(a.astype(float), 1, 'constant', constant_values=np.inf)
    if maxwindow is None: maxwindow = a.size - 1

    mins = []
    for i in range(3, maxwindow + 1):
        for j,w in enumerate(windowstride(a, i)):
            if (w[0] > w[1]) and (w[-2] < w[-1]):
                if (w[1:-1]==w[1]).all():
                    mins.append((j, j + i - 2))

    mins.sort()
    return mins

测试它:

test03=np.array([2,2,10,4,4,4,5,6,7,2,6,5,5,7,7,1,1])
local_min(test03)

输出:

[(0, 2), (3, 6), (9, 10), (11, 13), (15, 17)]

不是最有效的算法,但至少它很短。我敢肯定是的 O(n^2) ,因为 1/2*(n^2 + n) 在windows上进行迭代。这只是部分矢量化,因此可能有一种方法可以改进它。

编辑

为了澄清,输出是包含局部最小值运行的切片的索引。他们跑完一段是故意的(有人只是试图在编辑中“修正”)。您可以使用输出在输入数组中的最小值切片上进行迭代,如下所示:

for s in local_mins(test03):
    print(test03[slice(*s)])

输出:

[2 2]
[4 4 4]
[2]
[5 5]
[1 1]

纯numpy解决方案(修订答案):

import numpy as np 

y = np.array([2,2,10,4,4,4,5,6,7,2,6,5,5,7,7,1,1])
x = np.r_[y[0]+1, y, y[-1]+1]   # pad edges, gives possibility for minima

ups,   = np.where(x[:-1] < x[1:])
downs, = np.where(x[:-1] > x[1:])

minend = ups[np.unique(np.searchsorted(ups, downs))]
minbeg = downs[::-1][np.unique(np.searchsorted(-downs[::-1], -ups[::-1]))][::-1]

minlen = minend - minbeg

for line in zip(minlen, minbeg, minend-1): print("set of %d minima %d - %d" % line)

这给

set of 2 minima 0 - 1
set of 3 minima 3 - 5
set of 1 minima 9 - 9
set of 2 minima 11 - 12
set of 2 minima 15 - 16

np.searchsorted(ups, downs) 在每次失败后找到第一个上升点。这是最小值的“真实”结尾。 对于极小值的开始,我们做了类似的事情,但现在是相反的顺序。

它正在为这个例子工作,但尚未完全测试。但我认为这是一个很好的起点。

你可以用 argrelmax ,只要没有多个连续的相等元素,那么首先需要对数组进行游程编码,然后使用argrelmax(或 argrelmin ):

import numpy as np
from scipy.signal import argrelmax
from itertools import groupby


def local_max_scipy(a):
    start = 0
    result = [[a[0] - 1, 0, 0]]  # this is to guarantee the left edge is included
    for k, g in groupby(a):
        length = sum(1 for _ in g)
        result.append([k, start, length])
        start += length
    result.append([a[-1] - 1, 0, 0])  # this is to guarantee the right edge is included

    arr = np.array(result)
    maxima, = argrelmax(arr[:, 0])
    return arr[maxima]


test03 = np.array([2, 2, 10, 4, 4, 4, 5, 6, 7, 2, 6, 5, 5, 7, 7, 1, 1])
output = local_max_scipy(test03)

for val, start, length in output:
    print(f'set of {length} maxima start:{start} end:{start + length}')

输出

set of 1 maxima start:2 end:3
set of 1 maxima start:8 end:9
set of 1 maxima start:10 end:11
set of 2 maxima start:13 end:15