Remake does not handle dual type promotion with dummy derivatives

[hersle]

Here is a simple throw-ball-straight-up-into-the-air example (2nd order ODE)

using ModelingToolkit, OrdinaryDiffEq, NonlinearSolve
using ModelingToolkit: t_nounits as t, D_nounits as D

# Throw ball straight up into the air
@variables y(t)
@named sys = ODESystem([
    D(D(y)) ~ -9.81
], t; initialization_eqs = [
    y^2 ~ 0 # initialize y = 0 in a way that builds an initialization problem
])
sys = structural_simplify(sys)

# Find initial throw velocity that reaches exactly 100 m after 1 s
dprob0 = ODEProblem(sys, [D(y) => NaN], (0.0, 1.0), []; guesses = [y => 0.0])
nprob = NonlinearProblem((ics, _) -> begin
    dprob = remake(dprob0, u0 = Dict(D(y) => ics[1]))
    dsol = solve(dprob, Tsit5())
    return [dsol[y][end] - 100.0]
end, [1.0])
nsol = solve(nprob, NewtonRaphson())

This crashes with

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{DifferentiationInterface.FixTail{…}, Float64}, Float64, 1})       
The type `Float64` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:265
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:900
  Float64(::Float16)
   @ Base float.jl:342
  ...

Stacktrace:
  [1] #1117
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1386 [inlined]
  [2] (::ModelingToolkit.var"#1117#1120")(::Pair{Any, Any})
    @ ModelingToolkit .\none:0
  [3] iterate
    @ .\generator.jl:48 [inlined]
  [4] grow_to!(dest::Dict{Any, Any}, itr::Base.Generator{Dict{Any, Any}, ModelingToolkit.var"#1117#1120"})
    @ Base .\abstractdict.jl:587
  [5] dict_with_eltype
    @ .\abstractdict.jl:641 [inlined]
  [6] Dict
    @ .\dict.jl:117 [inlined]
  [7] ModelingToolkit.InitializationProblem{…}(sys::ODESystem, t::Float64, u0map::Dict{…}, parammap::Dict{…}; guesses::Dict{…}, check_length::Bool, warn_initialize_determined::Bool, initialization_eqs::Vector{…}, fully_determined::Nothing, check_units::Bool, use_scc::Bool, kwargs::@Kwargs{})    
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1385
  [8] InitializationProblem
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1287 [inlined]
  [9] #_#1112
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1265 [inlined]
 [10] InitializationProblem
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1264 [inlined]
 [11] #InitializationProblem#1110
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1253 [inlined]
 [12] InitializationProblem
    @ C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:1252 [inlined]
 [13] maybe_build_initialization_problem(sys::ODESystem, op::Dict{…}, u0map::Dict{…}, pmap::Dict{…}, t::Float64, defs::Dict{…}, guesses::Dict{…}, missing_unknowns::Set{…}; implicit_dae::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\problem_utils.jl:546
 [14] remake_initialization_data(sys::ODESystem, odefn::Function, u0::Dict{…}, t0::Float64, p::Missing, newu0::Vector{…}, newp::MTKParameters{…})   
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\nonlinear\initializesystem.jl:343
 [15] remake(prob::ODEProblem{…}; f::Missing, u0::Dict{…}, tspan::Missing, p::Missing, kwargs::Missing, interpret_symbolicmap::Bool, build_initializeprob::Bool, use_defaults::Bool, lazy_initialization::Nothing, _kwargs::@Kwargs{})
    @ SciMLBase C:\Users\herma\.julia\packages\SciMLBase\XzPx0\src\remake.jl:230
  ...

[hersle]

I think the issue is here:

• fullmap contains Differential(t)(y(t)) (without dummy derivative)
• unknowns(isys) and/or observed(isys) contains y_t(t) (with dummy derivative)
• These two are not found to match, so u0T does not promote to the type of D(y)

[hersle]

EDIT: skip this, simpler example below

The PR above (so far) fixes the original example. But here is one that still crashes:

using ModelingToolkit, OrdinaryDiffEq, NonlinearSolve
using ModelingToolkit: t_nounits as t, D_nounits as D

@variables a(t) E(t) ρ(t) ρr(t) ρm(t) ρQ(t) ϕ(t)
@parameters ρr0 ρm0
@named sys = ODESystem([
    ρr ~ ρr0 / a^4
    ρm ~ ρm0 / a^3
    ρQ ~ (D(ϕ) / a)^2 / 2
    ρ ~ ρr + ρm + ρQ
    E ~ D(a) / a
    D(a) ~ √(8π/3 * ρ * a^4)
    D(D(ϕ)) ~ -2*E*D(ϕ)
], t; initialization_eqs = [
    a ~ D(a) * t
])
sys = structural_simplify(sys)

aidx = ModelingToolkit.variable_index(sys, a)
callback = ContinuousCallback((u, _, _) -> (val = u[aidx]; val - 1.0), terminate!) # stop when a = 1
dprob0 = ODEProblem(sys, [ϕ => 1.0, D(ϕ) => NaN], (1e-5, 4.0), [ρr0 => 1e-5 * 3/8π, ρm0 => 0.3 * 3/8π]; guesses = [a => 0.1], callback)

function f(ics, _)
    dprob = remake(dprob0; u0 = Dict(D(ϕ) => ics[1]))
    dsol = solve(dprob, Tsit5())
    return [dsol[E][end] - 1.0]
end
nprob = NonlinearProblem(f, [1.0])
nsol = solve(nprob, TrustRegion())

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{DifferentiationInterface.FixTail{…}, Float64}, Float64, 1})
The type `Float64` exists, but no method is defined for this combination of argument types when trying to construct it. 

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:265
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:900
  Float64(::Float16)
   @ Base float.jl:342
  ...

Stacktrace:
  [1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{DifferentiationInterface.FixTail{…}, Float64}, Float64, 1})
    @ Base .\number.jl:7
  [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{…}, Float64, 1}, i::Int64)
    @ Base .\array.jl:987
  [3] set_parameter!(p::MTKParameters{…}, val::ForwardDiff.Dual{…}, pidx::ModelingToolkit.ParameterIndex{…})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\parameter_buffer.jl:374
  [4] set_parameter!(sys::NonlinearProblem{…}, val::ForwardDiff.Dual{…}, idx::ModelingToolkit.ParameterIndex{…})        
    @ SymbolicIndexingInterface C:\Users\herma\.julia\packages\SymbolicIndexingInterface\MHJDc\src\value_provider_interface.jl:67
  [5] (::SymbolicIndexingInterface.SetParameterIndex{…})(prob::NonlinearProblem{…}, val::ForwardDiff.Dual{…})
    @ SymbolicIndexingInterface C:\Users\herma\.julia\packages\SymbolicIndexingInterface\MHJDc\src\parameter_indexing.jl:678
  [6] (::SymbolicIndexingInterface.ParameterHookWrapper{…})(prob::NonlinearProblem{…}, args::ForwardDiff.Dual{…})       
    @ SymbolicIndexingInterface C:\Users\herma\.julia\packages\SymbolicIndexingInterface\MHJDc\src\parameter_indexing.jl:646
  [7] (::SymbolicIndexingInterface.var"#44#45"{…})(s!::SymbolicIndexingInterface.ParameterHookWrapper{…}, v::ForwardDiff.Dual{…})
    @ SymbolicIndexingInterface C:\Users\herma\.julia\packages\SymbolicIndexingInterface\MHJDc\src\parameter_indexing.jl:695
  [8] (::Base.var"#4#5"{…})(a::Tuple{…})
    @ Base .\generator.jl:37
  [9] iterate
    @ .\generator.jl:48 [inlined]
 [10] collect_to!(dest::Vector{…}, itr::Base.Generator{…}, offs::Int64, st::Tuple{…})
    @ Base .\array.jl:849
 [11] collect_to_with_first!
    @ .\array.jl:827 [inlined]
 [12] collect(itr::Base.Generator{Base.Iterators.Zip{Tuple{…}}, Base.var"#4#5"{SymbolicIndexingInterface.var"#44#45"{…}}})
    @ Base .\array.jl:801
 [13] map
    @ .\abstractarray.jl:3495 [inlined]
 [14] (::SymbolicIndexingInterface.MultipleSetters{…})(prob::NonlinearProblem{…}, val::Vector{…})
    @ SymbolicIndexingInterface C:\Users\herma\.julia\packages\SymbolicIndexingInterface\MHJDc\src\parameter_indexing.jl:695
 [15] (::ModelingToolkit.UpdateInitializeprob{…})(initializeprob::NonlinearProblem{…}, prob::OrdinaryDiffEqCore.ODEIntegrator{…})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\problem_utils.jl:459
 [16] get_initial_values(prob::ODEProblem{…}, valp::OrdinaryDiffEqCore.ODEIntegrator{…}, f::Function, alg::SciMLBase.OverrideInit{…}, iip::Val{…}; nlsolve_alg::NonlinearSolvePolyAlgorithm{…}, abstol::Float64, reltol::Float64, kwargs::@Kwargs{})
    @ SciMLBase C:\Users\herma\.julia\packages\SciMLBase\XzPx0\src\initialization.jl:234
 [17] get_initial_values
    @ C:\Users\herma\.julia\packages\SciMLBase\XzPx0\src\initialization.jl:221 [inlined]
 [18] _initialize_dae!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, prob::ODEProblem{…}, alg::SciMLBase.OverrideInit{…}, isinplace::Val{…})
    @ OrdinaryDiffEqCore C:\Users\herma\.julia\packages\OrdinaryDiffEqCore\3Talm\src\initialize_dae.jl:160
 [19] _initialize_dae!
    @ C:\Users\herma\.julia\packages\OrdinaryDiffEqCore\3Talm\src\initialize_dae.jl:50 [inlined]
 [20] initialize_dae! (repeats 2 times)
    @ C:\Users\herma\.julia\packages\OrdinaryDiffEqCore\3Talm\src\initialize_dae.jl:40 [inlined]
 [21] __init(prob::ODEProblem{…}, alg::Tsit5{…}, 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::ContinuousCallback{…}, dense::Bool, calck::Bool, dt::ForwardDiff.Dual{…}, dtmin::ForwardDiff.Dual{…}, dtmax::ForwardDiff.Dual{…}, 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 C:\Users\herma\.julia\packages\OrdinaryDiffEqCore\3Talm\src\solve.jl:560
 [22] __init (repeats 5 times)
    @ C:\Users\herma\.julia\packages\OrdinaryDiffEqCore\3Talm\src\solve.jl:11 [inlined]
 [23] #__solve#62
    @ C:\Users\herma\.julia\packages\OrdinaryDiffEqCore\3Talm\src\solve.jl:6 [inlined]
 [24] __solve
    @ C:\Users\herma\.julia\packages\OrdinaryDiffEqCore\3Talm\src\solve.jl:1 [inlined]
 [25] solve_call(_prob::ODEProblem{…}, args::Tsit5{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{}) 
    @ DiffEqBase C:\Users\herma\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:634
 [26] solve_call(_prob::ODEProblem{…}, args::Tsit5{…})      
    @ DiffEqBase C:\Users\herma\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:591
 [27] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::MTKParameters{…}, args::Tsit5{…}; kwargs::@Kwargs{})
    @ DiffEqBase C:\Users\herma\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:1122
 [28] solve_up
    @ C:\Users\herma\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:1101 [inlined]
 [29] solve(prob::ODEProblem{…}, args::Tsit5{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
    @ DiffEqBase C:\Users\herma\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:1038
 [30] solve(prob::ODEProblem{…}, args::Tsit5{…})
    @ DiffEqBase C:\Users\herma\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:1028
 [31] f(ics::Vector{ForwardDiff.Dual{ForwardDiff.Tag{…}, Float64, 1}}, ::SciMLBase.NullParameters)
    @ Main c:\Users\herma\.julia\dev\SymBoltz\bug.jl:25     
  ...

Debugging indicates the crash happens when it tries to set t in the initialization problem to Dual(1.0e-5,0.0).

[hersle]

Looks like it is somehow related to using a callback together with autodiff. Simpler example that crashes similarly to above:

using ModelingToolkit, OrdinaryDiffEq, NonlinearSolve, ForwardDiff
using ModelingToolkit: t_nounits as t, D_nounits as D

# Throw ball straight up into the air
@variables y(t) v(t)
@named sys = ODESystem([
    D(D(y)) ~ -9.81 # acceleration
    v ~ D(y) # velocity
], t)
sys = structural_simplify(sys)

yidx = ModelingToolkit.variable_index(sys, y)
callback = ContinuousCallback((u, _, _) -> (u[yidx] - 10.0), terminate!) # stop when y = 10 m
dprob0 = ODEProblem(sys, [y => 0.0, D(y) => NaN], (0.0, 50.0), []; callback)

function final_height(v0)
    dprob = remake(dprob0; u0 = Dict(D(y) => v0))
    dsol = solve(dprob, Tsit5())
    return dsol[y][end]
end
final_height(50.0)
ForwardDiff.derivative(final_height, 50.0) # should be 0

• Replacing D(y) => v0 by v => v0 makes the code run with the expected result, because it skips the initialization problem.
• Removing ; callback from ODEProblem makes the code run (though with a different result, which is OK).

[hersle]

The above PR seems to fix this particular example.

I will open another issue for my other issue, which seems unrelated to dummy derivatives upon further investigation.

[AayushSabharwal]

This has already been fixed and tested