Recherche de maxima / minima locaux avec Numpy dans un tableau numpy 1D

Pouvez-vous suggérer une fonction de module de numpy / scipy qui peut trouver des maxima / minima locaux dans un tableau numpy 1D? De toute évidence, l'approche la plus simple qui soit consiste à jeter un coup d'œil aux voisins les plus proches, mais j'aimerais avoir une solution acceptée qui fasse partie de la distribution numpy.

Si vous recherchez toutes les entrées du tableau 1d aplus petites que leurs voisins, vous pouvez essayer

numpy.r_[True, a[1:] < a[:-1]] & numpy.r_[a[:-1] < a[1:], True]

Vous pouvez également lisser votre tableau avant cette étape en utilisant numpy.convolve().

Je ne pense pas qu'il y ait une fonction dédiée à cela.

Dans SciPy> = 0,11

import numpy as np
from scipy.signal import argrelextrema

x = np.random.random(12)

# for local maxima
argrelextrema(x, np.greater)

# for local minima
argrelextrema(x, np.less)

Produit

>>> x
array([ 0.56660112,  0.76309473,  0.69597908,  0.38260156,  0.24346445,
    0.56021785,  0.24109326,  0.41884061,  0.35461957,  0.54398472,
    0.59572658,  0.92377974])
>>> argrelextrema(x, np.greater)
(array([1, 5, 7]),)
>>> argrelextrema(x, np.less)
(array([4, 6, 8]),)

Notez que ce sont les indices de x qui sont locaux max / min. Pour obtenir les valeurs, essayez:

>>> x[argrelextrema(x, np.greater)[0]]

scipy.signalfournit également argrelmaxet argrelminpour trouver respectivement des maxima et des minima.

Pour les courbes sans trop de bruit, je recommande le petit extrait de code suivant:

from numpy import *

# example data with some peaks:
x = linspace(0,4,1e3)
data = .2*sin(10*x)+ exp(-abs(2-x)**2)

# that's the line, you need:
a = diff(sign(diff(data))).nonzero()[0] + 1 # local min+max
b = (diff(sign(diff(data))) > 0).nonzero()[0] + 1 # local min
c = (diff(sign(diff(data))) < 0).nonzero()[0] + 1 # local max


# graphical output...
from pylab import *
plot(x,data)
plot(x[b], data[b], "o", label="min")
plot(x[c], data[c], "o", label="max")
legend()
show()

Le +1est important, car diffréduit le numéro d'index d'origine.

Une autre approche (plus de mots, moins de code) qui peut aider:

Les emplacements des maxima et minima locaux sont également les emplacements des passages par zéro de la première dérivée. Il est généralement beaucoup plus facile de trouver des passages à zéro que de trouver directement les maxima et minima locaux.

Malheureusement, la première dérivée a tendance à "amplifier" le bruit, de sorte que lorsqu'un bruit significatif est présent dans les données d'origine, la première dérivée n'est mieux utilisée qu'après que les données d'origine ont subi un certain degré de lissage.

Étant donné que le lissage est, dans le sens le plus simple, un filtre passe-bas, le lissage est souvent mieux (enfin, le plus facilement) effectué en utilisant un noyau de convolution, et "façonner" ce noyau peut fournir une quantité surprenante de capacités de préservation / amélioration des fonctionnalités . Le processus de recherche d'un noyau optimal peut être automatisé en utilisant une variété de moyens, mais le meilleur peut être une simple force brute (très rapide pour trouver de petits noyaux). Un bon noyau déformera (comme prévu) massivement les données originales, mais n'affectera PAS l'emplacement des pics / vallées d'intérêt.

Heureusement, assez souvent, un noyau approprié peut être créé via un simple SWAG ("supposé instruit"). La largeur du noyau de lissage doit être un peu plus large que le pic "intéressant" le plus large attendu dans les données d'origine, et sa forme ressemblera à ce pic (une ondelette à échelle unique). Pour les noyaux préservant la moyenne (ce que tout bon filtre de lissage devrait être), la somme des éléments du noyau doit être exactement égale à 1,00, et le noyau doit être symétrique par rapport à son centre (ce qui signifie qu'il aura un nombre impair d'éléments.

Étant donné un noyau de lissage optimal (ou un petit nombre de noyaux optimisés pour différents contenus de données), le degré de lissage devient un facteur d'échelle pour (le «gain») du noyau de convolution.

La détermination du degré "correct" (optimal) de lissage (gain du noyau de convolution) peut même être automatisée: Comparez l'écart type des premières données dérivées avec l'écart type des données lissées. Comment le rapport des deux écarts types change avec les changements du degré de came de lissage être utilisé pour prédire les valeurs de lissage efficaces. Quelques exécutions manuelles de données (qui sont vraiment représentatives) devraient suffire.

Toutes les solutions précédentes publiées ci-dessus calculent la première dérivée, mais elles ne la traitent pas comme une mesure statistique, et les solutions ci-dessus n'essaient pas non plus d'effectuer un lissage préservant / améliorant la fonctionnalité (pour aider les pics subtils à "sauter au-dessus" du bruit).

Enfin, la mauvaise nouvelle: trouver de "vrais" pics devient une douleur royale lorsque le bruit a également des caractéristiques qui ressemblent à de vrais pics (bande passante qui se chevauchent). La solution suivante plus complexe consiste généralement à utiliser un noyau de convolution plus long (une «ouverture de noyau plus large») qui prend en compte la relation entre les pics «réels» adjacents (tels que les taux minimum ou maximum d'occurrence de pic), ou d'utiliser plusieurs la convolution passe en utilisant des noyaux de largeurs différentes (mais seulement si elle est plus rapide: c'est une vérité mathématique fondamentale que les convolutions linéaires exécutées en séquence peuvent toujours être convolutionnées ensemble en une seule convolution). Mais il est souvent beaucoup plus facile de trouver d'abord une séquence de noyaux utiles (de largeurs variables) et de les convoluer ensemble que de trouver directement le noyau final en une seule étape.

Espérons que cela fournisse suffisamment d'informations pour permettre à Google (et peut-être un bon texte de statistiques) de combler les lacunes. J'aurais vraiment aimé avoir le temps de fournir un exemple travaillé, ou un lien vers un. Si quelqu'un en trouve un en ligne, veuillez le poster ici!

Depuis la version 1.1 de SciPy, vous pouvez également utiliser find_peaks . Voici deux exemples tirés de la documentation elle-même.

En utilisant l' heightargument, on peut sélectionner tous les maxima au-dessus d'un certain seuil (dans cet exemple, tous les maxima non négatifs; cela peut être très utile si l'on doit traiter une ligne de base bruyante; si vous voulez trouver des minima, multipliez simplement votre entrée par -1):

import matplotlib.pyplot as plt
from scipy.misc import electrocardiogram
from scipy.signal import find_peaks
import numpy as np

x = electrocardiogram()[2000:4000]
peaks, _ = find_peaks(x, height=0)
plt.plot(x)
plt.plot(peaks, x[peaks], "x")
plt.plot(np.zeros_like(x), "--", color="gray")
plt.show()

Un autre argument extrêmement utile est distance, qui définit la distance minimale entre deux pics:

peaks, _ = find_peaks(x, distance=150)
# difference between peaks is >= 150
print(np.diff(peaks))
# prints [186 180 177 171 177 169 167 164 158 162 172]

plt.plot(x)
plt.plot(peaks, x[peaks], "x")
plt.show()

Pourquoi ne pas utiliser la fonction intégrée signal.find_peaks_cwt de Scipy pour faire le travail?

from scipy import signal
import numpy as np

#generate junk data (numpy 1D arr)
xs = np.arange(0, np.pi, 0.05)
data = np.sin(xs)

# maxima : use builtin function to find (max) peaks
max_peakind = signal.find_peaks_cwt(data, np.arange(1,10))

# inverse  (in order to find minima)
inv_data = 1/data
# minima : use builtin function fo find (min) peaks (use inversed data)
min_peakind = signal.find_peaks_cwt(inv_data, np.arange(1,10))

#show results
print "maxima",  data[max_peakind]
print "minima",  data[min_peakind]

résultats:

maxima [ 0.9995736]
minima [ 0.09146464]

Cordialement

Mise à jour: je n'étais pas satisfait du dégradé, donc je l'ai trouvé plus fiable à utiliser numpy.diff. S'il vous plaît laissez-moi savoir s'il fait ce que vous voulez.

En ce qui concerne la question du bruit, le problème mathématique est de localiser les maxima / minima si nous voulons regarder le bruit, nous pouvons utiliser quelque chose comme convolve qui a été mentionné précédemment.

import numpy as np
from matplotlib import pyplot

a=np.array([10.3,2,0.9,4,5,6,7,34,2,5,25,3,-26,-20,-29],dtype=np.float)

gradients=np.diff(a)
print gradients


maxima_num=0
minima_num=0
max_locations=[]
min_locations=[]
count=0
for i in gradients[:-1]:
        count+=1

    if ((cmp(i,0)>0) & (cmp(gradients[count],0)<0) & (i != gradients[count])):
        maxima_num+=1
        max_locations.append(count)     

    if ((cmp(i,0)<0) & (cmp(gradients[count],0)>0) & (i != gradients[count])):
        minima_num+=1
        min_locations.append(count)


turning_points = {'maxima_number':maxima_num,'minima_number':minima_num,'maxima_locations':max_locations,'minima_locations':min_locations}  

print turning_points

pyplot.plot(a)
pyplot.show()

Alors que cette question est vraiment ancienne. Je pense qu'il existe une approche beaucoup plus simple dans numpy (une ligne unique).

import numpy as np

list = [1,3,9,5,2,5,6,9,7]

np.diff(np.sign(np.diff(list))) #the one liner

#output
array([ 0, -2,  0,  2,  0,  0, -2])

Pour trouver un max ou un min local, nous voulons essentiellement savoir quand la différence entre les valeurs de la liste (3-1, 9-3 ...) passe de positive à négative (max) ou négative à positive (min). Par conséquent, nous trouvons d'abord la différence. Ensuite, nous trouvons le signe, puis nous trouvons les changements de signe en reprenant la différence. (Un peu comme une première et une deuxième dérivée dans le calcul, seulement nous avons des données discrètes et n'avons pas de fonction continue.)

La sortie dans mon exemple ne contient pas les extrema (les première et dernière valeurs de la liste). De plus, tout comme le calcul, si la deuxième dérivée est négative, vous avez max, et si elle est positive, vous avez un min.

Ainsi, nous avons le match suivant:

[1,  3,  9,  5,  2,  5,  6,  9,  7]
    [0, -2,  0,  2,  0,  0, -2]
        Max     Min         Max

Aucune de ces solutions n'a fonctionné pour moi car je voulais également trouver des pics au centre des valeurs répétitives. par exemple, dans

ar = np.array([0,1,2,2,2,1,3,3,3,2,5,0])

la réponse devrait être

array([ 3,  7, 10], dtype=int64)

Je l'ai fait en utilisant une boucle. Je sais que ce n'est pas super propre, mais ça fait le travail.

def findLocalMaxima(ar):
# find local maxima of array, including centers of repeating elements    
maxInd = np.zeros_like(ar)
peakVar = -np.inf
i = -1
while i < len(ar)-1:
#for i in range(len(ar)):
    i += 1
    if peakVar < ar[i]:
        peakVar = ar[i]
        for j in range(i,len(ar)):
            if peakVar < ar[j]:
                break
            elif peakVar == ar[j]:
                continue
            elif peakVar > ar[j]:
                peakInd = i + np.floor(abs(i-j)/2)
                maxInd[peakInd.astype(int)] = 1
                i = j
                break
    peakVar = ar[i]
maxInd = np.where(maxInd)[0]
return maxInd 

import numpy as np
x=np.array([6,3,5,2,1,4,9,7,8])
y=np.array([2,1,3,5,3,9,8,10,7])
sortId=np.argsort(x)
x=x[sortId]
y=y[sortId]
minm = np.array([])
maxm = np.array([])
i = 0
while i < length-1:
    if i < length - 1:
        while i < length-1 and y[i+1] >= y[i]:
            i+=1

        if i != 0 and i < length-1:
            maxm = np.append(maxm,i)

        i+=1

    if i < length - 1:
        while i < length-1 and y[i+1] <= y[i]:
            i+=1

        if i < length-1:
            minm = np.append(minm,i)
        i+=1


print minm
print maxm

minmet maxmcontiennent des indices de minima et de maxima, respectivement. Pour un énorme ensemble de données, cela donnera beaucoup de maximas / minimas, donc dans ce cas, lissez d'abord la courbe, puis appliquez cet algorithme.

Une autre solution utilisant essentiellement un opérateur dilaté:

import numpy as np
from scipy.ndimage import rank_filter

def find_local_maxima(x):
   x_dilate = rank_filter(x, -1, size=3)
   return x_dilate == x

et pour les minima:

def find_local_minima(x):
   x_erode = rank_filter(x, -0, size=3)
   return x_erode == x

De plus, scipy.ndimagevous pouvez remplacer rank_filter(x, -1, size=3)par grey_dilationet rank_filter(x, 0, size=3)par grey_erosion. Cela ne nécessitera pas de tri local, donc c'est légèrement plus rapide.

Un autre:

def local_maxima_mask(vec):
    """
    Get a mask of all points in vec which are local maxima
    :param vec: A real-valued vector
    :return: A boolean mask of the same size where True elements correspond to maxima. 
    """
    mask = np.zeros(vec.shape, dtype=np.bool)
    greater_than_the_last = np.diff(vec)>0  # N-1
    mask[1:] = greater_than_the_last
    mask[:-1] &= ~greater_than_the_last
    return mask