如何在Python中将此峰值查找for循环向量化？

Question

如何在Python中将此峰值查找for循环向量化？

pythonnumpyscipysignal-processingvectorization

6

基本上我正在编写一个峰值查找函数，需要能够在基准测试中击败scipy.argrelextrema。这是我使用的数据和代码链接：

https://drive.google.com/open?id=1U-_xQRWPoyUXhQUhFgnM3ByGw-1VImKB

如果此链接已过期，您可以在Dukascopy银行的在线历史数据下载器中找到数据。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

data = pd.read_csv('EUR_USD.csv')
data.columns = ['Date', 'open', 'high', 'low', 'close','volume']

data.Date = pd.to_datetime(data.Date, format='%d.%m.%Y %H:%M:%S.%f')

data = data.set_index(data.Date)

data = data[['open', 'high', 'low', 'close']]

data = data.drop_duplicates(keep=False)

price = data.close.values

def fft_detect(price, p=0.4):

    trans = np.fft.rfft(price)
    trans[round(p*len(trans)):] = 0
    inv = np.fft.irfft(trans)
    dy = np.gradient(inv)
    peaks_idx = np.where(np.diff(np.sign(dy)) == -2)[0] + 1
    valleys_idx = np.where(np.diff(np.sign(dy)) == 2)[0] + 1

    patt_idx = list(peaks_idx) + list(valleys_idx)
    patt_idx.sort()

    label = [x for x in np.diff(np.sign(dy)) if x != 0]

    # Look for Better Peaks

    l = 2

    new_inds = []

    for i in range(0,len(patt_idx[:-1])):

        search = np.arange(patt_idx[i]-(l+1),patt_idx[i]+(l+1))

        if label[i] == -2:
            idx = price[search].argmax()
        elif label[i] == 2:
            idx = price[search].argmin()

        new_max = search[idx]
        new_inds.append(new_max)

    plt.plot(price)
    plt.plot(inv)
    plt.scatter(patt_idx,price[patt_idx])
    plt.scatter(new_inds,price[new_inds],c='g')
    plt.show()

    return peaks_idx, price[peaks_idx]

它基本上使用快速傅里叶变换（FFT）平滑数据，然后取导数以找到平滑数据的最小和最大索引，然后在未平滑的数据中找到相应的峰值。由于某些平滑效果，有时它找到的峰值不理想，因此我运行此for循环来搜索由指定的边界之间每个索引的更高或更低点。我需要帮助将此for循环向量化！我不知道如何做。没有for循环，我的代码比scipy.argrelextrema快约50％，但是for循环会减慢速度。因此，如果我能找到一种将其向量化的方法，它将是一个非常快速且非常有效的scipy.argrelextrema替代方法。这两个图像分别表示没有和有for循环的数据。

- ddm-j

1

嗨，发一些数据怎么样？那样试验起来更容易。 - mortysporty

@mortysporty 在我的编辑中添加了数据导入和链接。 - ddm-j

太棒了 :) 我怀疑我能解决它，但我会试一试！ - mortysporty

2个回答

2

这里有一个替代方案... 它使用列表推导式，通常比for循环更快

l = 2

# Define the bounds beforehand, its marginally faster than doing it in the loop
upper = np.array(patt_idx) + l + 1
lower = np.array(patt_idx) - l - 1

# List comprehension...
new_inds = [price[low:hi].argmax() + low if lab == -2 else 
            price[low:hi].argmin() + low 
            for low, hi, lab in zip(lower, upper, label)]

# Find maximum within each interval
new_max = price[new_inds]
new_global_max = np.max(new_max)

- mortysporty

谢谢，我已经对你的方法进行了基准测试，它比我的快大约30%。结果：第二种方法：0.00032 28.42% 加速比。 - ddm-j

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- kevinkayaks · Accepted Answer

这可能就是你需要的。虽然不完美，但希望能得到你想要的并向你展示如何进行向量化操作。欢迎听取你的任何改进意见。

label = np.array(label[:-1]) # not sure why this is 1 unit longer than search.shape[0]? 

# the idea is to make the index matrix you're for looping over row by row all in one go. 
# This part is sloppy and you can improve this generation. 

search = np.vstack((np.arange(patt_idx[i]-(l+1),patt_idx[i]+(l+1)) for i in range(0,len(patt_idx[:-1])))) # you can refine this. 

# then you can make the price matrix

price = price[search]

# and you can swap the sign of elements so you only need to do argmin instead of both argmin and argmax 

price[label==-2] = - price[label==-2]

# now find the indices of the minimum price on each row 

idx = np.argmin(price,axis=1)

# and then extract the refined indices from the search matrix 

new_inds = search[np.arange(idx.shape[0]),idx] # this too can be cleaner. 
# not sure what's going on here so that search[:,idx] doesn't work for me
# probably just a misunderstanding

我发现这个操作可以复制您的结果，但我没有计时。我怀疑搜索生成速度相当慢，但可能仍然比您的 for 循环要快。

编辑：

以下是更好的生成“search”的方法：

patt_idx = np.array(patt_idx)
starts = patt_idx[:-1]-(l+1)
stops = patt_idx[:-1]+(l+1)
ds = stops-starts
s0 = stops.shape[0]
s1 = ds[0]
search = np.reshape(np.repeat(stops - ds.cumsum(), ds) + np.arange(ds.sum()),(s0,s1))