Skip to content

type promotion bug with parameter dependent initialization_eqs #3246

@SebastianM-C

Description

@SebastianM-C

Describe the bug 🐞

If a parameter changes its type (e.g. becomes a dual number) via setsym_oop and we need parameters in initialization_eqs, then the initialization problem does not correctly promote the type of prob.f.initialization_data.initializeprob.f.resid_prototype, leading to

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{:tag, Float64}, Float64, 1})

during initialization.

Expected behavior

The tests below should pass.

Minimal Reproducible Example 👇

@testset "type promotion with parameter dependent initialization_eqs" begin
    @variables x(t)=1 y(t)=1
    @parameters a=1
    @named sys = ODESystem([D(x) ~ 0, D(y)~x+a], t; initialization_eqs = [y ~ a])

    ssys = structural_simplify(sys)
    prob = ODEProblem(ssys, [], (0, 1), [])

    @test SciMLBase.successful_retcode(solve(prob))

    seta = setsym_oop(prob, [a])
    (newu0, newp) = seta(prob, ForwardDiff.Dual{ForwardDiff.Tag{:tag, Float64}}.([1.0], 1))
    newprob = remake(prob, u0=newu0, p=newp)

    @test SciMLBase.successful_retcode(solve(newprob))
end

Error & Stacktrace ⚠️

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{:tag, 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(::UInt8)
   @ Base float.jl:245
  ...

Stacktrace:
  [1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{:tag, Float64}, Float64, 1})
    @ Base ./number.jl:7
  [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{:tag, Float64}, Float64, 1}, i::Int64)
    @ Base ./array.jl:976
  [3] macro expansion
    @ ~/.julia/packages/SymbolicUtils/jf8aQ/src/code.jl:430 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/Symbolics/YbNrd/src/build_function.jl:546 [inlined]
  [5] macro expansion
    @ ~/.julia/packages/SymbolicUtils/jf8aQ/src/code.jl:387 [inlined]
  [6] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:163 [inlined]
  [7] macro expansion
    @ ./none:0 [inlined]
  [8] generated_callfunc
    @ ./none:0 [inlined]
  [9] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::Nothing, ::Vector{…})
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:150
 [10] (::ModelingToolkit.var"#f#749"{})(du::Vector{…}, u::Nothing, p::MTKParameters{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/0O7FS/src/systems/nonlinear/nonlinearsystem.jl:320
 [11] (::NonlinearFunction{…})(::Vector{…}, ::Vararg{…})
    @ SciMLBase ~/.julia/packages/SciMLBase/VAClc/src/scimlfunctions.jl:2359
 [12] build_null_solution(prob::NonlinearLeastSquaresProblem{…}, args::Nothing; abstol::Float64, kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:729
 [13] solve_call(_prob::NonlinearLeastSquaresProblem{…}, args::Nothing; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:624
 [14] solve_up(prob::NonlinearLeastSquaresProblem{…}, sensealg::Nothing, u0::Nothing, p::MTKParameters{…}, args::Nothing; kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1105
 [15] solve(prob::NonlinearLeastSquaresProblem{…}, args::Nothing; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1036
 [16] get_initial_values(prob::ODEProblem{…}, valp::OrdinaryDiffEqCore.ODEIntegrator{…}, f::Function, alg::SciMLBase.OverrideInit{…}, iip::Val{…}; nlsolve_alg::Nothing, abstol::Float64, reltol::Float64, kwargs::@Kwargs{})
    @ SciMLBase ~/.julia/packages/SciMLBase/VAClc/src/initialization.jl:209
 [17] _initialize_dae!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, prob::ODEProblem{…}, alg::SciMLBase.OverrideInit{…}, isinplace::Val{…})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/gALQU/src/initialize_dae.jl:159
 [18] _initialize_dae!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, prob::ODEProblem{…}, alg::OrdinaryDiffEqCore.DefaultInit, x::Val{…})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/gALQU/src/initialize_dae.jl:50
 [19] initialize_dae!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, initializealg::OrdinaryDiffEqCore.DefaultInit)
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/gALQU/src/initialize_dae.jl:40
 [20] __init(prob::ODEProblem{…}, alg::OrdinaryDiffEqCore.CompositeAlgorithm{…}, 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_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEqCore.DefaultInit, kwargs::@Kwargs{…})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/gALQU/src/solve.jl:514
 [21] __init (repeats 5 times)
    @ ~/.julia/packages/OrdinaryDiffEqCore/gALQU/src/solve.jl:11 [inlined]
 [22] #__solve#62
    @ ~/.julia/packages/OrdinaryDiffEqCore/gALQU/src/solve.jl:6 [inlined]
 [23] __solve
    @ ~/.julia/packages/OrdinaryDiffEqCore/gALQU/src/solve.jl:1 [inlined]
 [24] solve_call(_prob::ODEProblem{…}, args::OrdinaryDiffEqCore.CompositeAlgorithm{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:632
 [25] solve_call
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:589 [inlined]
 [26] #solve_up#53
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1120 [inlined]
 [27] solve_up
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1099 [inlined]
 [28] #solve#51
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1036 [inlined]
 [29] solve
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1026 [inlined]
 [30] #__solve#3
    @ ~/.julia/packages/OrdinaryDiffEqDefault/2qVWT/src/default_alg.jl:46 [inlined]
 [31] __solve
    @ ~/.julia/packages/OrdinaryDiffEqDefault/2qVWT/src/default_alg.jl:45 [inlined]
 [32] #__solve#72
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1436 [inlined]
 [33] __solve
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1427 [inlined]
 [34] #solve_call#44
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:632 [inlined]
 [35] solve_call(::ODEProblem{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:589
 [36] solve_up(::ODEProblem{…}, ::Nothing, ::Vector{…}, ::MTKParameters{…}; kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1105
 [37] solve_up
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1099 [inlined]
 [38] solve(::ODEProblem{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1036
 [39] solve(::ODEProblem{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1026
 [40] top-level scope
    @ REPL[98]:1
Some type information was truncated. Use `show(err)` to see complete types.

Environment (please complete the following information):

  • Output of using Pkg; Pkg.status()
Status `~/juliasim/dev/Project.toml`
  [2b5f629d] DiffEqBase v6.160.0
  [961ee093] ModelingToolkit v9.54.0
  • Output of using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)
  • Output of versioninfo()
Julia Version 1.11.1
Commit 8f5b7ca12ad (2024-10-16 10:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 32 × Intel(R) Core(TM) i9-14900K
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, alderlake)
Threads: 32 default, 0 interactive, 16 GC (on 32 virtual cores)

Additional context

Add any other context about the problem here.

Metadata

Metadata

Labels

bugSomething isn't working

Type

No type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions