JuliaDiff / JuliaDiff/ForwardDiff.jl

`InexactError` when trying to differentiate `ODEProblem` using `LabelledArrays`

Open
#533 0 comments 0 reactions 0 assignees View on GitHub

Nobody has claimed this yet.

Dominant language
Julia
Stars
1k
Forks
160
PR merge metrics
No merged PRs in 30d

Description

It looks like the source of the problem is that recursive_unitless_bottom_eltype from RecursiveArrayTools is inferring an Int64 type rather than a Float64 type from the Dual of a LabelledArray. This is happening somewhere inside ODEProblem's __init.

Note that this works fine if you switch from LVector to Vector for the parameter/initial-condition arrays.

using LabelledArrays
using ModelingToolkit
using DifferentialEquations
using ForwardDiff


pars = @parameters  β α
@variables t E(t)
D = Differential(t)

eqs = [D(E) ~ α*t + β]
sys = ODESystem(eqs)

u0 = LVector(E = 1.0)
p = LVector(β = 0.8,
            α = 0.75)
tspan = (0.0,1.0)
sys = structural_simplify(sys)
prob = ODEProblem(sys,u0,tspan,p)

function loss(inputs)
    N = length(prob.p)
    _p = inputs[1:N]
    _u0 = inputs[N:end]
    tmp_sol = solve(prob,Tsit5(),p=_p,u0=_u0,saveat=0.1)
    sum(abs2,tmp_sol[E])
end

#Substitute a subset of the parameters/initial-conditions into the full set
inputter(par_vec) = LVector(β = p[1], α = par_vec[1], E = par_vec[2])
p0 = LVector(α = 1.0, E = 0.5)
dummy_loss(par_vec) = sum(abs2,par_vec)

ForwardDiff.gradient(x->dummy_loss(inputter(x)),p0) #works
ForwardDiff.gradient(dummy_loss,inputter(p0)) #Works
ForwardDiff.gradient(loss,inputter(p0)) #Works
ForwardDiff.gradient(x->loss(inputter(x)),p0) # Error

Error message:

ERROR: InexactError: Int64(161//1000)
Stacktrace:
  [1] Integer
    @ ./rational.jl:109 [inlined]
  [2] convert(#unused#::Type{Int64}, x::Rational{Int64})
    @ Base /Applications/Julia-1.6.app/Contents/Resources/julia/share/julia/base/number.jl:7
  [3] OrdinaryDiffEq.Tsit5ConstantCache(T::Type, T2::Type)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/a5qvy/src/tableaus/low_order_rk_tableaus.jl:610
  [4] alg_cache(alg::Tsit5, u::Vector{Real}, rate_prototype::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#65#66", Float64}, Float64, 2}}, uEltypeNoUnits::Type, uBottomEltypeNoUnits::Type, tTypeNoUnits::Type, uprev::Vector{Real}, uprev2::Vector{Real}, f::ODEFunction{true, ModelingToolkit.var"#f#165"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xbabc66d3, 0x633a919b, 0xe4683138, 0x0e254e2d, 0x36925fe8)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x7eb6766a, 0x02a244e1, 0x54678655, 0x9f7f1370, 0x31962eb1)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, t::Float64, dt::Float64, reltol::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#65#66", Float64}, Float64, 2}}, p::Vector{Real}, calck::Bool, #unused#::Val{true})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/a5qvy/src/caches/low_order_rk_caches.jl:348
  [5] __init(prob::ODEProblem{Vector{Real}, Tuple{Float64, Float64}, true, Vector{Real}, ODEFunction{true, ModelingToolkit.var"#f#165"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xbabc66d3, 0x633a919b, 0xe4683138, 0x0e254e2d, 0x36925fe8)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x7eb6766a, 0x02a244e1, 0x54678655, 0x9f7f1370, 0x31962eb1)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Float64, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/a5qvy/src/solve.jl:295
  [6] __solve(::ODEProblem{Vector{Real}, Tuple{Float64, Float64}, true, Vector{Real}, ODEFunction{true, ModelingToolkit.var"#f#165"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xbabc66d3, 0x633a919b, 0xe4683138, 0x0e254e2d, 0x36925fe8)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x7eb6766a, 0x02a244e1, 0x54678655, 0x9f7f1370, 0x31962eb1)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Tsit5; kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/a5qvy/src/solve.jl:4
  [7] solve_call(_prob::ODEProblem{Vector{Real}, Tuple{Float64, Float64}, true, Vector{Real}, ODEFunction{true, ModelingToolkit.var"#f#165"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xbabc66d3, 0x633a919b, 0xe4683138, 0x0e254e2d, 0x36925fe8)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x7eb6766a, 0x02a244e1, 0x54678655, 0x9f7f1370, 0x31962eb1)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5; merge_callbacks::Bool, kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/oe7VF/src/solve.jl:61
  [8] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, ModelingToolkit.var"#f#165"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xbabc66d3, 0x633a919b, 0xe4683138, 0x0e254e2d, 0x36925fe8)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x7eb6766a, 0x02a244e1, 0x54678655, 0x9f7f1370, 0x31962eb1)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Real}, p::Vector{Real}, args::Tsit5; kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/oe7VF/src/solve.jl:82
  [9] #solve#59
    @ ~/.julia/packages/DiffEqBase/oe7VF/src/solve.jl:70 [inlined]
 [10] loss(inputs::LArray{Real, 1, Vector{Real}, (:β, :α, :E)})
    @ Main ./REPL[227]:5
 [11] (::var"#65#66")(x::LArray{ForwardDiff.Dual{ForwardDiff.Tag{var"#65#66", Float64}, Float64, 2}, 1, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#65#66", Float64}, Float64, 2}}, (:α, :E)})
    @ Main ./REPL[234]:1
 [12] vector_mode_dual_eval
    @ ~/.julia/packages/ForwardDiff/QOqCN/src/apiutils.jl:37 [inlined]
 [13] vector_mode_gradient(f::var"#65#66", x::LArray{Float64, 1, Vector{Float64}, (:α, :E)}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#65#66", Float64}, Float64, 2, LArray{ForwardDiff.Dual{ForwardDiff.Tag{var"#65#66", Float64}, Float64, 2}, 1, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#65#66", Float64}, Float64, 2}}, (:α, :E)}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:106
 [14] gradient(f::Function, x::LArray{Float64, 1, Vector{Float64}, (:α, :E)}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#65#66", Float64}, Float64, 2, LArray{ForwardDiff.Dual{ForwardDiff.Tag{var"#65#66", Float64}, Float64, 2}, 1, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#65#66", Float64}, Float64, 2}}, (:α, :E)}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:19
 [15] gradient(f::Function, x::LArray{Float64, 1, Vector{Float64}, (:α, :E)}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#65#66", Float64}, Float64, 2, LArray{ForwardDiff.Dual{ForwardDiff.Tag{var"#65#66", Float64}, Float64, 2}, 1, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#65#66", Float64}, Float64, 2}}, (:α, :E)}}) (repeats 2 times)
    @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:17
 [16] top-level scope
    @ REPL[234]:1

Contributor guide

No contributing guide indexed for this repository

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

Research direction

The report points to recursive_unitless_bottom_eltype in RecursiveArrayTools and ODEProblem's init; begin by reproducing the supplied LVector/Vector comparison and follow the stack trace into OrdinaryDiffEq's Tsit5 cache. Done means the final ForwardDiff.gradient call completes without InexactError while the Vector case remains working.

Written by the indexing model from the issue text.

Assessment

Tech stack
julia
Domain
devtools
Issue type
Bug
Difficulty
4/5
Estimated time
3-5 days
Activity status
Stale
Clarity
Mostly clear
Newbie friendliness
35/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.