Finding singulars/sets of local maxima/minima in a 1D-NumPy array (once again)

Question

I would like to have a function that can detect where the local maxima/minima are in an array (even if there is a set of local maxima/minima). Example:

Given the array

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

I would like to have an output like:

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]

As you can see from the example, not only are the singular values detected but, also, sets of local maxima/minima.

I know in this question there are a lot of good answers and ideas, but none of them do the job described: some of them simply ignore the extreme points of the array and all ignore the sets of local minima/maxima.

Before asking the question, I wrote a function by myself that does exactly what I described above (the function is at the end of this question: local_min(a). With the test I did, it works properly).

Question: However, I am also sure that is NOT the best way to work with Python. Are there builtin functions, APIs, libraries, etc. that I can use? Any other function suggestion? A one-line instruction? A full vectored solution?

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

Note: I tried to enrich the code with some comments to let understand what I do. I know that the function that I propose is not clean and just prints the results that can be stored and returned at the end. It was written to give an example. The algorithm I propose should be O(n).

UPDATE:

Somebody was suggesting to import from scipy.signal import argrelextrema and use the function like:

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

To have something like that is what I am really looking for. However, it doesn't work properly when the sets of local minima/maxima have more than two values. For example:

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

The output is:

[ 0  2  4  8 10 13 14 16]

Of course in test03[4] I have a minimum and not a maximum. How do I fix this behavior? (I don't know if this is another question or if this is the right place where to ask it.)

Answer 1

A full vectored solution:

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]

EDIT

Unfortunately, This detects also maxima as soon as they are at least 3 items large, since they are seen as flat local minima. A numpy patch will be ugly this way.

To solve this problem I propose 2 other solutions, with numpy, then with numba.

Whith numpy using 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)    

A direct solution compatible with 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)  

Answer 2

I think another function from scipy.signal would be interesting.

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 has lots of options and might be quite useful, especially for noisy signals.

Update

The function is really powerful and versatile. You can set several parameters for peak minimal width, height, distance from each other and so on. As example:

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])})

See documentation for more examples.

Answer 3

There can be multiple ways to solve this. One approach listed here. You can create a custom function, and use the maximums to handle edge cases while finding 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))

Let me know if this does the trick for you. The idea is simple, you want to iterate through the list once and keep storing minima as you see them. Handle the edges by padding with maximum values on either end. (or by padding the last end, and using the max value for initial comparison)

Answer 4

Here's an answer based on restriding the array into an iterable of windows:

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

Testing it out:

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

Output:

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

Not the most efficient algorithm, but at least it's short. I'm pretty sure it's O(n^2), since there's roughly 1/2*(n^2 + n) windows to iterate over. This is only partially vectorized, so there may be a way to improve it.

Edit

To clarify, the output is the indices of the slices that contain the runs of local minimum values. The fact that they go one past the end of the run is intentional (someone just tried to "fix" that in an edit). You can use the output to iterate over the slices of minimum values in your input array like this:

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

Output:

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