feat: add circular dependency checking, substitution limit

Author: AayushSabharwal

src/systems/problem_utils.jl

@@ -250,6 +250,21 @@ function add_parameter_dependencies!(sys::AbstractSystem, varmap::AbstractDict)
     add_observed_equations!(varmap, parameter_dependencies(sys))
 end
 
+struct UnexpectedSymbolicValueInVarmap <: Exception
+    sym::Any
+    val::Any
+end
+
+function Base.showerror(io::IO, err::UnexpectedSymbolicValueInVarmap)
+    println(io,
+        """
+        Found symbolic value $(err.val) for variable $(err.sym). You may be missing an \
+        initial condition or have cyclic initial conditions. If this is intended, pass \
+        `symbolic_u0 = true`. In case the initial conditions are not cyclic but \
+        require more substitutions to resolve, increase `substitution_limit`.
+        """)
+end
+
 """
     $(TYPEDSIGNATURES)
 
@@ -278,6 +293,12 @@ function better_varmap_to_vars(varmap::AbstractDict, vars::Vector;
         isempty(missing_vars) || throw(MissingVariablesError(missing_vars))
     end
     vals = map(x -> varmap[x], vars)
+    if !allow_symbolic
+        for (sym, val) in zip(vars, vals)
+            symbolic_type(val) == NotSymbolic() && continue
+            throw(UnexpectedSymbolicValueInVarmap(sym, val))
+        end
+    end
 
     if container_type <: Union{AbstractDict, Tuple, Nothing, SciMLBase.NullParameters}
         container_type = Array
@@ -298,17 +319,83 @@ function better_varmap_to_vars(varmap::AbstractDict, vars::Vector;
     end
 end
 
+struct SubstitutionCycleError <: Exception
+    required_variables::Vector{Any}
+    cycle::Vector{Any}
+    rules::Vector{Any}
+end
+
+function Base.showerror(io::IO, err::SubstitutionCycleError)
+    println(io, "Detected a cyclic initial condition!")
+    println(io, "Variables being initialized: ", err.required_variables)
+    println(io)
+    println(io, "Variables in cycle: ", err.cycle)
+    println(io)
+    println(io, "Substitution rules:")
+    for rule in err.rules
+        println(io, rule)
+    end
+    println(io,
+        "Pass `check_cyclic_dependency = false` to disable this check if you think this is an error.")
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Check if any of the substitution rules in `varmap` lead to cycles involving
+variables in `vars`. Throw an error if such a cycle is found.
+"""
+function check_substitution_cycles(varmap::AbstractDict, vars)
+    # ordered set so that `vars` are the first `k` in the list
+    allvars = OrderedSet(vars)
+    union!(allvars, keys(varmap))
+    allvars = collect(allvars)
+    var_to_idx = Dict(allvars .=> eachindex(allvars))
+    graph = SimpleDiGraph(length(allvars))
+
+    buffer = Set()
+    for (k, v) in varmap
+        kidx = var_to_idx[k]
+        if symbolic_type(v) != NotSymbolic()
+            vars!(buffer, v)
+            for var in buffer
+                haskey(var_to_idx, var) || continue
+                add_edge!(graph, kidx, var_to_idx[var])
+            end
+        elseif v isa AbstractArray
+            for val in v
+                vars!(buffer, val)
+            end
+            for var in buffer
+                haskey(var_to_idx, var) || continue
+                add_edge!(graph, kidx, var_to_idx[var])
+            end
+        end
+        empty!(buffer)
+    end
+
+    # detect at most 100 cycles involving at most `length(varmap)` vertices
+    cycles = Graphs.simplecycles_limited_length(graph, length(varmap), 100)
+    # only count those which contain variables in `vars`
+    filter!(Base.Fix1(any, <=(length(vars))), cycles)
+
+    if !isempty(cycles)
+        cyclevars = allvars[cycles[1]]
+        throw(SubstitutionCycleError(vars, cyclevars, [v => varmap[v] for v in cyclevars]))
+    end
+end
+
 """
     $(TYPEDSIGNATURES)
 
 Performs symbolic substitution on the values in `varmap` for the keys in `vars`, using
 `varmap` itself as the set of substitution rules. If an entry in `vars` is not a key
 in `varmap`, it is ignored.
 """
-function evaluate_varmap!(varmap::AbstractDict, vars)
+function evaluate_varmap!(varmap::AbstractDict, vars; limit = 100)
     for k in vars
         haskey(varmap, k) || continue
-        varmap[k] = fixpoint_sub(varmap[k], varmap)
+        varmap[k] = fixpoint_sub(varmap[k], varmap; maxiters = limit)
     end
 end
 
@@ -407,6 +494,10 @@ Keyword arguments:
   length of `u0` vector for consistency. If `false`, do not check with equations. This is
   forwarded to `check_eqs_u0`
 - `symbolic_u0` allows the returned `u0` to be an array of symbolics.
+- `check_cyclic_dependency`: Whether to check the operating point for cyclic
+  substitution rules involving unknowns/parameters.
+- `substitution_limit`: The number times to substitute initial conditions into each
+  other to attempt to arrive at a numeric value.
 
 All other keyword arguments are passed as-is to `constructor`.
 """
@@ -416,7 +507,9 @@ function process_SciMLProblem(
         warn_initialize_determined = true, initialization_eqs = [],
         eval_expression = false, eval_module = @__MODULE__, fully_determined = false,
         check_initialization_units = false, tofloat = true, use_union = false,
-        u0_constructor = identity, du0map = nothing, check_length = true, symbolic_u0 = false, kwargs...)
+        u0_constructor = identity, du0map = nothing, check_length = true,
+        symbolic_u0 = false, check_cyclic_dependency = true,
+        substitution_limit = length(observed(sys)), kwargs...)
     dvs = unknowns(sys)
     ps = parameters(sys)
     iv = has_iv(sys) ? get_iv(sys) : nothing
@@ -466,7 +559,7 @@ function process_SciMLProblem(
             initializeprob = ModelingToolkit.InitializationProblem(
                 sys, t, u0map, pmap; guesses, warn_initialize_determined,
                 initialization_eqs, eval_expression, eval_module, fully_determined,
-                check_units = check_initialization_units)
+                check_units = check_initialization_units, check_cyclic_dependency)
             initializeprobmap = getu(initializeprob, unknowns(sys))
 
             punknowns = [p
@@ -503,7 +596,8 @@ function process_SciMLProblem(
     add_observed!(sys, op)
     add_parameter_dependencies!(sys, op)
 
-    evaluate_varmap!(op, dvs)
+    check_cyclic_dependency && check_substitution_cycles(op, dvs)
+    evaluate_varmap!(op, dvs; limit = substitution_limit)
 
     u0 = better_varmap_to_vars(
         op, dvs; tofloat = true, use_union = false,
@@ -515,7 +609,8 @@ function process_SciMLProblem(
 
     check_eqs_u0(eqs, dvs, u0; check_length, kwargs...)
 
-    evaluate_varmap!(op, ps)
+    check_cyclic_dependency && check_substitution_cycles(op, ps)
+    evaluate_varmap!(op, ps; limit = substitution_limit)
     if is_split(sys)
         p = MTKParameters(sys, op)
     else
@@ -527,6 +622,7 @@ function process_SciMLProblem(
         du0map = to_varmap(du0map, ddvs)
         merge!(op, du0map)
 
+        check_cyclic_dependency && check_substitution_cycles(op, ddvs)
         du0 = varmap_to_vars(du0map, ddvs; toterm = identity,
             tofloat = true)
         kwargs = merge(kwargs, (; ddvs))

test/initial_values.jl

@@ -149,3 +149,18 @@ end
     pmap = [c1 => 5.0, c2 => 1.0, c3 => 1.2]
     oprob = ODEProblem(osys, u0map, (0.0, 10.0), pmap)
 end
+
+@testset "Cyclic dependency checking and substitution limits" begin
+    @variables x(t) y(t)
+    @mtkbuild sys = ODESystem([D(x) ~ x, D(y) ~ y], t)
+    @test_throws ModelingToolkit.SubstitutionCycleError ODEProblem(
+        sys, [x => 2y, y => 3x], (0.0, 1.0); build_initializeprob = false)
+    @test_throws ModelingToolkit.UnexpectedSymbolicValueInVarmap ODEProblem(
+        sys, [x => 2y, y => 3x], (0.0, 1.0);
+        build_initializeprob = false, check_cyclic_dependency = false)
+
+    @parameters p q
+    @mtkbuild sys = ODESystem([D(x) ~ x * p, D(y) ~ y * q], t)
+    @test_throws ModelingToolkit.SubstitutionCycleError ODEProblem(
+        sys, [x => 1, y => 2], (0.0, 1.0), [p => 2q, q => 3p]; build_initializeprob = false)
+end