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

レポートは RecursiveArrayTools の recursive_unitless_bottom_eltype と ODEProblem の init を指摘している。まず、提示された LVector/Vector の比較を再現し、スタックトレースを OrdinaryDiffEq の Tsit5 キャッシュまで追ってください。Vector のケースが引き続き動作したまま、最終的な ForwardDiff.gradient 呼び出しが InexactError なしで完了すれば完了です。

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