Type-instability arising from initial condition

[fgittins]

Type-instability arising from initial condition

[lhapp27]

I just encountered a similar problem, however for SecondOrderODEProblem. Also, similarly to your second example, this occurs already before, and therefore independently of, passing the solution further into another call of ODEProblem. Please note that I am literally using this pacakge for 1 day, so I cannot rule out incorrect usage from my own side.

MWE:

function fun!(ddu, du, u, p, r)
    ddu[1] = -exp(-r^2)*u[1]
end
function main()
    du0 = [1.0]
    u0 = [0.0]
    tspan = (0.0,20.0)
    prob = SecondOrderODEProblem(fun!,du0,u0,tspan)
    sol = solve(prob, Tsit5())
end

Output of @code_warntype:

@code_warntype main()
MethodInstance for main()
  from main() @ Main REPL[3]:1
Arguments
  #self#::Core.Const(main)
Locals
  sol::Any
  prob::Any
  tspan::Tuple{Float64, Float64}
  u0::Vector{Float64}
  du0::Vector{Float64}
Body::Any
1 ─      (du0 = Base.vect(1.0))
│        (u0 = Base.vect(0.0))
│        (tspan = Core.tuple(0.0, 20.0))
│        (prob = Main.SecondOrderODEProblem(Main.fun!, du0, u0, tspan::Core.Const((0.0, 20.0))))
│   %5 = prob::Any
│   %6 = Main.Tsit5()::Core.Const(Tsit5(; stage_limiter! = trivial_limiter!, step_limiter! = trivial_limiter!, thread = static(false),))
│   %7 = Main.solve(%5, %6)::Any
│        (sol = %7)
└──      return %7

Further info:
julia version 1.9.4

packages:

using Pkg; Pkg.status()
  [be6e5d0e] Antique v0.3.1
  [6e4b80f9] BenchmarkTools v1.5.0
  [34da2185] Compat v4.14.0
  [5a033b19] CurveFit v0.6.0
  [a93c6f00] DataFrames v1.6.1
  [864edb3b] DataStructures v0.18.18
  [31c24e10] Distributions v0.25.107
  [ffbed154] DocStringExtensions v0.9.3
  [6912e4f1] Espresso v0.6.1
  [442a2c76] FastGaussQuadrature v1.0.2
  [5789e2e9] FileIO v1.16.3
  [c58ffaec] FortranFiles v0.6.0
  [92c85e6c] GSL v1.0.1
  [c27321d9] Glob v1.3.1
  [19dc6840] HCubature v1.6.0
  [34004b35] HypergeometricFunctions v0.3.23
  [a98d9a8b] Interpolations v0.15.1
  [4138dd39] JLD v0.13.4
  [033835bb] JLD2 v0.4.46
  [98e50ef6] JuliaFormatter v1.0.54
  [b964fa9f] LaTeXStrings v1.3.1
  [2fda8390] LsqFit v0.15.0
  [23992714] MAT v0.10.6
  [eff96d63] Measurements v2.11.0
  [d41bc354] NLSolversBase v7.8.3
⌅ [8913a72c] NonlinearSolve v3.4.0
  [6fe1bfb0] OffsetArrays v1.13.0
  [67456a42] OhMyThreads v0.5.0
  [429524aa] Optim v1.9.3
⌅ [1dea7af3] OrdinaryDiffEq v6.66.0
  [d96e819e] Parameters v0.12.3
  [793d2195] PartialWaveFunctions v0.2.0
  [a03496cd] PlotlyBase v0.8.19
  [f0f68f2c] PlotlyJS v0.18.13
⌅ [91a5bcdd] Plots v1.39.0
  [c46f51b8] ProfileView v1.7.2
  [1fd47b50] QuadGK v2.9.4
  [295af30f] Revise v3.5.14
  [f2b01f46] Roots v2.1.5
  [efcf1570] Setfield v1.1.1
  [276daf66] SpecialFunctions v2.3.1
  [90137ffa] StaticArrays v1.9.3
  [82ae8749] StatsAPI v1.7.0
  [2913bbd2] StatsBase v0.34.2
  [fd094767] Suppressor v0.2.6
  [bd369af6] Tables v1.11.1
  [0796e94c] Tokenize v0.5.28
  [770da0de] UpdateJulia v0.4.4
  [8ba89e20] Distributed
  [de0858da] Printf
  [1a1011a3] SharedArrays
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

[fgittins]

The type-instability for your ODE is resolved if the function is declared isinplace or not (see the documentation):

function main()
   du0 = [1.0]
   u0 = [0.0]
   tspan = (0.0,20.0)
   prob = SecondOrderODEProblem{true}(fun!,du0,u0,tspan)
   sol = solve(prob, Tsit5())
end

With this change, all types are now fully inferred:

julia> @code_warntype main()
MethodInstance for main()
  from main() @ Main REPL[3]:1
Arguments
  #self#::Core.Const(main)
Locals
  sol::ODESolution{Float64, 2, Vector{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}}, ODEProblem{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, DynamicalODEFunction{true, SciMLBase.FullSpecialize, ODEFunction{true, SciMLBase.FullSpecialize, typeof(fun!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#313#315", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SecondOrderODEProblem{true}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DynamicalODEFunction{true, SciMLBase.FullSpecialize, ODEFunction{true, SciMLBase.FullSpecialize, typeof(fun!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#313#315", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}, Vector{Float64}, Vector{Vector{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}}, Nothing, OrdinaryDiffEq.Tsit5Cache{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Nothing}, SciMLBase.DEStats, Nothing}
  prob::ODEProblem{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, DynamicalODEFunction{true, SciMLBase.FullSpecialize, ODEFunction{true, SciMLBase.FullSpecialize, typeof(fun!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#313#315", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SecondOrderODEProblem{true}}
  tspan::Tuple{Float64, Float64}
  u0::Vector{Float64}
  du0::Vector{Float64}
Body::ODESolution{Float64, 2, Vector{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}}, ODEProblem{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, DynamicalODEFunction{true, SciMLBase.FullSpecialize, ODEFunction{true, SciMLBase.FullSpecialize, typeof(fun!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#313#315", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SecondOrderODEProblem{true}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DynamicalODEFunction{true, SciMLBase.FullSpecialize, ODEFunction{true, SciMLBase.FullSpecialize, typeof(fun!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#313#315", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}, Vector{Float64}, Vector{Vector{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}}, Nothing, OrdinaryDiffEq.Tsit5Cache{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Nothing}, SciMLBase.DEStats, Nothing}
1 ─       (du0 = Base.vect(1.0))
│         (u0 = Base.vect(0.0))
│         (tspan = Core.tuple(0.0, 20.0))
│   %4  = Core.apply_type(Main.SecondOrderODEProblem, true)::Core.Const(SecondOrderODEProblem{true})
│   %5  = Main.fun!::Core.Const(fun!)
│   %6  = du0::Vector{Float64}
│   %7  = u0::Vector{Float64}
│         (prob = (%4)(%5, %6, %7, tspan::Core.Const((0.0, 20.0))))
│   %9  = prob::ODEProblem{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, DynamicalODEFunction{true, SciMLBase.FullSpecialize, ODEFunction{true, SciMLBase.FullSpecialize, typeof(fun!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#313#315", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SecondOrderODEProblem{true}}
│   %10 = Main.Tsit5()::Core.Const(Tsit5(; stage_limiter! = trivial_limiter!, step_limiter! = trivial_limiter!, thread = static(false),))
│   %11 = Main.solve(%9, %10)::ODESolution{Float64, 2, Vector{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}}, ODEProblem{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, DynamicalODEFunction{true, SciMLBase.FullSpecialize, ODEFunction{true, SciMLBase.FullSpecialize, typeof(fun!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#313#315", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SecondOrderODEProblem{true}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DynamicalODEFunction{true, SciMLBase.FullSpecialize, ODEFunction{true, SciMLBase.FullSpecialize, typeof(fun!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#313#315", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}, Vector{Float64}, Vector{Vector{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}}, Nothing, OrdinaryDiffEq.Tsit5Cache{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Nothing}, SciMLBase.DEStats, Nothing}
│         (sol = %11)
└──       return %11

[lhapp27]

Ah thank you, I am sorry to have highjacked your issue. I really thought it's a similar problem.

[fgittins]

I have found a way to remove the type-instability. If I change the initial condition passed to prob2 to be 2.0 rather than 2,

prob2 = ODEProblem{false}(perturbation, 2.0, rspan, (n, K, sol1))

then the types are fully inferred.

In retrospect, this should have been obvious. I was passing an Int64 as the initial value, when it needed to be a Float64. This is also why the other example I considered above is type-unstable.

However, I am surprised that Julia's compiler was unable to infer from the ODE and the other arguments that the initial condition needed to be promoted to a Float64. If this isn't possible, I am happy to close this issue.

[ChrisRackauckas]

It should be possible. It is odd that promote_type just gives up here.