CODEWITHSUNDEEP
pythonarraysalgorithmnumpy
Taenyfan
Taenyfan

Reputation: 13

Improving algorithm to find local maxima of 1D numpy array

I am trying to find the local maxima of the function f(x) = (sin(x)/x)^2. For approximate solutions, I initialised two variables x and y and first plotted a graph to have a visual representation. 

x = np.linspace(-20.0, 20.0, num=int((20+20)/0.01)) 
y = np.power(np.sin(x)/x, 2)
plt.plot(x, y, 'r.', markersize= 1)
plt.show()

This shows graph.

I then tried to create an algorithm to find the Maxima:

def returnMaxima(num, x, y):
    """
    number, np.array, np.array -> list
    num: number of maxima needed | x: x 1D array | y: y 1D array
    returns [[x1,y1], [x2,y2]...] in descending order of y
    """
    allMaximaPoints = [] # stores all Maxima points
    reqMaximaPoints = [] # stores num Maxima points
    for i in range(y.size): 
        # for first y value
        if i == 0: 
            if y[i] > y[i+1]:
                allMaximaPoints += [[x[i], y[i]], ]
        # for last y value
        elif i == y.size - 1:
            if y[i] > y[i-1]:
                allMaximaPoints += [[x[i], y[i]], ]
        # for rest y values
        else: 
            if y[i] > y[i-1] and y[i] > y[i+1]:
                allMaximaPoints += [[x[i], y[i]], ]
    # extract largest maximas from allMaximaPoints
    while num > 0: 
        reqMaximaPoints += [max(allMaximaPoints, key=lambda item:item[1]),]
        del allMaximaPoints[allMaximaPoints.index(max(allMaximaPoints, key=lambda item:item[1]))]
        num -= 1
    return reqMaximaPoints

When I tried returnMaxima(2, x, y) I get [[-4.4961240310077528, 0.04719010162459622], [4.4961240310077528, 0.04719010162459622]].

This is incorrect as it skipped the local maxima at x = 0. I suspect it is because the y[i-1] and y[i+1] values adjacent to the maxima at y[i] at x=0 is approximately equal to y[i] causing the code 

else:
    if y[i] > y[i-1] and y[i] > y[i+1]:
        allMaximaPoints += [[x[i], y[i]], ]

to not account for that point. This is because when I changedx = np.linspace(-20.0, 20.0, num=int((20+20)/0.01)) to say x = np.linspace(-20.0, 20.0, num=int((20+20)/0.1)) i.e. larger steps in x, the local maxima at x=0 is correctly found. However, even if I changed the > signs in the above code to >=, that maxima at x=0 is still not counted.

Why is that so? How should I improve my code to get the correct results? Thanks!

Upvotes: 0

Views: 205

Answers (2)

Alok Singhal
Alok Singhal

Reputation: 96101

You might be better off using something like scipy.signal.find_peaks_cwt. Something like:

indices = scipy.signal.find_peaks_cwt(y, [1, 2, 3, 4], noise_perc=50)
plt.plot(x, y)
plt.plot(x[indices], y[indices], 'r.')
plt.show()

Results in:

Peaks

Upvotes: 1

dtauxe
dtauxe

Reputation: 153

Short answer: x is never 0 (and you probably don't want it to be either).

We can test this by running

In [53]: 0 in x
Out[53]: False

This is actually fortunate as x=0 would cause a runtime warning and produce a "nan" (Not a Number) value due to the division by 0. The problem is that the function is undefined at 0.

Otherwise, you are essentially correct about why it doesn't work. If you swap the > for >= you'll see that the output is different:

In [55]: returnMaxima(2, x, y)
Out[55]:
[[-0.0050012503125778096, 0.99999166252624228],
[0.0050012503125778096, 0.99999166252624228]]

This is actually the "correct" output, as these are the largest values and are at the values of x closest to 0.

Upvotes: 0

Related Questions