Cannot initialize well-defined ODE with defaults

[hersle]

The well-defined ODE problem

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

@parameters y0
@variables x(t) y(t) z(t)
@named sys = ODESystem([
    D(x) ~ 0
    y ~ y0 / x
    D(z) ~ y
], t; defaults = [
    y0 => 1
    x => 1
    z => y
    #y => y0 / x # error message suggests this should be added, but it only results in an infinite loop
])
ssys = structural_simplify(sys)
prob = ODEProblem(ssys, [], (0.0, 1.0), [])

errors with

ERROR: LoadError: ArgumentError: Any[y(t)] are either missing from the variable map or missing from the system's unknowns/parameters list.
Stacktrace:
  [1] throw_missingvars_in_sys(vars::Vector{Any})
    @ ModelingToolkit C:\Users\herma\.julia\packages\ModelingToolkit\kAiI0\src\utils.jl:633
  [2] promote_to_concrete(vs::Vector{Any}; tofloat::Bool, use_union::Bool)
    @ ModelingToolkit C:\Users\herma\.julia\packages\ModelingToolkit\kAiI0\src\utils.jl:650
  [3] varmap_to_vars(varmap::Vector{…}, varlist::Vector{…}; defaults::Dict{…}, check::Bool, toterm::Function, promotetoconcrete::Nothing, tofloat::Bool, use_union::Bool)
    @ ModelingToolkit C:\Users\herma\.julia\packages\ModelingToolkit\kAiI0\src\variables.jl:165
  [4] get_u0(sys::ODESystem, u0map::Vector{Any}, parammap::SciMLBase.NullParameters; symbolic_u0::Bool, toterm::Function) 
    @ ModelingToolkit C:\Users\herma\.julia\packages\ModelingToolkit\kAiI0\src\systems\diffeqs\abstractodesystem.jl:736   
  [5] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{…}, parammap::Vector{…}; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, eval_module::Module, use_union::Bool, tofloat::Bool, symbolic_u0::Bool, u0_constructor::typeof(identity), guesses::Dict{…}, t::Float64, warn_initialize_determined::Bool, build_initializeprob::Bool, initialization_eqs::Vector{…}, kwargs::@Kwargs{…}) 
    @ ModelingToolkit C:\Users\herma\.julia\packages\ModelingToolkit\kAiI0\src\systems\diffeqs\abstractodesystem.jl:846   
  [6] process_DEProblem
    @ C:\Users\herma\.julia\packages\ModelingToolkit\kAiI0\src\systems\diffeqs\abstractodesystem.jl:741 [inlined]
  [7] (ODEProblem{…})(sys::ODESystem, u0map::Vector{…}, tspan::Tuple{…}, parammap::Vector{…}; callback::Nothing, check_length::Bool, warn_initialize_determined::Bool, eval_expression::Bool, eval_module::Module, kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\herma\.julia\packages\ModelingToolkit\kAiI0\src\systems\diffeqs\abstractodesystem.jl:998   
  [8] (ODEProblem{…})(sys::ODESystem, u0map::Vector{…}, tspan::Tuple{…}, parammap::Vector{…})
    @ ModelingToolkit C:\Users\herma\.julia\packages\ModelingToolkit\kAiI0\src\systems\diffeqs\abstractodesystem.jl:986   
  [9] (ODEProblem{true})(::ODESystem, ::Vector{Any}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\herma\.julia\packages\ModelingToolkit\kAiI0\src\systems\diffeqs\abstractodesystem.jl:973   
 [10] (ODEProblem{true})(::ODESystem, ::Vector{Any}, ::Vararg{Any})
    @ ModelingToolkit C:\Users\herma\.julia\packages\ModelingToolkit\kAiI0\src\systems\diffeqs\abstractodesystem.jl:972   
 [11] ODEProblem(::ODESystem, ::Vector{Any}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\herma\.julia\packages\ModelingToolkit\kAiI0\src\systems\diffeqs\abstractodesystem.jl:962   
 [12] ODEProblem(::ODESystem, ::Vector{Any}, ::Vararg{Any})  
    @ ModelingToolkit C:\Users\herma\.julia\packages\ModelingToolkit\kAiI0\src\systems\diffeqs\abstractodesystem.jl:961   
 [13] top-level scope
    @ C:\Users\herma\.julia\dev\Symboltz\bug.jl:17
in expression starting at C:\Users\herma\.julia\dev\Symboltz\bug.jl:17
Some type information was truncated. Use `show(err)` to see complete types.

[hersle]

The error suggests including a default for y, as well. If I do that by uncommenting the commented line, the program seems to enter an infinite loop. By intercepting the program with CTRL+C, it seems like it is stuck somewhere inside

  ...
  [7] fast_substitute(expr::SymbolicUtils.BasicSymbolic{Real}, subs::Dict{Any, Any}; operator::Type)
    @ Symbolics C:\Users\herma\.julia\packages\Symbolics\Nk48t\src\variable.jl:474
  [8] fixpoint_sub(x::SymbolicUtils.BasicSymbolic{Real}, dict::Dict{Any, Any}; operator::Type)
    @ Symbolics C:\Users\herma\.julia\packages\Symbolics\Nk48t\src\variable.jl:432
  [9] fixpoint_sub
    @ C:\Users\herma\.julia\packages\Symbolics\Nk48t\src\variable.jl:428 [inlined]
 [10] _varmap_to_vars(varmap::Dict{…}, varlist::Vector{…}; defaults::Dict{…}, check::Bool, toterm::Function, initialization_phase::Bool)
    @ ModelingToolkit C:\Users\herma\.julia\packages\ModelingToolkit\kAiI0\src\variables.jl:201
  ...

[hersle]

Debugging the original example (without the default for y) shows that get_u0() correctly finds the required defaults for the initialization:

However, as shown, the detected default for y is flipped from y => y0 / x to y0 / x => y. I think this is the culprit.
This flipping seems to be caused by this line. Why is this flipping done?

[ChrisRackauckas]

I think this fix should become a part of #2850. Or expanding observed in get_u0. Not too tricky to fix but indeed it looks like an edge case we missed.

[ChrisRackauckas]

I would investigate whether that flipping is still necessary. I'm not sure it is.

[hersle]

It seems to have been introduced here.

[ChrisRackauckas]

If you have x => p; y => p, then sending p => y everywhere makes the initialization problem pick up the implicit x => y.

[hersle]

Wouldn't such a problem initialize fine without using the implicit x => y? If x => p; y => p and p must be specified, then x and y are already known.

[ChrisRackauckas]

It wouldn't tear. But I don't see how this would be involved here since the flip only occurs on parameters.

[hersle]

Sorry, you're right. I misread the flipping change. I thought the commit changed the standard behavior from eq.rhs => eq.lhs to eq.lhs => eq.rhs, but the standard behavior has been eq.rhs => eq.lhs all along.

I guess the real question is when defaults from observed equations should be eq.rhs => eq.lhs (standard), and when they should be eq.lhs => eq.rhs (would fix my example).