using QuantEcon
using BasisMatrices
using ContinuousDPs
using PyPlot
using Random
alpha = 0.4
beta = 0.96
mu = 0
sigma = 0.1;
f(s, x) = log(x)
g(s, x, e) = (s - x)^alpha * e;
shock_size = 250
shocks = exp.(mu .+ sigma * randn(shock_size))
weights = fill(1/shock_size, shock_size);
grid_max = 4.
n = 30
s_min, s_max = 1e-5, grid_max
basis = Basis(ChebParams(n, s_min, s_max))
1 dimensional Basis on the hypercube formed by (1.0e-5,) × (4.0,). Basis families are Cheb
x_lb(s) = s_min
x_ub(s) = s;
ab = alpha * beta
c1 = log(1 - ab) / (1 - beta)
c2 = (mu + alpha * log(ab)) / (1 - alpha)
c3 = 1 / (1 - beta)
c4 = 1 / (1 - ab)
# True optimal policy
c_star(y) = (1 - alpha * beta) * y
# True value function
v_star(y) = c1 + c2 * (c3 - c4) + c4 * log(y);
cdp = ContinuousDP(f, g, beta, shocks, weights, x_lb, x_ub, basis);
@code_warntype ContinuousDP(f, g, beta, shocks, weights, x_lb, x_ub, basis)
Body::ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)} │ 60 1 ─ %1 = invoke ContinuousDPs.Interp(_9::Basis{1,Tuple{ChebParams{Float64}}})::ContinuousDPs.Interp{1,Array{Float64,1},Array{Float64,2},LinearAlgebra.LU{Float64,Array{Float64,2}}} │╻╷ Type61 │ %2 = %new(ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)}, f, g, discount, shocks, weights, x_lb, x_ub, %1)::ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)} │ 62 └── return %2
C_star = cdp.interp.Phi \ v_star.(cdp.interp.S)
Tv = Array{Float64}(undef, cdp.interp.length)
C = copy(C_star)
bellman_operator!(cdp, C, Tv);
@code_warntype bellman_operator!(cdp, C, Tv)
Body::Array{Float64,1} │╻ getproperty198 1 ─ %1 = (Base.getfield)(cdp, :interp)::ContinuousDPs.Interp{1,Array{Float64,1},TM,TL} where TL<:LinearAlgebra.Factorization where TM<:(AbstractArray{T,2} where T) ││ │ %2 = (Base.getfield)(%1, :S)::Array{Float64,1} │ │ %3 = invoke ContinuousDPs.s_wise_max!(_2::ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)}, %2::Array{Float64,1}, _3::Array{Float64,1}, _4::Array{Float64,1})::Array{Float64,1} │╻ getproperty199 │ %4 = (Base.getfield)(cdp, :interp)::ContinuousDPs.Interp{1,Array{Float64,1},TM,TL} where TL<:LinearAlgebra.Factorization where TM<:(AbstractArray{T,2} where T) ││ │ %5 = (Base.getfield)(%4, :Phi_lu)::LinearAlgebra.Factorization │ │ (ContinuousDPs.ldiv!)(C, %5, %3) │ 200 └── return C
grid_size = 200
grid_y = collect(range(s_min, stop=s_max, length=grid_size))
V_approx = funeval(C, cdp.interp.basis, grid_y);
fig, ax = subplots(figsize=(9, 5))
ax[:set_ylim](-35, -24)
ax[:plot](grid_y, V_approx, lw=2, alpha=0.6, label=L"$Tv^*$")
ax[:plot](grid_y, v_star.(grid_y), lw=2, alpha=0.6, label=L"$v^*$")
ax[:legend](loc="lower right")
show()
@time bellman_operator!(cdp, C, Tv)
@time bellman_operator!(cdp, C, Tv)
@time bellman_operator!(cdp, C, Tv);
0.106662 seconds (1.21 M allocations: 54.788 MiB, 19.19% gc time) 0.083827 seconds (1.13 M allocations: 51.089 MiB, 11.62% gc time) 0.091252 seconds (1.20 M allocations: 54.406 MiB, 13.58% gc time)
s = 2.
@code_warntype ContinuousDPs._s_wise_max(cdp, s, C)
Body::Tuple{Float64,Float64} │╻ getproperty151 1 ── %1 = (Base.getfield)(cdp, :shocks)::Array{Float64,1} │╻ size │ %2 = (Base.arraysize)(%1, 1)::Int64 ││╻ Type │ %3 = $(Expr(:foreigncall, :(:jl_alloc_array_2d), Array{Float64,2}, svec(Any, Int64, Int64), :(:ccall), 3, Array{Float64,2}, :(%2), 1, 1, :(%2)))::Array{Float64,2} │ 152 │ %4 = %new(getfield(ContinuousDPs, Symbol("#objective#3")){ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)},Float64,Array{Float64,1},Array{Float64,2}}, cdp, s, C, %3)::getfield(ContinuousDPs, Symbol("#objective#3")){ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)},Float64,Array{Float64,1},Array{Float64,2}} │ 161 │ %5 = Optim.optimize::Core.Compiler.Const(Optim.optimize, false) │╻ x_lb │ %6 = Main.s_min::Any │ │ %7 = (isa)(%6, Float64)::Bool │ └─── goto #12 if not %7 │ 2 ── %9 = π (%6, Float64) ││╻ sqrt │ %10 = (Base.Math.sqrt_llvm)(2.220446049250313e-16)::Float64 ││╻ #optimize#78 │ %11 = π (1, Int64) │││ │ %12 = π (%11, Int64) │││ │ %13 = π (false, Bool) │││ │ %14 = π (%12, Int64) │││ │ invoke Core.kwfunc(Optim.optimize::Any) │││ │ %16 = Optim.optimize::typeof(Optim.optimize) │││╻╷╷╷╷ #optimize │ %17 = (Base.slt_int)(0, 1)::Bool ││││┃│││ isempty └─── goto #4 if not %17 │││││┃││ iterate 3 ── goto #5 ││││││┃│ iterate 4 ── invoke Base.getindex(()::Tuple{}, 1::Int64) │││││││┃ iterate └─── $(Expr(:unreachable)) │││││││ 5 ┄─ goto #6 │││││╻ iterate 6 ── goto #7 │││││ 7 ── goto #8 ││││ 8 ── %25 = invoke Optim.:(#optimize#71)(%10::Float64, 2.220446049250313e-16::Float64, 1000::Int64, false::Bool, %13::Bool, nothing::Nothing, %14::Int64, false::Bool, %16::typeof(Optim.optimize), %4::getfield(ContinuousDPs, Symbol("#objective#3")){ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)},Float64,Array{Float64,1},Array{Float64,2}}, %9::Float64, _3::Float64, $(QuoteNode(Optim.Brent()))::Optim.Brent)::Optim.UnivariateOptimizationResults{Float64,Float64,_1,_2,Float64,Optim.Brent} where _2 where _1 ││││ └─── goto #9 │││ 9 ── goto #10 ││ 10 ─ goto #11 │ 11 ─ goto #13 │ 12 ─ %30 = (%5)(%4, %6, s)::Optim.UnivariateOptimizationResults{Float64,Float64,_1,_2,Float64,Optim.Brent} where _2 where _1 │ └─── goto #13 │ 13 ┄ %32 = φ (#11 => %25, #12 => %30)::Optim.UnivariateOptimizationResults{Float64,Float64,_1,_2,Float64,Optim.Brent} where _2 where _1 │╻ getproperty162 │ %33 = (Base.getfield)(%32, :minimum)::Any │ │ (Core.typeassert)(%33, ContinuousDPs.Float64) │ │ %35 = π (%33, Float64) │╻ - │ %36 = (Base.neg_float)(%35)::Float64 │╻ getproperty163 │ %37 = (Base.getfield)(%32, :minimizer)::Any │ │ (Core.typeassert)(%37, ContinuousDPs.Float64) │ │ %39 = π (%37, Float64) │ 164 │ %40 = (Core.tuple)(%36, %39)::Tuple{Float64,Float64} │ └─── return %40
@time ContinuousDPs._s_wise_max(cdp, s, C)
@time ContinuousDPs._s_wise_max(cdp, s, C)
@time ContinuousDPs._s_wise_max(cdp, s, C);
0.002056 seconds (28.31 k allocations: 1.281 MiB) 0.002404 seconds (28.29 k allocations: 1.278 MiB) 0.002289 seconds (28.29 k allocations: 1.278 MiB)
v_init_func(s) = 5 * log(s)
w = v_init_func.(grid_y)
n = 35
fig, ax = subplots(figsize=(9, 6))
ax[:set_ylim](-50, 10)
ax[:set_xlim](minimum(grid_y), maximum(grid_y))
lb = "initial condition"
jet = ColorMap("jet")
ax[:plot](grid_y, w, color=jet(0), lw=2, alpha=0.6, label=lb)
S = cdp.interp.S
V = v_init_func.(S)
for i in 1:n
C = cdp.interp.Phi \ V
bellman_operator!(cdp, C, V)
w = funeval(C, cdp.interp.basis, grid_y)
ax[:plot](grid_y, w, color=jet(i / n), lw=2, alpha=0.6)
end
lb = "true value function"
ax[:plot](grid_y, v_star.(grid_y), "k-", lw=2, alpha=0.8, label=lb)
ax[:legend](loc="lower right")
show()
res = solve(cdp);
Compute iterate 6 with error 1.0658141036401503e-13 Converged in 6 steps
@code_warntype solve(cdp)
Body::ContinuousDPs.CDPSolveResult{PFI,1,Array{Float64,1},Array{Float64,1}} │ 296 1 ─ %1 = ContinuousDPs.PFI::Core.Compiler.Const(PFI, false) │╻ solve │ %2 = (Base.Math.sqrt_llvm)(2.220446049250313e-16)::Float64 ││ │ %3 = invoke ContinuousDPs.:(#solve#7)(%2::Float64, 500::Int64, 2::Int64, 50::Int64, _1::Function, _2::ContinuousDP{1,Array{Float64,1},Array{Float64,1},typeof(f),typeof(g),typeof(x_lb),typeof(x_ub)}, %1::Type{PFI})::ContinuousDPs.CDPSolveResult{PFI,1,Array{Float64,1},Array{Float64,1}} │ └── return %3
@time res = solve(cdp)
@time res = solve(cdp)
@time res = solve(cdp);
Compute iterate 6 with error 1.0658141036401503e-13 Converged in 6 steps 0.909546 seconds (11.46 M allocations: 571.261 MiB, 11.57% gc time) Compute iterate 6 with error 1.0658141036401503e-13 Converged in 6 steps 0.914953 seconds (11.46 M allocations: 571.270 MiB, 11.24% gc time) Compute iterate 6 with error 1.0658141036401503e-13 Converged in 6 steps 0.910241 seconds (11.46 M allocations: 571.270 MiB, 11.05% gc time)
set_eval_nodes!(res, grid_y);
@code_warntype set_eval_nodes!(res, grid_y)
Body::ContinuousDPs.CDPSolveResult{PFI,1,Array{Float64,1},Array{Float64,1}} │╻╷ set_eval_nodes!137 1 ─ %1 = (Base.getfield)(s_nodes_coord, 1, true)::Array{Float64,1} ││╻ setproperty! │ (Base.setfield!)(res, :eval_nodes, %1) ││╻ setproperty! │ (Base.setfield!)(res, :eval_nodes_coord, s_nodes_coord) │╻ set_eval_nodes! │ %4 = invoke ContinuousDPs.evaluate!(_2::ContinuousDPs.CDPSolveResult{PFI,1,Array{Float64,1},Array{Float64,1}})::ContinuousDPs.CDPSolveResult{PFI,1,Array{Float64,1},Array{Float64,1}} │ └── return %4
@time set_eval_nodes!(res, grid_y)
@time set_eval_nodes!(res, grid_y)
@time set_eval_nodes!(res, grid_y);
0.575157 seconds (7.74 M allocations: 349.491 MiB, 11.62% gc time) 0.573510 seconds (7.74 M allocations: 349.491 MiB, 10.07% gc time) 0.592195 seconds (7.74 M allocations: 349.491 MiB, 9.38% gc time)
fig, ax = subplots(figsize=(9, 5))
ax[:set_ylim](-35, -24)
ax[:plot](grid_y, res.V, lw=2, alpha=0.6, label="approximate value function")
ax[:plot](grid_y, v_star.(grid_y), lw=2, alpha=0.6, label="true value function")
ax[:legend](loc="lower right")
show()
fig, ax = subplots(figsize=(9, 5))
ax[:plot](grid_y, res.X, lw=2, alpha=0.6, label="approximate policy function")
ax[:plot](grid_y, c_star.(grid_y), lw=2, alpha=0.6, label="true policy function")
ax[:legend](loc="lower right")
show()
fig, ax = subplots(figsize=(9, 5))
ax[:plot](grid_y, res.resid, lw=2, alpha=0.6, label="residual")
show()
s_init = 0.1
ts_length = 100
y = simulate(res, s_init, ts_length)
fig, ax = subplots(figsize=(9, 6))
ax[:plot](1:ts_length, y, lw=2, alpha=0.6, label="beta = $(cdp.discount)" )
ax[:legend](loc="lower right")
show()
fig, ax = subplots(figsize=(9, 6))
for beta in (0.9, 0.94, 0.98)
cdp.discount = beta
res = solve(cdp, verbose=0)
set_eval_nodes!(res, grid_y)
y = simulate(res, s_init, ts_length)
ax[:plot](1:ts_length, y, lw=2, alpha=0.6, label="beta = $(cdp.discount)" )
end
ax[:legend](loc="lower right")
show()
@code_warntype simulate!(MersenneTwister(0), Array{Float64}(undef, ts_length), res, s_init)
Body::Array{Float64,1} │╻ size335 1 ── %1 = (Base.arraysize)(s_path, 1)::Int64 │╻ getproperty336 │ %2 = (Base.getfield)(res, :cdp)::ContinuousDP{1,Array{Float64,1},Array{Float64,1},Tf,Tg,Tlb,Tub} where Tub<:Function where Tlb<:Function where Tg<:Function where Tf<:Function ││ │ %3 = (Base.getfield)(%2, :weights)::Array{Float64,1} │╻ cumsum │ invoke Core.kwfunc(Base.cumsum::Any) ││ │ %5 = Base.cumsum::typeof(cumsum) ││╻╷╷╷╷ #cumsum │ %6 = (Base.slt_int)(0, 1)::Bool │││┃│││ isempty └─── goto #3 if not %6 ││││┃││ iterate 2 ── goto #4 │││││┃│ iterate 3 ── invoke Base.getindex(()::Tuple{}, 1::Int64) ││││││┃ iterate └─── $(Expr(:unreachable)) ││││││ 4 ┄─ goto #5 ││││╻ iterate 5 ── goto #6 ││││ 6 ── goto #7 │││ 7 ── %14 = invoke Base.:(#cumsum#572)(1::Int64, %5::Function, %3::Array{Float64,1})::Array{Float64,1} │││ └─── goto #8 ││ 8 ── goto #9 │╻ -337 9 ── %17 = (Base.sub_int)(%1, 1)::Int64 ││╻╷╷╷ rand │ %18 = $(Expr(:foreigncall, :(:jl_alloc_array_1d), Array{Float64,1}, svec(Any, Int64), :(:ccall), 2, Array{Float64,1}, :(%17), :(%17)))::Array{Float64,1} │││╻╷ rand! │ %19 = (Base.arraylen)(%18)::Int64 ││││╻ rand! │ %20 = (Base.mul_int)(8, %19)::Int64 │││││╻ _rand! │ %21 = (Base.arraylen)(%18)::Int64 ││││││╻ * │ %22 = (Base.mul_int)(8, %21)::Int64 ││││││╻ <= │ %23 = (Base.sle_int)(%20, %22)::Bool ││││││ └─── goto #11 if not %23 ││││││ 10 ─ goto #12 │ 11 ─ nothing ││││││ 12 ┄ %27 = φ (#10 => true, #11 => false)::Bool ││││││ └─── goto #14 if not %27 ││││││╻ macro expansion 13 ─ %29 = $(Expr(:gc_preserve_begin, :(%18))) │││││││╻╷ pointer │ %30 = $(Expr(:foreigncall, :(:jl_array_ptr), Ptr{Float64}, svec(Any), :(:ccall), 1, :(%18)))::Ptr{Float64} │││││││╻ Type │ %31 = %new(Random.UnsafeView{Float64}, %30, %19)::Random.UnsafeView{Float64} │││││││ │ invoke Random.rand!(_2::MersenneTwister, %31::Random.UnsafeView{Float64}, $(QuoteNode(Random.SamplerTrivial{Random.CloseOpen01{Float64},Float64}(Random.CloseOpen01{Float64}())))::Random.SamplerTrivial{Random.CloseOpen01{Float64},Float64}) │││││││ │ $(Expr(:gc_preserve_end, :(%29))) │││││││ └─── goto #15 ││││││╻ Type 14 ─ %35 = %new(Core.AssertionError, "sizeof(Float64) * n64 <= sizeof(T) * length(A) && isbitstype(T)")::AssertionError ││││││ │ (Base.throw)(%35) ││││││ └─── $(Expr(:unreachable)) │││││ 15 ┄ goto #16 ││││ 16 ─ goto #17 │││ 17 ─ goto #18 ││ 18 ─ goto #19 │╻ -338 19 ─ %42 = (Base.sub_int)(%1, 1)::Int64 ││╻ Type │ %43 = $(Expr(:foreigncall, :(:jl_alloc_array_1d), Array{Int64,1}, svec(Any, Int64), :(:ccall), 2, Array{Int64,1}, :(%42), :(%42)))::Array{Int64,1} │╻ -339 │ %44 = (Base.sub_int)(%1, 1)::Int64 ││╻╷╷╷ Type │ %45 = (Base.sle_int)(1, %44)::Bool │││╻ unitrange_last │ (Base.sub_int)(%44, 1) ││││ │ %47 = (Base.ifelse)(%45, %44, 0)::Int64 ││╻╷╷ isempty │ %48 = (Base.slt_int)(%47, 1)::Bool ││ └─── goto #21 if not %48 ││ 20 ─ goto #22 ││ 21 ─ goto #22 │ 22 ┄ %52 = φ (#20 => true, #21 => false)::Bool │ │ %53 = φ (#21 => 1)::Int64 │ │ %54 = φ (#21 => 1)::Int64 │ │ %55 = (Base.not_int)(%52)::Bool │ └─── goto #38 if not %55 │ 23 ┄ %57 = φ (#22 => %53, #37 => %92)::Int64 │ │ %58 = φ (#22 => %54, #37 => %93)::Int64 │╻ getindex340 │ %59 = (Base.arrayref)(true, %18, %57)::Float64 ││╻╷╷╷╷ #searchsortedlast#5 │ %60 = (Base.arraysize)(%14, 1)::Int64 │││╻╷╷╷ searchsortedlast │ %61 = (Base.slt_int)(%60, 0)::Bool ││││┃│││││ axes │ %62 = (Base.ifelse)(%61, 0, %60)::Int64 ││││╻╷ searchsortedlast └─── %63 = (Base.add_int)(%62, 1)::Int64 ││││╻ searchsortedlast 24 ┄ %64 = φ (#23 => 0, #28 => %78)::Int64 │││││ │ %65 = φ (#23 => %63, #28 => %79)::Int64 │││││╻ - │ %66 = (Base.sub_int)(%65, 1)::Int64 │││││╻ < │ %67 = (Base.slt_int)(%64, %66)::Bool │││││ └─── goto #29 if not %67 │││││╻ + 25 ─ %69 = (Base.add_int)(%64, %65)::Int64 ││││││╻ >>> │ %70 = (Base.lshr_int)(%69, 0x0000000000000001)::Int64 ││││││╻ << │ %71 = (Base.shl_int)(%69, 0xffffffffffffffff)::Int64 ││││││ │ %72 = (Base.ifelse)(true, %70, %71)::Int64 │││││╻ getindex │ %73 = (Base.arrayref)(false, %14, %72)::Float64 ││││││╻ isless │ %74 = (Base.fpislt)(%59, %73)::Bool │││││ └─── goto #27 if not %74 │││││ 26 ─ goto #28 │ 27 ─ nothing │││││ 28 ┄ %78 = φ (#26 => %64, #27 => %72)::Int64 │││││ │ %79 = φ (#26 => %72, #27 => %65)::Int64 │││││ └─── goto #24 │││││ 29 ─ goto #30 ││││ 30 ─ goto #31 │││ 31 ─ goto #32 ││ 32 ─ goto #33 │╻ + 33 ─ %85 = (Base.add_int)(%64, 1)::Int64 │╻ setindex! │ (Base.arrayset)(true, %43, %85, %57) ││╻ == │ %87 = (%58 === %47)::Bool ││ └─── goto #35 if not %87 ││ 34 ─ goto #36 ││╻ + 35 ─ %90 = (Base.add_int)(%58, 1)::Int64 │╻ iterate └─── goto #36 │ 36 ┄ %92 = φ (#35 => %90)::Int64 │ │ %93 = φ (#35 => %90)::Int64 │ │ %94 = φ (#34 => true, #35 => false)::Bool │ │ %95 = (Base.not_int)(%94)::Bool │ └─── goto #38 if not %95 │ 37 ─ goto #23 │╻ getproperty343 38 ─ %98 = (Base.getfield)(res, :eval_nodes_coord)::Tuple{Array{Float64,1}} ││╻ getindex │ %99 = (Base.getfield)(%98, 1, true)::Array{Float64,1} ││╻ Type │ %100 = invoke LinParams{Array{Float64,1}}(%99::Array{Float64,1}, 0::Int64)::LinParams{Array{Float64,1}} ││ │ %101 = (Core.tuple)(%100)::Tuple{LinParams{Array{Float64,1}}} ││╻╷╷ _Basis │ %102 = %new(Basis{1,Tuple{LinParams{Array{Float64,1}}}}, %101)::Basis{1,Tuple{LinParams{Array{Float64,1}}}} │╻ getproperty344 │ %103 = (Base.getfield)(res, :X)::Array{Float64,1} ││╻ Type │ %104 = invoke BasisMatrices.BasisMatrix(BasisMatrices.Nothing::Type{Nothing}, %102::Basis{1,Tuple{LinParams{Array{Float64,1}}}}, $(QuoteNode(Tensor()))::Tensor)::BasisMatrix{Tensor,SparseArrays.SparseMatrixCSC{Float64,Int64}} ││╻ Type │ %105 = invoke BasisMatrices.get_coefs(%102::Basis{1,Tuple{LinParams{Array{Float64,1}}}}, %104::BasisMatrix{Tensor,SparseArrays.SparseMatrixCSC{Float64,Int64}}, %103::Array{Float64,1})::Array{Float64,1} │╻╷ axes346 │ %106 = (Base.arraysize)(s_path, 1)::Int64 ││╻╷╷╷ map │ %107 = (Base.slt_int)(%106, 0)::Bool │││┃││ Type │ (Base.ifelse)(%107, 0, %106) │╻ getproperty347 │ %109 = (Base.getfield)(res, :cdp)::ContinuousDP{1,Array{Float64,1},Array{Float64,1},Tf,Tg,Tlb,Tub} where Tub<:Function where Tlb<:Function where Tg<:Function where Tf<:Function ││ │ %110 = (Base.getfield)(%109, :shocks)::Array{Float64,1} ││╻ size │ %111 = (Base.arraysize)(%110, 1)::Int64 │││╻╷╷╷ Type │ %112 = (Base.slt_int)(%111, 0)::Bool ││││┃│ Type │ (Base.ifelse)(%112, 0, %111) │╻ setindex!348 │ (Base.arrayset)(true, s_path, s_init, 1) │╻ -349 │ %115 = (Base.sub_int)(%1, 1)::Int64 ││╻╷╷╷ Type │ %116 = (Base.sle_int)(1, %115)::Bool │││╻ unitrange_last │ (Base.sub_int)(%115, 1) ││││ │ %118 = (Base.ifelse)(%116, %115, 0)::Int64 ││╻╷╷ isempty │ %119 = (Base.slt_int)(%118, 1)::Bool ││ └─── goto #40 if not %119 ││ 39 ─ goto #41 ││ 40 ─ goto #41 │ 41 ┄ %123 = φ (#39 => true, #40 => false)::Bool │ │ %124 = φ (#40 => 1)::Int64 │ │ %125 = φ (#40 => 1)::Int64 │ │ %126 = (Base.not_int)(%123)::Bool │ └─── goto #47 if not %126 │ 42 ┄ %128 = φ (#41 => %124, #46 => %149)::Int64 │ │ %129 = φ (#41 => %125, #46 => %150)::Int64 │╻ getindex350 │ %130 = (Base.arrayref)(true, s_path, %128)::Float64 │╻╷╷╷╷╷╷ Interpoland351 │ %131 = $(Expr(:foreigncall, :(:jl_alloc_array_2d), Array{Float64,2}, svec(Any, Int64, Int64), :(:ccall), 3, Array{Float64,2}, 1, 1, 1, 1))::Array{Float64,2} ││┃│││ Interpoland │ %132 = invoke Base.fill!(%131::Array{Float64,2}, %130::Float64)::Array{Float64,2} │││┃ funeval │ %133 = invoke BasisMatrices.funeval(%105::Array{Float64,1}, %102::Basis{1,Tuple{LinParams{Array{Float64,1}}}}, %132::Array{Float64,2}, 0::Int64)::Array{Float64,2} ││││╻ getindex │ %134 = (Base.arrayref)(true, %133, 1)::Float64 │╻ getindex352 │ %135 = (Base.arrayref)(true, %43, %128)::Int64 │╻ getproperty │ %136 = (Base.getfield)(res, :cdp)::ContinuousDP{1,Array{Float64,1},Array{Float64,1},Tf,Tg,Tlb,Tub} where Tub<:Function where Tlb<:Function where Tg<:Function where Tf<:Function ││ │ %137 = (Base.getfield)(%136, :shocks)::Array{Float64,1} │╻ getindex │ %138 = (Base.arrayref)(true, %137, %135)::Float64 │╻ +353 │ %139 = (Base.add_int)(%128, 1)::Int64 │╻ getproperty │ %140 = (Base.getfield)(res, :cdp)::ContinuousDP{1,Array{Float64,1},Array{Float64,1},Tf,Tg,Tlb,Tub} where Tub<:Function where Tlb<:Function where Tg<:Function where Tf<:Function ││ │ %141 = (Base.getfield)(%140, :g)::Function │ │ %142 = (%141)(%130, %134, %138)::Any │ │ (Base.setindex!)(s_path, %142, %139) ││╻ == │ %144 = (%129 === %118)::Bool ││ └─── goto #44 if not %144 ││ 43 ─ goto #45 ││╻ + 44 ─ %147 = (Base.add_int)(%129, 1)::Int64 │╻ iterate └─── goto #45 │ 45 ┄ %149 = φ (#44 => %147)::Int64 │ │ %150 = φ (#44 => %147)::Int64 │ │ %151 = φ (#43 => true, #44 => false)::Bool │ │ %152 = (Base.not_int)(%151)::Bool │ └─── goto #47 if not %152 │ 46 ─ goto #42 │ 356 47 ─ return s_path
@time simulate!(MersenneTwister(0), Array{Float64}(undef, ts_length), res, s_init)
@time simulate!(MersenneTwister(0), Array{Float64}(undef, ts_length), res, s_init)
@time simulate!(MersenneTwister(0), Array{Float64}(undef, ts_length), res, s_init);
0.000365 seconds (3.56 k allocations: 291.109 KiB) 0.000315 seconds (3.55 k allocations: 289.953 KiB) 0.000329 seconds (3.55 k allocations: 289.953 KiB)