将Matlab转换为Python-加速循环 [英] Translating Matlab to Python - Speeding up a loop
问题描述
我一直在将一些代码从Matlab转换为Python，我们将其用于分析实验室中的数据.我们有两个时间戳列表，我们想用一个来预示另一个时间戳:对于第一个列表中的每个元素，我们在第二个列表中寻找时间间隔精确的时间戳.如果有的话，我们将它们放在单独的列表中.
I have been translating some code from Matlab to Python that we use to analyse data in our lab. We have two lists of time stamps and we want to use one to herald the other: for every element in the first list we look for time stamps in the second list that have a precise separation in time. In case there are, we place these in a separate list.
这是我正在使用的带有随机数据的Matlab代码类型的可运行示例.因为我不太熟悉Matlab，所以它可能非常粗糙.在下面的 Ctrigger 是触发列表，而 Csignal 是我们要预示的信号列表.对于 Ctrigger 的每个元素，我们查看 Csignal 中是否有元素位于以 offset 为中心的窗口内，且宽度为 gate .所选事件将放置在 Hsignal 中.
Here is an runnable example of the kind of Matlab code I am using, with random data. It is probably VERY crude, as I am not well versed in Matlab. In the following Ctrigger is the trigger list, and Csignal is the signal list that we want to herald. For every element of Ctrigger we look if there are elements in Csignal that are within a window centred on offset, and with width gate. The selected events will be placed in Hsignal.
% Matlab code
Ctrigger = linspace(0, 3000000, (3000000-1)/3);
length_t = length(Ctrigger);
Bsignal = linspace(0, 3000000, (3000000-1)/10);
length_s = length(Bsignal);
noise = reshape(20*rand(length_s,1)-10,[1,length_s]);
Csignal = Bsignal + noise;
offset = 3;
gate = 1;
Hsignal=zeros(length_s,1);
marker = 1;
tic
for j=1:length_t-1
m = marker;
tstart=Ctrigger(j)+offset-gate/2;
tstop=Ctrigger(j)+offset+gate/2;
while(m <= length_s-1)
if(Csignal(m)<tstart)
marker=m;
m=m+1;
end
if(Csignal(m)>=tstart && Csignal(m)<=tstop)
Hsignal(m)=Csignal(m);
m = m+1;
end
if(Csignal(m)>tstop)
break;
end
end
end
toc
Hsignal=Hsignal(Hsignal~=0);
Hsignal = unique(Hsignal);
大约有90'000个事件被选择放置在 Hsignal 中，而Matlab大约需要0.05秒来运行它.我已经介绍了 marker 计数器，因为两个列表 Csignal 和 Ctrigger 区域已经按时排序了. marker 设置在一个先驱窗口的开始处:当我移至下一个触发器时，我不会在所有 Csignal 中再次查看，而只会在该窗口的开始处进行查看.为了避免重复计算，我在最后删除了重复项.
Roughly 90'000 events are selected to be placed in Hsignal, and Matlab takes about 0.05 seconds to run this. I have introduced the marker counter because the two lists Csignal and Ctrigger area already ordered in time. marker is set at the start of one heralding window: when I move to the next trigger I will not look again in all of Csignal, but only from the start of that window. To avoid a double count, I remove the duplicates at the end.
如果您想对代码有所了解，下面是输入和输出的简化版本:
If you want to have an idea of the code, here is a simplified version of the input and output:
Ctrigger = [1, 10, 11, 20, 30, 40, 50, 60]
Csignal = [4, 11, 13, 17, 25, 34, 41, 42, 50, 57, 65]
print(Hsignal)
# [4, 11, 13, 41, 42]
现在，我已经从Matlab复制了此代码，只是对其进行了少许调整以适合python.遵循一些建议，我首先声明包含主算法的函数，然后调用它:
Now, I have copied this code from Matlab, just slightly adjusting it to fit into python. Following some advice I first declare the function that contains the main algorithm, and then call it:
# Python code
def main(list1, list2, list3, delay, window):
marker = 1
for j in range(len(list1)):
m = marker
t_star = list1[j] + delay - window/2
t_sto = list1[j] + delay + window/2
while m < len(list2):
if (list2[m] < t_star):
marker = m
m = m + 1
elif (list2[m] >= t_star and list2[m] <= t_sto):
list3[m] = list2[m]
m = m + 1
elif (list2[m] > t_sto):
break
Ctrigger = range(0, 3000000, 3)
length_t = len(Ctrigger)
Bsignal = range(0, 3000000, 10)
length_s = len(Bsignal)
noise = 1e-05*np.asarray(random.sample(range(-1000000,1000000), int(length_s)))
Csignal = list(np.sort(np.asarray(Bsignal) + noise))
offset = 3
gate = 1
length_t = len(Ctrigger)
length_s = len(Csignal)
Hsignal = list(np.zeros(len(Ctrigger)))
start = time.time()
main(Ctrigger, Csignal, Hsignal, offset, gate)
end = time.time()
Hsignal = np.sort(np.asarray(list(set(Hsignal))))
print(end-start)
类似地，在 Hsignal 中放置了大约90'000个元素.关键问题是python大约需要1.1秒来运行它！我什至尝试了这种替代方法，它消除了一些循环(这里我仍然使用数组，因为我必须将元素添加到整个列表中):
Similarly, about 90'000 elements are placed in Hsignal. The key problem is that python takes about 1.1 seconds to run this! I have even tried with this alternative, that removes some loops (here I still use arrays, as I have to add elements to an entire list):
start = time.time()
result = list()
for event in Ctrigger:
c = Csignal - event - offset
d = Csignal[abs(c) <= gate/2]
result.append(list(d))
flat = [item for sublist in result for item in sublist]
flat = np.sort(np.asarray(list(set(flat))))
end = time.time()
print(end-start)
但更糟的是，将近10分钟.
but it's even worse, almost 10 minutes.
我真的不明白问题出在哪里.对于我的应用程序， Ctrigger 长100e06，而 Csignal 长20e06.在matlab中，相同的代码需要1.06秒，而在python中则需要10分钟以上.似乎同时删除循环并加快处理速度并非易事.
I can't really understand where the problem is. For my application Ctrigger is 100e06 long, and Csignal around 20e06. In matlab the same code takes 1.06 seconds, against more than 10 minutes in python. It also seems that it is not straightforward to remove the loops and speeding the process at the same time.
编辑我:我已经介绍了我正在使用的Matlab代码以及一个可执行示例.我还列出了 Hsignal 列表，而 Ctrigger 和 Csignal 仍然是数组.结果:0.05s和6.5s
EDIT I: I have introduced the Matlab code I am using, as well as an executable example. I also made Hsignal a list, while Ctrigger and Csignal are still arrays. Result: 0.05s vs 6.5s
编辑II:现在，我仅使用RiccardoBucco建议的列表.结果:0.05s和1.5s
EDIT II: now I only use lists, as suggested by RiccardoBucco. Result: 0.05s vs 1.5s
编辑III:而不是附加到 Hsignal 上，我先声明它，然后更改单个元素，我注意到这带来了小的提高(即使看起来将 Hsignal 保留为数组会更快！).然后，我用主算法声明了一个函数.结果:0.05s和1.1s
EDIT III: instead of appending to Hsignal I am declaring it first, then changing individual elements, which I noticed brought a small speed up (even though it seems that keeping Hsignal as an array is faster!). Then I declared a function with the main algorithm. Result: 0.05s vs 1.1s
推荐答案
如何将运行时间降至6ms
您已经看到Python循环非常慢.默认情况下，没有jit编译器可以像Matlab中那样加速循环.因此，您有以下可能性:
How to get the runtime down to 6ms
As you already have seen Python loops are extremely slow. Per default there is no jit-Compiler which speeds up loops as in Matlab. So you have following possibilities:
- Vectorize your code in Numpy, if possible.
- Use Cython to compile the function
- Use Numba to compile the function
在下面的示例中，我使用Numba，因为在这种情况下使用起来非常简单.
In the following example I use Numba, because it is really simple to use in such cases.
示例
import numpy as np
import numba as nb
@nb.njit()
def main_nb(Ctrigger, Csignal, offset, gate):
Hsignal = np.zeros(Ctrigger.shape[0])
marker = 1
for j in range(Ctrigger.shape[0]):
m = marker
t_star = Ctrigger[j] + offset - gate/2
t_sto = Ctrigger[j] + offset + gate/2
while m < Csignal.shape[0]:
if (Csignal[m] < t_star):
marker = m
m = m + 1
elif (Csignal[m] >= t_star and Csignal[m] <= t_sto):
Hsignal[m] = Csignal[m]
m = m + 1
elif (Csignal[m] > t_sto):
break
return Hsignal
还请注意避免使用列表.像在Matlab中一样使用简单的数组.
Also note to avoid Lists if possible. Use simple arrays like you would do in Matlab.
时间
import time
#Use simple numpy arrays if possible, not lists
Ctrigger = np.arange(0, 3000000, 3)
length_t = Ctrigger.shape[0]
Bsignal = np.arange(0, 3000000, 10)
noise = 1e-05*np.random.rand(Bsignal.shape[0])
Csignal = np.sort(np.asarray(Bsignal) + noise)
offset = 3
gate = 1
start = time.time()
Hsignal=main(Ctrigger, Csignal, offset, gate)
print("Pure Python takes:" +str(time.time()-start))
#Pure Python takes:6.049151659011841
#First call takes longer (compilation overhead)
#The same may be the case in matlab
start = time.time()
Hsignal=main_nb(Ctrigger, Csignal, offset, gate)
print("First Numba run takes:" +str(time.time()-start))
#First Numba run takes:0.16272664070129395
start = time.time()
Hsignal=main_nb(Ctrigger, Csignal, offset, gate)
print("All further Numba calls run takes:" +str(time.time()-start))
#All further Numba calls run takes:0.006016731262207031
Hsignal = np.unique(Hsignal)
这篇关于将Matlab转换为Python-加速循环的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
