# 一段简单代码的优化问题

``````using LinearAlgebra
using SpecialFunctions
using BenchmarkTools

function main(θ::Array{Float64,1}, ε::Array{Float64,1}, t::Array{Float64,1}, h::Float64)
α = θ[1]
b = θ[2]
σ = zeros(size(ε))
ω = cumprod(vcat(1, (1 .- (α.+1)./(1:length(ε)-1))))

for i = 1:length(t)
w_t = ω[1:i]
σ = σ .- σ[1]
y_t = σ[i:-1:1]
σ[i] = h^(-α)*dot(w_t, y_t)
end
σ = @. b*(σ + ε*t^(-α) / gamma(1-α))
return σ
end

α = 0.3
h = 0.0001
t = collect(0:h:2)
θ = [α, 5.0]
ε = 1.0 .* ones(size(t))

@btime main(θ, ε, t, h)`
``````

``````3.983 s (115938 allocations: 5.97 GiB)
``````

``````using LinearAlgebra
using SpecialFunctions
using BenchmarkTools

mutable struct Params
#
α :: Float64
b :: Float64
#
ω :: Array{Float64,1}
#
h :: Float64
t :: Array{Float64,1}
ε :: Array{Float64,1}
σ :: Array{Float64,1}
end

"""
Initial parameters.
"""
function initialparameters(params::Array{Float64,1}, ε::Array{Float64,1}, t::Array{Float64,1}, h::Float64)
α = params[1]
b = params[2]
σ = zeros(size(ε))
ω = cumprod(vcat(1, (1 .- (α.+1)./(1:length(ε)-1))))
Params(α, b, ω, h, t, ε, σ)
end

"""
Kernel function for computing σ(t = tn)
"""
function computesigma!(p::Params, i::Int64)
w_t = p.ω[1:i]
p.σ = p.σ .- p.σ[1]
y_t = p.σ[i:-1:1]
p.σ[i] = p.h^(-p.α)*dot(w_t, y_t)
end

"""
Main function
"""
function main(θ::Array{Float64,1}, ε::Array{Float64,1}, t::Array{Float64,1}, h::Float64)
p = initialparameters(θ, ε, t, h)
for i in 1:length(p.t)
computesigma!(p, i)
end
p.σ = @. p.b*(p.σ + p.ε*p.t^(-p.α) / gamma(1-p.α))
p.σ
end

α = 0.3
h = 0.0001
t = collect(0:h:2)
θ = [α, 5.0]
ε = 1.0 .* ones(size(t))

@btime main(θ, ε, t, h)
``````

``````3.960 s (115937 allocations: 5.97 GiB)
``````

``````- p.σ = @. p.b*(p.σ + p.ε*p.t^(-p.α) / gamma(1-p.α))
+ @. p.σ = p.b*(p.σ + p.ε*p.t^(-p.α) / gamma(1-p.α))
``````

``````julia> function main(θ::Array{Float64,1}, ε::Array{Float64,1}, t::Array{Float64,1}, h::Float64)
α = θ[1]
b = θ[2]
σ = zeros(size(ε))
ω = cumprod(vcat(1, (1 .- (α.+1)./(1:length(ε)-1))))

@inbounds for i = 1:length(t)
w_t = @view ω[1:i]
σ .-= σ[1]
y_t = @view σ[i:-1:1]
σ[i] = h^(-α)*dot(w_t, y_t)
end
@. σ = b*(σ + ε*t^(-α) / gamma(1-α))
return σ
end
main (generic function with 1 method)

julia> @btime main(θ, ε, t, h);
492.309 ms (25 allocations: 626.14 KiB)
``````