在1D-NumPy阵列中查找单个/一组局部最大值/最小值(再次)

Question

我想有一个函数可以检测数组中局部最大值/最小值的位置(即使有一组局部最大值/最小值).例:

鉴于阵列

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

我想有一个输出像:

set of 2 local minima => array[0]:array[1]
set of 3 local minima => array[3]:array[5]
local minima, i = 9
set of 2 local minima => array[11]:array[12]
set of 2 local minima => array[15]:array[16]

从示例中可以看出,不仅检测到奇异值,还检测局部最大值/最小值.

我知道在这个问题中有很多好的答案和想法,但是他们都没有完成所描述的工作:他们中的一些人只是忽略了数组的极端点而忽略了局部最小值/最大值的集合.

在提出这个问题之前,我自己编写了一个函数,它完全按照上面的描述进行(函数在这个问题的最后:local_min(a).通过我做的测试,它可以正常工作).

问题:但是,我也确信这不是使用Python的最佳方式.我可以使用内置函数,API,库等吗？还有其他功能建议吗？一行指令？一个完整的矢量解决方案？

def local_min(a):
    candidate_min=0
    for i in range(len(a)):

        # Controlling the first left element
        if i==0 and len(a)>=1:
            # If the first element is a singular local minima
            if a[0]<a[1]:
                print("local minima, i = 0")
            # If the element is a candidate to be part of a set of local minima
            elif a[0]==a[1]:
                candidate_min=1
        # Controlling the last right element
        if i == (len(a)-1) and len(a)>=1:
            if candidate_min > 0:
                if a[len(a)-1]==a[len(a)-2]:
                    print("set of " + str(candidate_min+1)+ " local minima => array["+str(i-candidate_min)+"]:array["+str(i)+"]")
            if a[len(a)-1]<a[len(a)-2]:
                print("local minima, i = " + str(len(a)-1))
        # Controlling the other values in the middle of the array
        if i>0 and i<len(a)-1 and len(a)>2:
            # If a singular local minima
            if (a[i]<a[i-1] and a[i]<a[i+1]):
                print("local minima, i = " + str(i))
                # print(str(a[i-1])+" > " + str(a[i]) + " < "+str(a[i+1])) #debug
            # If it was found a set of candidate local minima
            if candidate_min >0:
                # The candidate set IS a set of local minima
                if a[i] < a[i+1]:
                    print("set of " + str(candidate_min+1)+ " local minima => array["+str(i-candidate_min)+"]:array["+str(i)+"]")
                    candidate_min = 0
                # The candidate set IS NOT a set of local minima
                elif a[i] > a[i+1]:
                    candidate_min = 0
                # The set of local minima is growing
                elif a[i] == a[i+1]:
                    candidate_min = candidate_min + 1
                # It never should arrive in the last else
                else:
                    print("Something strange happen")
                    return -1
            # If there is a set of candidate local minima (first value found)
            if (a[i]<a[i-1] and a[i]==a[i+1]):
                candidate_min = candidate_min + 1

注意:我尝试用一​​些注释来丰富代码,以了解我想要做什么.我知道我建议的功能不干净,只打印可以存储并在最后返回的结果.它是为了举个例子而写的.我建议的算法应该是O(n).

更新:

有人建议导入from scipy.signal import argrelextrema并使用如下功能:

def local_min_scipy(a):
    minima = argrelextrema(a, np.less_equal)[0]
    return minima

def local_max_scipy(a):
    minima = argrelextrema(a, np.greater_equal)[0]
    return minima

拥有这样的东西是我真正想要的.但是,当局部最小值/最大值集合具有两个以上的值时,它无法正常工作.例如:

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

print(local_max_scipy(test03))

输出是:

[ 0  2  4  8 10 13 14 16]

当然,test03[4]我有一个最小而不是最大.我该如何解决这个问题？(我不知道这是否是另一个问题,或者这是否是正确的地方.)

Answer 1

完整的矢量解决方案:

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补丁会这样丑陋.

为了解决这个问题,我提出了另外两个解决方案,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加速兼容的直接解决方案:

#@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)  

Answer 2

我认为另一个函数 fromscipy.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])})

有关更多示例，请参阅文档。