Cannot solve ODE with Duals if it has initialization and callback

[hersle]

Here is a simple ODE that stops when y ~ 0.5:

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

@variables x(t) y(t)
stop!(integrator, _, _, _) = terminate!(integrator)
@named sys = ODESystem([
    D(x) ~ 1.0
    D(y) ~ 1.0
], t; initialization_eqs = [
    y ~ 0.0
], continuous_events = [
    [y ~ 0.5] => (stop!, [y], [], [], nothing)
])
sys = structural_simplify(sys)
prob0 = ODEProblem(sys, [x => NaN], (0.0, 1.0), [])

# final_x(x0) is equivalent to x0 + 0.5
function final_x(x0)
    prob = remake(prob0; u0 = [x => x0])
    sol = solve(prob)
    return sol[x][end]
end
final_x(0.3) # should be 0.8
ForwardDiff.derivative(final_x, 0.3) # should be 1.0

On #master, the last line fails with

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(final_x), 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{typeof(final_x), Float64}, Float64, 1})
    @ Base .\number.jl:7
  [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(final_x), Float64}, 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:373
  [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
    @ .\array.jl:791 [inlined]
 [11] map
    @ .\abstractarray.jl:3495 [inlined]
 [12] MultipleSetters
    @ C:\Users\herma\.julia\packages\SymbolicIndexingInterface\MHJDc\src\parameter_indexing.jl:695 [inlined]
 [13] (::ModelingToolkit.UpdateInitializeprob{…})(initializeprob::NonlinearProblem{…}, prob::ODEProblem{…})
    @ ModelingToolkit C:\Users\herma\.julia\dev\ModelingToolkit\src\systems\problem_utils.jl:458
 [14] get_initial_values(prob::ODEProblem{…}, valp::ODEProblem{…}, f::Function, alg::SciMLBase.OverrideInit{…}, iip::Val{…}; nlsolve_alg::Nothing, abstol::Nothing, reltol::Nothing, kwargs::@Kwargs{})
    @ SciMLBase C:\Users\herma\.julia\packages\SciMLBase\XzPx0\src\initialization.jl:234
 [15] get_initial_values(prob::ODEProblem{…}, valp::ODEProblem{…}, f::Function, alg::SciMLBase.OverrideInit{…}, iip::Val{…})
    @ SciMLBase C:\Users\herma\.julia\packages\SciMLBase\XzPx0\src\initialization.jl:221
 [16] remake(prob::ODEProblem{…}; f::Function, u0::Vector{…}, tspan::Tuple{…}, p::MTKParameters{…}, 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:263
 [17] get_concrete_problem(prob::ODEProblem{…}, isadapt::Bool; kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\herma\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:1223
 [18] get_concrete_problem
    @ C:\Users\herma\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:1209 [inlined]
 [19] solve_up(::ODEProblem{…}, ::Nothing, ::Vector{…}, ::MTKParameters{…}; kwargs::@Kwargs{})
    @ DiffEqBase C:\Users\herma\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:1105
 [20] solve_up
    @ C:\Users\herma\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:1101 [inlined]
 [21] solve(::ODEProblem{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
    @ DiffEqBase C:\Users\herma\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:1038
 [22] solve(::ODEProblem{…})
    @ DiffEqBase C:\Users\herma\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:1028
 [23] final_x(x0::ForwardDiff.Dual{ForwardDiff.Tag{typeof(final_x), Float64}, Float64, 1})
    @ Main c:\Users\herma\.julia\dev\SymBoltz\bug.jl:94
 [24] derivative(f::typeof(final_x), x::Float64)
    @ ForwardDiff C:\Users\herma\.julia\packages\ForwardDiff\UBbGT\src\derivative.jl:14

However, the code runs if

• initialization_eqs = [y ~ 0.0] is replaced with defaults = [y => 0.0], or
• continuous_events = [...] is removed altogether (with different results, of course).

[SebastianM-C]

This looks similar to #3318