python - 将 Matlab 转换为 Python - 加速循环

Question

我一直在将一些代码从 Matlab 翻译成 Python，我们用来分析实验室中的数据。我们有两个时间戳列表，我们想用一个来预示另一个:对于第一个列表中的每个元素，我们在第二个列表中寻找在时间上有精确分隔的时间戳。如果有，我们将它们放在单独的列表中。

这是我正在使用的那种带有随机数据的 Matlab 代码的可运行示例。它可能非常粗糙，因为我不精通 Matlab。下面的 Ctrigger 是触发列表，而 Csignal 是我们想要预告的信号列表。对于 Ctrigger 的每个元素，我们查看 Csignal 中是否有元素位于以偏移为中心的窗口内，并且具有宽度门。选定的事件将被放置在 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 秒来运行它。之所以引入标记计​​数器，是因为这两个列表Csignal 和Ctrigger 区域已经及时排序了。标记设置在一个预兆窗口的开始处:当我移动到下一个触发器时，我不会再次查看所有 Csignal，而只会从该窗口的开始处查看。为了避免重复计算，我在最后删除了重复项。

如果您想了解代码，这里是输入和输出的简化版本:

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。按照一些建议，我首先声明包含主算法的函数，然后调用它:

# 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)

类似地，大约 90'000 个元素放置在 Hsignal 中。关键问题是python运行这个大约需要1.1秒!我什至尝试过这种替代方法，它删除了一些循环(这里我仍然使用数组，因为我必须将元素添加到整个列表中):

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分钟。

我真的无法理解问题出在哪里。对于我的应用程序，Ctrigger 的长度为 100e06，而 Csignal 的长度约为 20e06。在 matlab 中，相同的代码需要 1.06 秒，而在 python 中则需要 10 多分钟。同时消除循环和加速过程似乎也并不简单。

编辑我:我已经介绍了我正在使用的 Matlab 代码，以及一个可执行示例。我也把 Hsignal 做成了一个列表，而 Ctrigger 和 Csignal 仍然是数组。结果:0.05s vs 6.5s

编辑二:现在我只使用列表，正如 RiccardoBucco 所建议的那样。结果:0.05s vs 1.5s

编辑三:而不是附加到 Hsignal 我首先声明它，然后更改单个元素，我注意到这带来了一个小速度(即使将 Hsignal 保留为数组似乎更快!)。然后我用主算法声明了一个函数。结果:0.05s vs 1.1s

Answer 1

可能会减慢您的算法速度的是使用 np.append在

Hsignal = np.append(Hsignal, Csignal[m])

您应该使用列表，而不是 NumPy 数组:

Ctrigger = [1, 10, 11, 20, 30, 40, 50, 60]
Csignal = [4, 11, 13, 17, 25, 34, 41, 42, 50, 57, 65]

offset = 2
gate = 2

Hsignal = []
marker = 0

for j in range(len(Ctrigger)):
    m = marker
    t_start = Ctrigger[j] + offset - gate/2
    t_stop = Ctrigger[j] + offset + gate/2   
    while m < len(Csignal):   
        if Csignal[m] < t_start:
            marker = m
            m = m + 1
        elif Csignal[m] <= t_stop:
            Hsignal.append(Csignal[m])
            m = m + 1
        else:
            break

Hsignal = sorted(set(Hsignal))

列表构建完成后，您可以将其转换为数组:

Hsignal = np.array(Hsignal)