Julia 的运行速度比 python 慢一个数量级

Question

我试图将 python 代码移植到 Julia 中进行尝试（下面给出了两个代码）。 Julia 在我的机器上的运行速度比 python 慢了大约 10 倍。我究竟做错了什么？我对 Julia 很陌生，所以感谢任何帮助。

这是python代码：

import matplotlib.pyplot as plt
from numba import jit
from numpy import random
import time

N=1000
kplus=2
kminus=1
T=20
T_th=10
sdt=1
frac = 0.5
threshold = frac*N

@jit(nopython=True)
def run(kp, km):
    base=np.ones(N)
    mb=np.arange(N)
    m=N
    th=0
    time_data = np.zeros(int(T/sdt))
    histogram=np.zeros(N+1)
    time_data[0]=N
    time_temp = sdt
    while th<T:
        if m==0:
            #print(th)
            break
        
        if th>time_temp:
            time_data[int(time_temp/sdt)] = m
            if th>T_th:
                histogram[int(m)] += 1
            #time_data[int(time_temp/sdt)] = N if m>threshold else 0
            time_temp = time_temp + 1*sdt
            
        kt=m*(kp+km)
        th=th+random.exponential(1/kt)
        ran=kt*random.rand()
        index=int(ran/(kp+km))
        rem=ran-index*(kp+km)
        #print(rem)
        if rem<km:
            base[mb[index]]=0
            tmp=mb[index]
            mb[index]=mb[m-1]
            mb[m-1]=tmp
            m=m-1
                
        else:
            pos=random.randint(N)
            if base[pos]==0:
                base[pos]=1
                mb[m]=pos
                m=m+1
                
            
    return time_data, histogram
    
    
num_runs = 1000
time_data_avg = np.zeros(int(T/sdt))
td_var=np.zeros(int(T/sdt))
hist=np.zeros(N+1)

for _ in range(num_runs):
    m,l = run(2,1)
    hist += l
    time_data_avg += m/num_runs
    td_var += m*m/num_runs
td_var -= time_data_avg**2

这是我编写的相应 Julia 代码：

using Random
using Distributions
using Plots

N=1000
kplus=2
kminus=1
T=20
sdt=1
frac = 0.5
threshold = frac*N

function run(kp,km)
    base=fill(1,N)
    mb=collect(1:N)
    m=N
    th=0
    time_data = fill(0,floor(Int, T/sdt))
    time_data[1]=N
    time_temp = sdt
    
    while th<T
        # println(th, ' ', m)
        if m==0
            println(th)
            break
        end
        
        if th>time_temp
            time_data[ceil(Int, time_temp/sdt)+1]=m
            time_temp += sdt
        end
        
        kt=m*(kp+km)
        th=th+rand(Exponential(1/kt))
        ran=kt*rand(Float64)
        index=floor(Int,ran/(kp+km))
        rem=ran-index*(kp+km)
        index=index+1
        
        if rem<km
            base[mb[index]]=0
            tmp=mb[index]
            mb[index]=mb[m]
            mb[m]=tmp
            m=m-1
        else
            pos=rand(1:N)
            if base[pos]==0
                base[pos]=1
                mb[m+1]=pos
                m=m+1
            end
        end
        
    end
    return time_data
end


function sample(num_runs)
    time_data_avg = fill(0.0, floor(Int, T/sdt))
    td_var=fill(0.0, floor(Int, T/sdt))
    for i in 1:num_runs
        m = run(2,1)
        time_data_avg .+= m/num_runs
        td_var .+= m.*(m/num_runs)
    end
    td_var .-= time_data_avg.^2
    
    return time_data_avg, td_var
end

@time begin
   tm,tv=sample(1000)
end

对于 python 代码，我使用 bash time 命令测量时间。我还确保 numba 没有并行化。

Answer 1

我无法直接在我的机器上比较你的代码，因为你的 Python 代码没有正确执行产生时间。但是，在 Julia 中，无需更改代码即可轻松解决上述 @DNF 和 @Przemyslaw Szufel 提到的问题。 只需将所有内容包装在 function 或 let 块中即可。然后你不需要在你的函数中重写任何东西。我进行了这样的测试，结果是：

您的 Julia 代码“按原样”计时

 22.927925 seconds (360.97 M allocations: 5.731 GiB, 3.36% gc time, 1.49% compilation time)

（可以立即看到代码有问题，因为您有 361M 的分配，总计 5.7GiB 的 RAM）

如果我将它包装在一个函数中，你的代码的时间安排

我所做的只是添加两行：

function fun()

在顶部和一个额外的

end

在底部。

时间是：

julia> fun()
  0.779523 seconds (5.00 k allocations: 16.144 MiB, 0.91% gc time)

（分配显着下降；这可以进一步优化，但我知道您想比较两种语言的完全相同的代码）

所以你得到 22.9 / 0.78 = 29.3 的加速。其中，鉴于您报告的内容应该比 Python 快 3 倍左右。

小的编码风格注释：我需要更改的是在 function sample(num_runs) 之前删除一个空行，因为通常 Julia REPL 假定传递两个空行结束一个定义（并且在将所有内容包装在一个函数中时，您不希望它的定义完成） .

话虽如此 - 这不是一个人通常会像其他人已经评论过的那样编写 Julia 程序的方式。我建议您查看 Julia 手册的 Performance tips 部分了解详细信息。

编辑

这是对代码的快速而肮脏的重写，它执行的分配更少并且速度更快：

function f()
    N=1000
    T=20
    sdt=1
    base=fill(1,N)
    mb=collect(1:N)
    time_data = fill(0,floor(Int, T/sdt))

    function run(kp,km)
        fill!(base, 1)
        mb .= 1:N
        fill!(time_data, 0)
        m=N
        th=0
        time_data[1]=N
        time_temp = sdt
        
        @inbounds while th<T
            # println(th, ' ', m)
            if m==0
                println(th)
                break
            end
            
            if th>time_temp
                time_data[ceil(Int, time_temp/sdt)+1]=m
                time_temp += sdt
            end
            
            kt=m*(kp+km)
            th=th+rand(Exponential(1/kt))
            ran=kt*rand(Float64)
            index=floor(Int,ran/(kp+km))
            rem=ran-index*(kp+km)
            index=index+1
            
            if rem<km
                base[mb[index]]=0
                tmp=mb[index]
                mb[index]=mb[m]
                mb[m]=tmp
                m=m-1
            else
                pos=rand(1:N)
                if base[pos]==0
                    base[pos]=1
                    mb[m+1]=pos
                    m=m+1
                end
            end
            
        end
        return time_data
    end

    function sample(num_runs)
        time_data_avg = fill(0.0, floor(Int, T/sdt))
        td_var=fill(0.0, floor(Int, T/sdt))
        for i in 1:num_runs
            m = run(2,1)
            time_data_avg .+= m ./ num_runs
            td_var .+= m.*(m ./ num_runs)
        end
        td_var .-= time_data_avg.^2
        
        return time_data_avg, td_var
    end

    @time begin
        tm,tv=sample(1000)
    end
end

它确实：

julia> f();
  0.739664 seconds (2 allocations: 448 bytes)

虽然使用您的原始代码相同的是：

julia> f();
  0.778454 seconds (5.00 k allocations: 16.144 MiB)

（所以时间上的差异不是很大但很明显，但分配明显更低，这意味着对 GC 的压力更小；我没有分析你的代码逻辑 - 我只是改变了内存管理）