Unexpected InexactError() while using DifferentialEquations.jl


#1

Playing around with DifferentialEquations.jl, I ran into a simple corner case. The code corresponds to a trivial ODE, which is a modification of the following (that works just fine).

function f(u,p,t)
    return p * u
end

u0=0.5
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan,1.01)
sol = solve(prob)

If I try to implement an in-place version, however, I get an InexactError()

function g(du,u,p,t)
    du = p*u
end

u0=0.5
tspan = (0.0,1.0)
prob = ODEProblem(g,u0,tspan,1.01)
sol = solve(prob)

Here is the full stack trace

InexactError()

Stacktrace:
 [1] convert(::Type{Int64}, ::Float64) at ./float.jl:679
 [2] zeros at ./array.jl:263 [inlined]
 [3] zeros at ./array.jl:264 [inlined]
 [4] zeros at ./array.jl:266 [inlined]
 [5] alg_cache(::OrdinaryDiffEq.Tsit5, ::Float64, ::Float64, ::Type{T} where T, ::Type{T} where T, ::Type{T} where T, ::Float64, ::Float64, ::#g, ::Float64, ::Float64, ::Float64, ::Float64, ::Bool, ::Type{Val{true}}) at /home/jrun/.julia/v0.6/OrdinaryDiffEq/src/caches/low_order_rk_caches.jl:399
 [6] #init#1331(::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Void, ::Bool, ::Void, ::Bool, ::Bool, ::Void, ::Bool, ::Bool, ::Float64, ::Bool, ::Rational{Int64}, ::Void, ::Void, ::Int64, ::Rational{Int64}, ::Int64, ::Int64, ::Rational{Int64}, ::Bool, ::Int64, ::Rational{Int64}, ::Rational{Int64}, ::Int64, ::Float64, ::Float64, ::DiffEqBase.#ODE_DEFAULT_NORM, ::DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN, ::DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::String, ::DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE, ::Void, ::Bool, ::Bool, ::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.ODEProblem{Float64,Float64,true,Float64,#g,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/jrun/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:240
 [7] (::DiffEqBase.#kw##init)(::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.ODEProblem{Float64,Float64,true,Float64,#g,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0
 [8] #solve#1330(::Array{Any,1}, ::Function, ::DiffEqBase.ODEProblem{Float64,Float64,true,Float64,#g,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/jrun/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:6
 [9] (::DiffEqBase.#kw##solve)(::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqBase.ODEProblem{Float64,Float64,true,Float64,#g,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0 (repeats 2 times)
 [10] #solve#2(::Bool, ::Array{Any,1}, ::Function, ::DiffEqBase.ODEProblem{Float64,Float64,true,Float64,#g,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::Void) at /home/jrun/.julia/v0.6/DifferentialEquations/src/default_solve.jl:14
 [11] (::DiffEqBase.#kw##solve)(::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqBase.ODEProblem{Float64,Float64,true,Float64,#g,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::Void) at ./<missing>:0
 [12] #solve#1(::Bool, ::Array{Any,1}, ::Function, ::DiffEqBase.ODEProblem{Float64,Float64,true,Float64,#g,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}) at /home/jrun/.julia/v0.6/DifferentialEquations/src/default_solve.jl:5
 [13] solve(::DiffEqBase.ODEProblem{Float64,Float64,true,Float64,#g,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}) at /home/jrun/.julia/v0.6/DifferentialEquations/src/default_solve.jl:2
 [14] include_string(::String, ::String) at ./loading.jl:522

Everything goes through (with minor syntax adjustments) if the ODE is 2-dimensional or higher dimensional.

Any suggestions as to what is the issue here?


#2

You cannot write in-place functions with immutable types. If you do this, notice that du isn’t actually updated but the original value is just shadowed. Since it’s not updated, this g function doesn’t actually do anything (try it). Because this isn’t actually doing anything, the cache construction isn’t setup to be compatible with immutable types.

Though we could probably use a better error message.


#3

Thanks, that clarifies things. I had forgotten about the importance of immutability for certain types in Julia, and was expecting a pass-by-reference behavior.

A better error message would have certainly been helpful, but I can see how it is a subtle issue.

What do you suggest is the best way to flag the issue? Make a DifferentialEquations.jl issue about the error message?


#4

Immutable structs in Julia are value-types so they don’t have references, which is why they can’t pass-by-reference (this is a requirement to make then fast). You can manually wrap it in a Ref but that’s usually not a good idea if you don’t need to. DiffEq has the out-of-place route because immutables can be much more efficient if treated without references.

Yup, please feel free to start that thread. There’s a few other places where we could use a more informative error as well.


#5

Every argument in Julia gets passed-by-value.