SciML / DifferentialEquations.jl

Multi-language suite for high-performance solvers of differential equations and scientific machine learning (SciML) components. Ordinary differential equations (ODEs), stochastic differential equations (SDEs), delay differential equations (DDEs), differential-algebraic equations (DAEs), and more in Julia.
https://docs.sciml.ai/DiffEqDocs/stable/
Other
2.87k stars 229 forks source link

Return type for automatic solver selection is not inferrable #655

Open imciner2 opened 4 years ago

imciner2 commented 4 years ago

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.

ChrisRackauckas commented 4 years ago

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.