DifferentialEquations.jl

DifferentialEquations.jl copied to clipboard

Using an automatic solver to solve this sample problem does not allow for the compiler to infer the type of sol1 or sol2 at all.

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!(u::Array{Float64,1})
    tspan = (0.0,50.0)
    prob1 = ODEProblem{true}(lorenz!,u,tspan)
    sol1 = solve(prob1)
    u1 = sol1[end]

    tspan = (50.0,100.0)
    prob2 = ODEProblem{true}(lorenz!,u1,tspan)
    sol2 = solve(prob2)
    return sol2[end]
end

Result of @code_warntype

julia> @code_warntype lorenzSim!([1.0;0.0;0.0])
Variables
  #self#::Core.Compiler.Const(lorenzSim!, false)
  u::Array{Float64,1}
  prob1::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}
  sol1::Any
  u1::Any
  tspan::Tuple{Float64,Float64}
  prob2::ODEProblem{_A,Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem} where _A
  sol2::Any

Body::Any
1 ─       (tspan = Core.tuple(0.0, 50.0))
│   %2  = Core.apply_type(Main.ODEProblem, true)::Core.Compiler.Const(ODEProblem{true,tType,isinplace,P,F,K,PT} where PT where K where F where P where isinplace where tType, false)
│         (prob1 = (%2)(Main.lorenz!, u, tspan::Core.Compiler.Const((0.0, 50.0), false)))
│         (sol1 = Main.solve(prob1::Core.Compiler.PartialStruct(ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, Any[Core.Compiler.Const(ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}(lorenz!, LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing), false), Array{Float64,1}, Core.Compiler.Const((0.0, 50.0), false), Core.Compiler.Const(DiffEqBase.NullParameters(), false), Core.Compiler.Const(Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}(), false), Core.Compiler.Const(DiffEqBase.StandardODEProblem(), false)])))
│   %5  = sol1::Any
│   %6  = Base.lastindex(sol1)::Any
│         (u1 = Base.getindex(%5, %6))
│         (tspan = Core.tuple(50.0, 100.0))
│   %9  = Core.apply_type(Main.ODEProblem, true)::Core.Compiler.Const(ODEProblem{true,tType,isinplace,P,F,K,PT} where PT where K where F where P where isinplace where tType, false)
│   %10 = u1::Any
│         (prob2 = (%9)(Main.lorenz!, %10, tspan::Core.Compiler.Const((50.0, 100.0), false)))
│         (sol2 = Main.solve(prob2))
│   %13 = sol2::Any
│   %14 = Base.lastindex(sol2)::Any
│   %15 = Base.getindex(%13, %14)::Any
└──       return %15

Changing it to use a solver such as Tsit5() makes the type inferrable, but using a solver such as Rosenbrock23() doesn't allow inference of the type of the solution either.

In theory this could infer with constant propagation. In practice this will be hard. The ODE default algorithm is inherently dependent on the values of keyword arguments, so without total constant specialization this won't occur.