`remake` doesn't evaluate interdependencies between `u0` and `p`

[jonathanfischer97]

Made a comment about this in Catalyst.jl:

"Was just about to make an issue related to this, it's a particular problem when using remake to update the initial conditions of an ODEProblem(rs, ...; remove_conserved = true).

• This is because the Γ conservation law constants are defined as parameters, which have symbolic dependencies on the initial conditions.
• Perhaps an overload is necessary for the internal SciMLBase.remake function here:

function _updated_u0_p_symmap(prob, u0, ::Val{true}, p, ::Val{true})
    isu0dep = any(symbolic_type(v) !== NotSymbolic() for (_, v) in u0)
    ispdep = any(symbolic_type(v) !== NotSymbolic() for (_, v) in p)

    if !isu0dep && !ispdep
        return remake_buffer(prob, state_values(prob), u0),
        remake_buffer(prob, parameter_values(prob), p)
    end
    if !isu0dep
        u0 = remake_buffer(prob, state_values(prob), u0)
        return _updated_u0_p_symmap(prob, u0, Val(false), p, Val(true))
    end
    if !ispdep
        p = remake_buffer(prob, parameter_values(prob), p)
        return _updated_u0_p_symmap(prob, u0, Val(true), p, Val(false))
    end

    varmap = merge(u0, p)
    u0 = anydict(k => symbolic_type(v) === NotSymbolic() ? v : symbolic_evaluate(v, varmap)
    for (k, v) in u0)
    p = anydict(k => symbolic_type(v) === NotSymbolic() ? v : symbolic_evaluate(v, varmap)
    for (k, v) in p)
    return remake_buffer(prob, state_values(prob), u0),
    remake_buffer(prob, parameter_values(prob), p)
end

• Here the main issue I've found are that the ispdep check doesn't find that Γ in the parameters is a symbolic type, because it is evaluated to numerical values in the ODEProblem. The only place where the symbolic expressions for Γ are exposed are in defaults field of the ODESystem.
• This obstructs the ostensible allowance for interdependencies between u0 and p, seen with the varmap = merge(u0, p) line.
• Currently my workaround is calling symbolic_evaluate prior to calling remake using the symbolic expressions from the defaults map:

function compute_Γ(tunable_var_value_map, constraint_syms, sys_defs::T) where T

 #* Merge the default system definitions with the new symbolic map
 merged_defs = merge(sys_defs, Dict(tunable_var_value_map))

 NTuple{4, eltype(tunable_var_value_map)}(constraint => symbolic_evaluate(merged_defs[constraint], merged_defs) for constraint in constraint_syms)
end

I can make a separate issue for this in this repo (not sure if it pertains most here, or ModelingToolkit.jl or SciMLBase.jl) if this is truly a problem and I'm not missing anything."

Originally posted by @jonathanfischer97 in https://github.com/SciML/Catalyst.jl/issues/806#issuecomment-2059826318

I should mention that I'm using ModelingToolkit.jl version 8.75 rather than 9.0, due to compat with Catalyst.jl, if that makes a difference.

[AayushSabharwal]

Could you give an MWE? We have tests specifically for interdependent u0 and p so this should work. If you're relying on it using defaults, you'll need to pass use_defaults = true

[SebastianM-C]

I'm not sure if this would work on MTK@v8.

[ChrisRackauckas]

Yes it's only v9

[AayushSabharwal]

Interdependent initialization should work even on v8, but thanks for pointing out that it won't respect defaults on MTKv8

[jonathanfischer97]

Yes, I initially tried use_defaults = true, but got an error related to symbolic_container.

[AayushSabharwal]

Try updating your SymbolicIndexingInterface version.

It won't start using defaults for MTKv8, but it'll at least not error the same way 😅

[hersle]

Here is a simple example that fails, which I think is related to this issue:

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

@parameters P
@variables x(t) y(t)
sys = complete(ODESystem([D(x) ~ P, D(y) ~ P], t, [x, y], [P]; defaults = [x => P + 1.0], name = :sys))
prob1 = ODEProblem(sys, unknowns(sys) .=> NaN, (0.0, 1.0), parameters(sys) .=> NaN)
prob2 = remake(prob1; u0 = [y => 2.0], p = [P => 1.0], use_defaults = true) # fails

It errors with

ERROR: MethodError: no method matching getindex(::Nothing, ::Int64)
Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/SymbolicUtils/qyMYa/src/code.jl:375 [inlined]
  [2] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:163 [inlined]
  [3] macro expansion
    @ ./none:0 [inlined]
  [4] generated_callfunc
    @ ./none:0 [inlined]
  [5] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Nothing, ::Vector{…}, ::Nothing)
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:150
  [6] (::ModelingToolkit.var"#709#722"{…})(u::Nothing, p::ModelingToolkit.MTKParameters{…}, t::Nothing)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/diffeqs/abstractodesystem.jl:464
  [7] (::SymbolicIndexingInterface.TimeDependentObservedFunction{…})(::SymbolicIndexingInterface.NotTimeseries, prob::SymbolicIndexingInterface.ProblemState{…})
    @ SymbolicIndexingInterface ~/.julia/packages/SymbolicIndexingInterface/SzZz3/src/state_indexing.jl:88
  [8] (::SymbolicIndexingInterface.TimeDependentObservedFunction{…})(prob::SymbolicIndexingInterface.ProblemState{…})
    @ SymbolicIndexingInterface ~/.julia/packages/SymbolicIndexingInterface/SzZz3/src/value_provider_interface.jl:144
  [9] (::SciMLBase.var"#667#671"{ODEProblem{…}, SymbolicIndexingInterface.ProblemState{…}})(::Pair{Any, Any})
    @ SciMLBase ./none:0
 [10] iterate
    @ ./generator.jl:47 [inlined]
 [11] Dict{Any, Any}(kv::Base.Generator{Dict{…}, SciMLBase.var"#667#671"{…}})
    @ Base ./dict.jl:83
 [12] anydict(d::Base.Generator{Dict{…}, SciMLBase.var"#667#671"{…}})
    @ SciMLBase ~/.julia/packages/SciMLBase/SDjaO/src/remake.jl:385
 [13] _updated_u0_p_symmap(prob::ODEProblem{…}, u0::Dict{…}, ::Val{…}, p::ModelingToolkit.MTKParameters{…}, ::Val{…})
    @ SciMLBase ~/.julia/packages/SciMLBase/SDjaO/src/remake.jl:479
 [14] _updated_u0_p_symmap(prob::ODEProblem{…}, u0::Dict{…}, ::Val{…}, p::Dict{…}, ::Val{…})
    @ SciMLBase ~/.julia/packages/SciMLBase/SDjaO/src/remake.jl:517
 [15] _updated_u0_p_internal(prob::ODEProblem{…}, u0::Vector{…}, p::Vector{…}; interpret_symbolicmap::Bool, use_defaults::Bool)
    @ SciMLBase ~/.julia/packages/SciMLBase/SDjaO/src/remake.jl:427
 [16] _updated_u0_p_internal
    @ ~/.julia/packages/SciMLBase/SDjaO/src/remake.jl:412 [inlined]
 [17] #updated_u0_p#688
    @ ~/.julia/packages/SciMLBase/SDjaO/src/remake.jl:548 [inlined]
 [18] updated_u0_p
    @ ~/.julia/packages/SciMLBase/SDjaO/src/remake.jl:529 [inlined]
 [19] remake(prob::ODEProblem{…}; f::Missing, u0::Vector{…}, tspan::Missing, p::Vector{…}, kwargs::Missing, interpret_symbolicmap::Bool, use_defaults::Bool, _kwargs::@Kwargs{})
    @ SciMLBase ~/.julia/packages/SciMLBase/SDjaO/src/remake.jl:95

It seems to me like _updated_u0_p_symmap(prob, u0, ::Val{true}, p, ::Val{false}) fails to substitute in P => 1.0 when evaluating x => P + 1.0.

This is with the latest

  [0c46a032] DifferentialEquations v7.13.0
  [961ee093] ModelingToolkit v9.15.0
  [0bca4576] SciMLBase v2.38.0
  [2efcf032] SymbolicIndexingInterface v0.3.21