I am working on a problem where I simulate an ODE as part of a larger optimization scheme (the ODE becomes part of the cost function computation in a derivative-free setting), and I was trying to track down some inefficiencies in the code that was causing it to have what seems like excessive runtime. I looked at the solver wrapper that I am using, and it appears that it is having some type instability when creating the ODEProblem and the return type from the solve routine.
I tried this on the simple Lorenz demo from the docs, and it is also showing type instability. Is there a better way to create these objects so that the compiler is able to assign concrete types to them and not Any?
Here is the code I used for the Lorenz sample:
using DifferentialEquations
function lorenz!(du,u,p,t)
du[1] = 10.0*(u[2]-u[1])
du[2] = u[1]*(28.0-u[3]) - u[2]
du[3] = u[1]*u[2] - (8/3)*u[3]
end
function lorenzSim!()
u0 = [1.0;0.0;0.0]
tspan = (0.0,50.0)
prob = ODEProblem(lorenz!,u0,tspan)
sol = solve(prob)
u0 = sol[end];
tspan = (50.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan)
sol = solve(prob)
return sol[end]
end
Here is the type warning output, note that prob, sol, and u0 are all Any.
julia> @code_warntype lorenzSim!([1.0;0.0;0.0])
Variables
#self#::Core.Compiler.Const(lorenzSim!, false)
u::Array{Float64,1}
u0::Any
tspan::Tuple{Float64,Float64}
prob::Any
sol::Any
Body::Any
1 ─ (u0 = u)
│ (tspan = Core.tuple(0.0, 50.0))
│ (prob = Main.ODEProblem(Main.lorenz!, u0::Array{Float64,1}, tspan::Core.Compiler.Const((0.0, 50.0), false)))
│ (sol = Main.solve(prob))
│ %5 = sol::Any
│ %6 = Base.lastindex(sol)::Any
│ (u0 = Base.getindex(%5, %6))
│ (tspan = Core.tuple(50.0, 100.0))
│ (prob = Main.ODEProblem(Main.lorenz!, u0, tspan::Core.Compiler.Const((50.0, 100.0), false)))
│ (sol = Main.solve(prob))
│ %11 = sol::Any
│ %12 = Base.lastindex(sol)::Any
│ %13 = Base.getindex(%11, %12)::Any
└── return %13