Tracers in new versions of ADNLPModels trigger errors with DifferentialEquations

[nathanemac]

When defining an ADNLPModel with DifferentialEquations

using ADNLPModels   # default version
using DifferentialEquations

function FH_ODE(dx, x, p, t)
  V, W = x
  I, μ, a, b, c = p
  dx[1] = (V - V^3 / 3 - W + I) / μ
  dx[2] = μ * (a * V - b * W + c)
end

function simulate(p)
  temp_prob = DifferentialEquations.remake(prob_FH, p=p)
  sol = DifferentialEquations.solve(
    temp_prob,
    DifferentialEquations.Vern9(),
  )
  return vec(sol)
end

function residual(p)
  F = simulate(p)
  F .-= data
  return F
end

function misfit(p)
  F = residual(p)
  return dot(F, F) / 2
end

u0 = [2.0; 0.0]
tspan = (0.0, 20.0)
savetime = 0.2
pars_FH = [0.5, 0.08, 1.0, 0.8, 0.7]
prob_FH = DifferentialEquations.ODEProblem(FH_ODE, u0, tspan, pars_FH)

The following line will fail if ADNLPModels is > 0.7.2

ADNLPModels.ADNLPModel(misfit, ones(5))

gives :

ERROR: TypeError: in typeassert, expected Float64, got a value of type SparseConnectivityTracer.HessianTracer{SparseConnectivityTracer.DictHessianPattern{Int64, BitSet, Dict{Int64, BitSet}, SparseConnectivityTracer.NotShared}}

And this line:

nequ = 202
ADNLPModels.ADNLSModel(residual, ones(5), nequ)

also fails:

ERROR: TypeError: in typeassert, expected Float64, got a value of type SparseConnectivityTracer.HessianTracer{SparseConnectivityTracer.DictHessianPattern{Int64, BitSet, Dict{Int64, BitSet}, SparseConnectivityTracer.NotShared}}

In both cases, the line causing the error is line 7 in my code:

function FH_ODE(dx, x, p, t)

Here is the complete error report for the ADNLPModel

ERROR: TypeError: in typeassert, expected Float64, got a value of type SparseConnectivityTracer.HessianTracer{SparseConnectivityTracer.DictHessianPattern{Int64, BitSet, Dict{Int64, BitSet}, SparseConnectivityTracer.NotShared}}
Stacktrace:
  [1] setindex!(A::Vector{…}, x::SparseConnectivityTracer.HessianTracer{…}, i::Int64)
    @ Base ./array.jl:987
  [2] FH_ODE(dx::Vector{Float64}, x::Vector{Float64}, p::Vector{SparseConnectivityTracer.HessianTracer{…}}, t::Float64)
    @ Main ~/Desktop/PhD/Recherche/iR2N/Coding/debug/debug RegularizedProblems/test_patch.jl:7
  [3] (::SciMLBase.Void{typeof(FH_ODE)})(::Vector{Float64}, ::Vararg{Any})
    @ SciMLBase ~/.julia/packages/SciMLBase/Zrnkf/src/utils.jl:486
  [4] (::FunctionWrappers.CallWrapper{…})(f::SciMLBase.Void{…}, arg1::Vector{…}, arg2::Vector{…}, arg3::Vector{…}, arg4::Float64)
    @ FunctionWrappers ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:65
  [5] macro expansion
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:137 [inlined]
  [6] do_ccall
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:125 [inlined]
  [7] FunctionWrapper
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:144 [inlined]
  [8] _call
    @ ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:12 [inlined]
  [9] FunctionWrappersWrapper
    @ ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:10 [inlined]
 [10] ODEFunction
    @ ~/.julia/packages/SciMLBase/Zrnkf/src/scimlfunctions.jl:2470 [inlined]
 [11] ode_determine_initdt(u0::Vector{…}, t::Float64, tdir::Float64, dtmax::Float64, abstol::Float64, reltol::Float64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), prob::ODEProblem{…}, integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/initdt.jl:53
 [12] auto_dt_reset!
    @ ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/integrators/integrator_interface.jl:450 [inlined]
 [13] handle_dt!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/solve.jl:647
 [14] __init(prob::ODEProblem{…}, alg::Vern9{…}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{…}; saveat::Tuple{}, 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::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{…}, abstol::Nothing, reltol::Nothing, qmin::Rational{…}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{…}, 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), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias::ODEAliasSpecifier, initializealg::OrdinaryDiffEqCore.DefaultInit, kwargs::@Kwargs{})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/solve.jl:609
 [15] __init (repeats 5 times)
    @ ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/solve.jl:11 [inlined]
 [16] #__solve#62
    @ ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/solve.jl:6 [inlined]
 [17] __solve
    @ ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/solve.jl:1 [inlined]
 [18] #solve_call#35
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:635 [inlined]
 [19] solve_call
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:592 [inlined]
 [20] #solve_up#44
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1142 [inlined]
 [21] solve_up
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1120 [inlined]
 [22] #solve#42
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1057 [inlined]
 [23] solve(prob::ODEProblem{…}, args::Vern9{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1047
 [24] simulate(p::Vector{SparseConnectivityTracer.HessianTracer{SparseConnectivityTracer.DictHessianPattern{…}}})
    @ Main ~/Desktop/PhD/Recherche/iR2N/Coding/debug/debug RegularizedProblems/test_patch.jl:13
 [25] residual(p::Vector{SparseConnectivityTracer.HessianTracer{SparseConnectivityTracer.DictHessianPattern{…}}})
    @ Main ~/Desktop/PhD/Recherche/iR2N/Coding/debug/debug RegularizedProblems/test_patch.jl:21
 [26] misfit(p::Vector{SparseConnectivityTracer.HessianTracer{SparseConnectivityTracer.DictHessianPattern{…}}})
    @ Main ~/Desktop/PhD/Recherche/iR2N/Coding/debug/debug RegularizedProblems/test_patch.jl:27
 [27] (::ADNLPModels.var"#lagrangian#55"{…})(x::Vector{…})
    @ ADNLPModels ~/.julia/packages/ADNLPModels/nR2jK/src/sparsity_pattern.jl:44
 [28] trace_function(::Type{…}, f::ADNLPModels.var"#lagrangian#55"{…}, x::Vector{…})
    @ SparseConnectivityTracer ~/.julia/packages/SparseConnectivityTracer/yg6Dn/src/trace_functions.jl:48
 [29] _hessian_sparsity(f::Function, x::Vector{Float64}, ::Type{SparseConnectivityTracer.HessianTracer{…}})
    @ SparseConnectivityTracer ~/.julia/packages/SparseConnectivityTracer/yg6Dn/src/trace_functions.jl:115
 [30] hessian_sparsity
    @ ~/.julia/packages/SparseConnectivityTracer/yg6Dn/src/adtypes_interface.jl:64 [inlined]
 [31] #compute_hessian_sparsity#54
    @ ~/.julia/packages/ADNLPModels/nR2jK/src/sparsity_pattern.jl:54 [inlined]
 [32] compute_hessian_sparsity
    @ ~/.julia/packages/ADNLPModels/nR2jK/src/sparsity_pattern.jl:35 [inlined]
 [33] macro expansion
    @ ~/.julia/packages/ADNLPModels/nR2jK/src/sparse_hessian.jl:31 [inlined]
 [34] macro expansion
    @ ./timing.jl:421 [inlined]
 [35] ADNLPModels.SparseADHessian(nvar::Int64, f::Function, ncon::Int64, c!::ADNLPModels.var"#2#3"; x0::Vector{…}, coloring_algorithm::SparseMatrixColorings.GreedyColoringAlgorithm{…}, detector::SparseConnectivityTracer.TracerSparsityDetector{…}, show_time::Bool, kwargs::@Kwargs{})
    @ ADNLPModels ~/.julia/packages/ADNLPModels/nR2jK/src/sparse_hessian.jl:30
 [36] macro expansion
    @ ~/.julia/packages/ADNLPModels/nR2jK/src/ad.jl:90 [inlined]
 [37] macro expansion
    @ ./timing.jl:421 [inlined]
 [38] ADNLPModels.ADModelBackend(nvar::Int64, f::typeof(misfit); backend::Symbol, matrix_free::Bool, show_time::Bool, gradient_backend::Type, hprod_backend::Type, hessian_backend::Type, kwargs::@Kwargs{…})
    @ ADNLPModels ~/.julia/packages/ADNLPModels/nR2jK/src/ad.jl:86
 [39] ADNLPModel(f::Function, x0::Vector{Float64}; name::String, minimize::Bool, kwargs::@Kwargs{})
    @ ADNLPModels ~/.julia/packages/ADNLPModels/nR2jK/src/nlp.jl:124
 [40] ADNLPModel(f::Function, x0::Vector{Float64})
    @ ADNLPModels ~/.julia/packages/ADNLPModels/nR2jK/src/nlp.jl:119
 [41] top-level scope
    @ ~/Desktop/PhD/Recherche/iR2N/Coding/debug/debug RegularizedProblems/test_patch.jl:37

And for ADNSLModel:

ERROR: TypeError: in typeassert, expected Float64, got a value of type SparseConnectivityTracer.HessianTracer{SparseConnectivityTracer.DictHessianPattern{Int64, BitSet, Dict{Int64, BitSet}, SparseConnectivityTracer.NotShared}}
Stacktrace:
  [1] setindex!(A::Vector{…}, x::SparseConnectivityTracer.HessianTracer{…}, i::Int64)
    @ Base ./array.jl:987
  [2] FH_ODE(dx::Vector{Float64}, x::Vector{Float64}, p::Vector{SparseConnectivityTracer.HessianTracer{…}}, t::Float64)
    @ Main ~/Desktop/PhD/Recherche/iR2N/Coding/debug/debug RegularizedProblems/test_patch.jl:7
  [3] (::SciMLBase.Void{typeof(FH_ODE)})(::Vector{Float64}, ::Vararg{Any})
    @ SciMLBase ~/.julia/packages/SciMLBase/Zrnkf/src/utils.jl:486
  [4] (::FunctionWrappers.CallWrapper{…})(f::SciMLBase.Void{…}, arg1::Vector{…}, arg2::Vector{…}, arg3::Vector{…}, arg4::Float64)
    @ FunctionWrappers ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:65
  [5] macro expansion
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:137 [inlined]
  [6] do_ccall
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:125 [inlined]
  [7] FunctionWrapper
    @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:144 [inlined]
  [8] _call
    @ ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:12 [inlined]
  [9] FunctionWrappersWrapper
    @ ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:10 [inlined]
 [10] ODEFunction
    @ ~/.julia/packages/SciMLBase/Zrnkf/src/scimlfunctions.jl:2470 [inlined]
 [11] ode_determine_initdt(u0::Vector{…}, t::Float64, tdir::Float64, dtmax::Float64, abstol::Float64, reltol::Float64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), prob::ODEProblem{…}, integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/initdt.jl:53
 [12] auto_dt_reset!
    @ ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/integrators/integrator_interface.jl:450 [inlined]
 [13] handle_dt!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/solve.jl:647
 [14] __init(prob::ODEProblem{…}, alg::Vern9{…}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{…}; saveat::Tuple{}, 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::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{…}, abstol::Nothing, reltol::Nothing, qmin::Rational{…}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{…}, 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), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias::ODEAliasSpecifier, initializealg::OrdinaryDiffEqCore.DefaultInit, kwargs::@Kwargs{})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/solve.jl:609
 [15] __init (repeats 5 times)
    @ ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/solve.jl:11 [inlined]
 [16] #__solve#62
    @ ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/solve.jl:6 [inlined]
 [17] __solve
    @ ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/solve.jl:1 [inlined]
 [18] #solve_call#35
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:635 [inlined]
 [19] solve_call
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:592 [inlined]
 [20] #solve_up#44
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1142 [inlined]
 [21] solve_up
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1120 [inlined]
 [22] #solve#42
    @ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1057 [inlined]
 [23] solve(prob::ODEProblem{…}, args::Vern9{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1047
 [24] simulate(p::Vector{SparseConnectivityTracer.HessianTracer{SparseConnectivityTracer.DictHessianPattern{…}}})
    @ Main ~/Desktop/PhD/Recherche/iR2N/Coding/debug/debug RegularizedProblems/test_patch.jl:13
 [25] residual(p::Vector{SparseConnectivityTracer.HessianTracer{SparseConnectivityTracer.DictHessianPattern{…}}})
    @ Main ~/Desktop/PhD/Recherche/iR2N/Coding/debug/debug RegularizedProblems/test_patch.jl:21
 [26] (::ADNLPModels.var"#F!#432"{…})(output::Vector{…}, x::Vector{…})
    @ ADNLPModels ~/.julia/packages/ADNLPModels/nR2jK/src/nls.jl:141
 [27] F
    @ ~/.julia/packages/ADNLPModels/nR2jK/src/ad.jl:233 [inlined]
 [28] (::ADNLPModels.var"#F#15"{Int64, ADNLPModels.var"#F#6#16"{…}})(x::Vector{SparseConnectivityTracer.HessianTracer{…}})
    @ ADNLPModels ~/.julia/packages/ADNLPModels/nR2jK/src/ad.jl:231
 [29] (::ADNLPModels.var"#7#17")(x::Vector{SparseConnectivityTracer.HessianTracer{…}})
    @ ADNLPModels ~/.julia/packages/ADNLPModels/nR2jK/src/ad.jl:236
 [30] (::ADNLPModels.var"#lagrangian#55"{…})(x::Vector{…})
    @ ADNLPModels ~/.julia/packages/ADNLPModels/nR2jK/src/sparsity_pattern.jl:44
 [31] trace_function(::Type{…}, f::ADNLPModels.var"#lagrangian#55"{…}, x::Vector{…})
    @ SparseConnectivityTracer ~/.julia/packages/SparseConnectivityTracer/yg6Dn/src/trace_functions.jl:48
 [32] _hessian_sparsity(f::Function, x::Vector{Float64}, ::Type{SparseConnectivityTracer.HessianTracer{…}})
    @ SparseConnectivityTracer ~/.julia/packages/SparseConnectivityTracer/yg6Dn/src/trace_functions.jl:115
 [33] hessian_sparsity
    @ ~/.julia/packages/SparseConnectivityTracer/yg6Dn/src/adtypes_interface.jl:64 [inlined]
 [34] #compute_hessian_sparsity#54
    @ ~/.julia/packages/ADNLPModels/nR2jK/src/sparsity_pattern.jl:54 [inlined]
 [35] compute_hessian_sparsity
    @ ~/.julia/packages/ADNLPModels/nR2jK/src/sparsity_pattern.jl:35 [inlined]
 [36] macro expansion
    @ ~/.julia/packages/ADNLPModels/nR2jK/src/sparse_hessian.jl:31 [inlined]
 [37] macro expansion
    @ ./timing.jl:421 [inlined]
 [38] ADNLPModels.SparseADHessian(nvar::Int64, f::Function, ncon::Int64, c!::ADNLPModels.var"#9#19"; x0::Vector{…}, coloring_algorithm::SparseMatrixColorings.GreedyColoringAlgorithm{…}, detector::SparseConnectivityTracer.TracerSparsityDetector{…}, show_time::Bool, kwargs::@Kwargs{})
    @ ADNLPModels ~/.julia/packages/ADNLPModels/nR2jK/src/sparse_hessian.jl:30
 [39] macro expansion
    @ ~/.julia/packages/ADNLPModels/nR2jK/src/ad.jl:266 [inlined]
 [40] macro expansion
    @ ./timing.jl:421 [inlined]
 [41] ADModelNLSBackend(nvar::Int64, F!::ADNLPModels.var"#F!#432"{…}, nequ::Int64; backend::Symbol, matrix_free::Bool, show_time::Bool, gradient_backend::Type, hprod_backend::Type, hessian_backend::Type, hprod_residual_backend::Type, jprod_residual_backend::Type, jtprod_residual_backend::Type, jacobian_residual_backend::Type, hessian_residual_backend::Type, kwargs::@Kwargs{…})
    @ ADNLPModels ~/.julia/packages/ADNLPModels/nR2jK/src/ad.jl:262
 [42] ADNLSModel!(F!::Function, x0::Vector{…}, nequ::Int64; linequ::Vector{…}, name::String, minimize::Bool, kwargs::@Kwargs{})
    @ ADNLPModels ~/.julia/packages/ADNLPModels/nR2jK/src/nls.jl:163
 [43] ADNLSModel!
    @ ~/.julia/packages/ADNLPModels/nR2jK/src/nls.jl:151 [inlined]
 [44] #ADNLSModel#431
    @ ~/.julia/packages/ADNLPModels/nR2jK/src/nls.jl:148 [inlined]
 [45] ADNLSModel(F::Function, x0::Vector{Float64}, nequ::Int64)
    @ ADNLPModels ~/.julia/packages/ADNLPModels/nR2jK/src/nls.jl:139
 [46] top-level scope
    @ ~/Desktop/PhD/Recherche/iR2N/Coding/debug/debug RegularizedProblems/test_patch.jl:40

[tmigot]

@nathanemac Do you use the Hessian and Jacobian or just matrix-vector products ? If so, you can pass the keyword matrix_free = true when instantiating the nlpmodels. This can be a workaround.

The issue is not related to ADNLPModels. The following clarify where the error comes from:

using ADTypes, SparseConnectivityTracer
detector = TracerSparsityDetector()
S = ADTypes.hessian_sparsity(misfit, x0, detector)

It looks like DifferentialEquations make some typing assumptions prevent the sparsity detection. @gdalle I have you seen this before?

[gdalle]

Yeah this is not an ADNLPModels issue, nor is it really a SparseConnectivityTracer issue. It's about which types the DiffEq problem expects. I'm not an expert but I see at least two solutions:

• also modify u0 when you remake the problem, using eltype(p), see eg https://docs.sciml.ai/DiffEqDocs/stable/basics/problem/#Modification-of-problem-types
• use PreallocationTools from the start so that the storage adapts to the required eltype, see eg Correctly initializes a SparseConnectivityTracer.GradientTracer cache SciML/PreallocationTools.jl#117

[nathanemac]

Thank you very much guys!
I'll investigate your options.