2013-07-28 4 views
3

Я хотел бы получить наиболее эффективный способ найти локальные максимумы в огромный набор точек данных, содержащий тысячи значений. В качестве входных данных используются два длинных списка с значениями x и y.Поиск локальных максимумов xy графика точек данных с numpy?

Рассмотрим простой пример:

xval = [-0.15, -0.02, 0.1, 0.22, 0.36, 0.43, 0.58, 0.67, 0.79, 0.86, 0.96 ] 

yval = [-0.09, 0.13, -0.01, -0.1, -0.05, 0.2, 0.56, 0.47, 0.35, 0.43, 0.69] 

enter image description here

Желаемая выход список с индексами точек пика, то здесь locMaxId = [1,6,10]. Сравнение ближайших соседей - это решение, но для значений 10 тыс.?

ответ

4

Вы можете позволить NumPy обрабатывать итерацию, т.е. векторизовать его:

def local_maxima(xval, yval): 
    xval = np.asarray(xval) 
    yval = np.asarray(yval) 

    sort_idx = np.argsort(xval) 
    yval = yval[sort_idx] 
    gradient = np.diff(yval) 
    maxima = np.diff((gradient > 0).view(np.int8)) 
    return np.concatenate((([0],) if gradient[0] < 0 else()) + 
          (np.where(maxima == -1)[0] + 1,) + 
          (([len(yval)-1],) if gradient[-1] > 0 else())) 

EDIT Таким образом, код первого вычисляет отклонение от каждой точки до NEX (gradient). Следующий шаг немного сложный ... Если вы делаете np.diff((gradient > 0), результирующий булевой массив равен True, где есть изменение от роста (> 0) до не растущего (<= 0). Сделав его подписанным int того же размера, что и логический массив, вы можете отличать от переходов от роста до не растущего (-1) к противоположному (+1). Принимая .view(np.int8) знакового целочисленного типа с тем же размером dtype, что и в булевом массиве, мы избегаем копирования данных, как это было бы, если бы мы делали менее хакерский .astype(int). Все, что осталось, - это забота о первой и последней точках и объединение всех точек в один массив. Одна вещь, которую я узнал сегодня, состоит в том, что если вы включаете пустой список в кортеж, который вы отправляете в np.concatenate, он выдается как пустой массив dtype np.float, и это в конечном итоге является dtype результата, следовательно, более сложная конкатенация пустые кортежи в приведенном выше коде.

Он работает:

In [2]: local_maxima(xval, yval) 
Out[2]: array([ 1, 6, 10], dtype=int64) 

И достаточно быстро:

In [3]: xval = np.random.rand(10000) 

In [4]: yval = np.random.rand(10000) 

In [5]: local_maxima(xval, yval) 
Out[5]: array([ 0, 2, 4, ..., 9991, 9995, 9998], dtype=int64) 

In [6]: %timeit local_maxima(xval, yval) 
1000 loops, best of 3: 1.16 ms per loop 

Кроме того, большую часть времени является преобразование данных из списков массивов и их сортировки. Если ваши данные уже отсортированы и хранятся в массивах, вы, вероятно, можете повысить производительность по сравнению с выше в 5 раз.

+0

Да, именно это я и искал! Возможно, вы могли бы добавить комментарии для неопытных пользователей :) – ptaeck

+0

@ptaeck Посмотрите, имеет ли смысл какой-либо смысл, я часто объясняю себя лучше на Python, чем на английском ... – Jaime

+0

Да, теперь это намного лучше! Если вы хотите, [этот] (http://bit.ly/11q0aSN) немного сложнее. – ptaeck