在 1D-NumPy 数组中查找奇异值/局部最大值/最小值集（再次）

Posted 2023-02-23

【中文标题】在 1D-NumPy 数组中查找奇异值/局部最大值/最小值集（再次）

【英文标题】：Finding singulars/sets of local maxima/minima in a 1D-NumPy array (once again)

【发布时间】：2019-04-27 05:39:38

【问题描述】：

我希望有一个函数可以检测数组中的局部最大值/最小值（即使有一组局部最大值/最小值）。示例：

给定数组

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]

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

我知道在this question 中有很多好的答案和想法，但没有一个能够完成所描述的工作：其中一些简单地忽略了数组的极值点，并且都忽略了局部最小值/最大值的集合。

在问这个问题之前，我自己写了一个函数，它完全符合我上面描述的功能（该函数在这个问题的末尾：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

注意：我尝试用一​​些 cmets 来丰富代码，让大家了解我的工作。我知道我建议的功能是 不干净，只打印可以存储和返回的结果 在末尾。它是为了举例而写的。我提出的算法应该是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] 我有一个最小值而不是最大值。如何解决此行为？ （我不知道这是否是另一个问题，或者这是否是提出问题的正确地方。）

【问题讨论】：

有趣的问题，快速搜索似乎表明没有预先构建的解决方案。但是，为此设计一个简约的解决方案应该足够简单。我可以想到两种方法。让我尝试实现一个，看看它是否像我认为的那样干净。

你想如何处理egdes？

【参考方案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) 在每次下跌之后找到第一个上涨。这是最小值的“真正”结束。 对于最小值的开始，我们以类似的方式进行操作，但现在以相反的顺序进行。

它适用于示例，但尚未经过全面测试。但我想说一个很好的起点。

【讨论】：

为什么在您的原始答案中使用消除 print？新版本的 Python？我正在使用 python 3.6，它与您的原始解决方案不打印任何内容

BUG：使用y = np.array([2,2,10,4,4,4,5,6,7,2,6,5,5,7,7,7,7,7,7,7,1,1]) 它还检测到最大值！所以解决方案不能这么紧凑

如果没有print，它会显示在交互式python2/3中。因此，我将print 视为代码噪声。在 python3 中，我发现输入括号仍然很烦人。

我已经完全修改了我的答案，解决了这个错误。

【参考方案2】：

你可以使用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

【参考方案3】：

我认为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。

【讨论】：

输出是什么意思？

发现局部峰的指数

我怎么知道它是只有一个值还是一组相等的值？有办法吗？还是我应该手动添加其他代码？

您是否在一个峰值中有 2 个以上的值时测试了解决方案？

test04 = mp.array([1,1,5,5,5,5,5,5,5,5,1,1,1,1,1,5,5,5,1,5,1,5,1])。对不起，如果我不能自己测试，我现在不能

【参考方案4】：

一个完整的矢量解决方案：

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 补丁会很丑。

为了解决这个问题，我提出了另外 2 个解决方案，先使用 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)  

【讨论】：

这很漂亮。感谢您添加 cmets，否则很难遵循。您能否详细说明创建拆分位置的最后第二行？

最后一行 list(zip(begin[:k],end[:k])) 占用了 Numba 解决方案总运行时间的大约 80-90%。返回一个简单的 numpy 数组会更快，例如。 out=np.empty((k,2),dtype=a.dtype)for i in range(k):out[i,0]=begin[i]out[i,1]=end[i]return out

【参考方案5】：

这是一个基于将数组重新跨入可迭代窗口的答案：

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) 窗口可以迭代。这只是部分矢量化的，因此可能有办法改进它。

编辑

为了澄清，输出是包含局部最小值运行的切片的索引。他们超过运行结束的事实是故意的（有人只是试图在编辑中“修复”它）。您可以使用输出来迭代输入数组中的最小值切片，如下所示：

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

输出：

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

【讨论】：

+1 它可以正常工作，但会成倍增加复杂性。然而，使用 Windows 的技巧是我不知道的，这就是为什么 +1 ！我仍然可以将想法用于其他上下文

@Leos313 是的，strides 可用于向量化各种迭代。实际上，我昨天才学会如何使用它们。事实证明，重排并不是解决这个问题的最佳选择，但是当锤子又亮又新时，“当你有一把锤子时，每个问题看起来都像钉子”。

【参考方案6】：

可以有多种方法来解决这个问题。此处列出的一种方法。 您可以创建一个自定义函数，并在查找 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))

如果这对你有用，请告诉我。这个想法很简单，您希望遍历列表一次并在看到它们时继续存储最小值。通过在两端填充最大值来处理边缘。 （或通过填充最后一端，并使用最大值进行初始比较）

【讨论】：

用test03 测试过（你可以在问题中找到）并且有些奇怪。例如输出的第 3 行是 0 minima found starting at location 6 and ending at location 5

+1 它可以正常工作，而且可以肯定的是，它比我提出的解决方案更好。它仍然可以与for 一起使用，并且 Python 可以与数组和矩阵一起使用（例如 Matlab）。但是，使用列表是个好主意！

当边缘不是最小时，我检测到一个错误。给我一些时间来纠正。

哎呀，不错的选择。这是同一件事，需要确保只为匹配和“增量”添加大于 0 的计数。 @BM