Performance repo

repo.jl

function get_f(version=:paper)
    κ = 0.3

    η = 2
    b1 = 0.0028388
    A_0 = 0.02722^(1/(1-κ))
    g_A0 = 0.0092/(1-0.3)
    σ_0 = 0.13418
    M_pre = 596.4
    δ_B = 0.01
    s = 3.0

    if version==:code
        η = 2
        b1 = 0.0028388
        A_0 = 0.02722^(1/(1-κ))
        g_A0 = 0.0092/(1-0.3)
        σ_0 = 0.13418
        M_pre = 596.4
        δ_B = 0.01
        s = 3.0
    end

    b2 = 2
    δ_u = 0.015
    L_0 = 6514
    L_∞ = 8600
    g_L_star = 0.035
    K_0 = 137
    δ_K = 0.1
    δ_A = 0.001
    g_σ0 = -0.0073
    δ_σ = 0.003
    a0 = 1.17
    a1 = 2
    a2 = 2.8
    g_Ψ_star = -0.005
    T_0 = 0.76
    M_0 = 808.9
    δ_M0 = 0.014
    δ_M∞ = 0.004
    δ_M_star = 0.01
    B_0 = 1.1
    η_forc = 3.8
    λ = η_forc / s
    EF_0 = -0.06
    EF_100 = 0.3
    σ_forc = 0.032
    σ_ocean = 0.007
    Δt = 1.
    ζ = 0.02

    L(t) = L_0 + (L_∞-L_0) * (1-exp(-g_L_star * t))
    A(t) = A_0 * exp(g_A0 * (1-exp(-δ_A * t))/δ_A)
    σ(t) = σ_0 * exp(g_σ0 * (1-exp(-δ_σ * t))/δ_σ)
    Ψ(t) = σ(t) / a2 * a0 * (1- (1-exp(g_Ψ_star * t))/a1)
    g_L(t) = (log(L(t+Δt))-log(L(t)))/Δt
    g_A(t) = (log(A(t+Δt))-log(A(t)))/Δt
    B(t) = B_0 * exp(-δ_B * t)
    EF(t) = EF_0 + 0.01 * (EF_100-EF_0) * min(t,100)
    δ_M(t) = δ_M∞ + (δ_M0-δ_M∞) * exp(-δ_M_star * t)
    ΔT(t) = max(0.7 + 0.02t - 0.00007 * t^2,0)

    function transition(state,state_new, x)
        c = x[1]
        Λ = x[2]
        k = state[1]
        M = state[2]
        τ = state[3]
        T = state[4]

        t = -log(1-τ)/ζ
        μ = (Λ/Ψ(t))^(1/a2)
        y = (1 - Ψ(t) * μ^a2) / (1+b1*T^b2) * k^κ
        k_new = ((1-δ_K)^Δt * k + y * Δt - c * Δt) * exp(-(g_A(t) + g_L(t)) * Δt)
        E = (1-μ) * σ(t) * A(t) * L(t) * k^κ + B(t)
        M_new = M_pre + (M - M_pre) * (1-δ_M(t))^Δt + E * Δt
        T_new = (1-σ_forc) * T + σ_forc * (η_forc * log(M_new/M_pre)/log(2.) + EF(t)) / λ - σ_ocean * ΔT(t)
        t_new = t + 1.
        τ_new = 1 - exp(-ζ*t_new)

        state_new[1] = k_new
        state_new[2] = M_new
        state_new[3] = τ_new
        state_new[4] = T_new
    end

	return transition
end

function main(n=10_000_000)
    f = get_f()

    state_old = [0.5, 700., 0., 0.] * 1.1
    state_new = ones(4)
    x = ones(2)

    f(state_old, state_new, x)

    @time for i=1:n
        f(state_old, state_new, x)
    end
end

# make sure everything is compiled
main(1)
#Real perf test
main()

# Use this to see type instabilities
# code_warntype(f, (Vector{Float64}, Vector{Float64}, Vector{Float64}))