对 R、Python 和 Julia 中的简单水文学模块进行基准测试

Question

我最近开始与一个团队合作，该团队的成员不同程度地喜欢 Python、R 和 Julia（按受欢迎程度排序）。我们正在将由不同团队成员构建的一组水文学模块翻译成一个公共库，我想在就我们应该为此选择哪种语言达成共识之前做一些基本的基准测试。问题是我对 Python 的了解很差，而且对 Julia 还很陌生，而且我相信由于糟糕的编码实践，我可能会导致结果出现偏差。

这些模块（代码如下）的结果是：Python（93 毫秒）、R（55 毫秒）和 Julia（0.5 毫秒）。代码无法矢量化（根据我的理解），所以我假设 Julia 最快，但不是 100 倍，而且我还假设 Python 会比 R 更快。

谁能指出一些效率方面的基本错误？

Python：

import numpy as np
def SM_routine_f(infiltration=None, PET=None, rcScl=0.1, rcExp=1.3, PETexp=2, Zr=1000, n=0.5, smhp=0.00, smwp=0.10, smfc=0.25, s0=0.5):
    nt = len(infiltration)
    SM_store = np.empty(nt+1)
    SM_store[0] = s0
    smhp_stor = Zr * smhp
    smwp_stor = Zr * smwp
    smfc_stor = Zr * smfc
    max_stor = Zr * n
    for i in range(nt):
        thisSMstor = SM_store[i]
        AET = np.where(thisSMstor > smhp_stor, (thisSMstor - smhp_stor) * PET[i] * (thisSMstor / max_stor) ** PETexp / max_stor, 0)
        thisSMstor -= AET
        deepDrainage = np.where(thisSMstor > smfc_stor, rcScl * thisSMstor * (thisSMstor / smfc_stor) ** rcExp - (thisSMstor - smfc_stor), 0)
        SM_store[i+1] = min(max_stor, thisSMstor + infiltration[i] - deepDrainage)
    return SM_store / n
exInfil = np.random.normal(loc=1, scale=1, size=365*20)
exPET = np.random.normal(loc=2, scale=1, size=365*20)
import timeit
timeit.timeit(lambda: SM_routine_f(infiltration=exInfil, PET=exPET), number=1)

R：

SM_routine_f = function(infiltration = NA, PET = NA, 
    rcScl = 0.1,    # water retention curve scalar
    rcExp = 1.3,    # water retention curve exponent
    PETexp = 2,     # exponent of PET decay
    Zr = 1000,      # root depth
    n = 0.5,        # soil porosity
    smhp = 0.00,    # soil moisture at hydroscopic point
    smwp = 0.10,    # soil moisture at wilting point
    smfc = 0.25,    # soil moisture field capacity
    s0 = .5)    # initial soil moisture 
    {
    nt = length(infiltration)
    SM_store = c(s0, rep(NA, nt))
    smhp_stor = Zr * smhp
    smwp_stor = Zr * smwp
    smfc_stor = Zr * smfc
    max_stor = Zr * n
    
    for(i in 1:length(infiltration))    {
        thisSMstor = SM_store[i]
        AET = ifelse(thisSMstor > smhp_stor, 
            min(thisSMstor - smhp_stor, PET[i] * (thisSMstor / max_stor) ^ PETexp),
            0)
        thisSMstor = thisSMstor - AET
        
        deepDrainage = ifelse(thisSMstor > smfc_stor,
            min(thisSMstor - smfc_stor, rcScl * thisSMstor * (thisSMstor / smfc_stor)^rcExp),
            0)
        
        SM_store[i+1] = min(max_stor, thisSMstor - min(c(thisSMstor, deepDrainage)) + infiltration[i]) 
        }
    return(SM_store / n)
    }

    # generating dummy data for testing, e.g. 20 years of climate data
exInfil = rnorm(365*20, mean=1, sd=1)
exPET = rnorm(365*20, mean=2, sd=1)

    # benchmarking
microbenchmark::microbenchmark(SM_routine_f(infiltration=exInfil, PET = exPET))

朱莉娅：

function SM_routine_f(infiltration::Vector{Float64}, PET::Vector{Float64})

rcScl = 0.1
rcExp = 1.3
PETexp = 2
Zr = 1000.0
n = 0.5
smhp = 0.00
smwp = 0.10
smfc = 0.25
s0 = 0.5
nt = length(infiltration)
SM_store = [s0; fill(missing,nt)]
smhp_stor = Zr * smhp
smwp_stor = Zr * smwp
smfc_stor = Zr * smfc
max_stor = Zr * n

for i in 1:length(infiltration)
    thisSMstor = SM_store[i]
    AET = (thisSMstor > smhp_stor) ? min(thisSMstor - smhp_stor, PET[i] * (thisSMstor / max_stor) ^ PETexp) : 0
    thisSMstor = thisSMstor - AET
    deepDrainage = (thisSMstor > smfc_stor) ?  min(thisSMstor - smfc_stor, rcScl * thisSMstor * (thisSMstor / smfc_stor)^rcExp) : 0
    SM_store[i+1] = min(max_stor, thisSMstor - min(thisSMstor, deepDrainage) + infiltration[i])
end
return(SM_store / n)
end

@benchmark SM_routine_f($exInfiltration, $exPET)

Answer 1

总的来说，100 倍大约是在将 Julia（或其他快速语言）与 Python/R 进行比较时所期望的结果。