Merge branch 'master' into myb/state_priority

Author: YingboMa

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <[email protected]>", "Chris Rackauckas <[email protected]> and contributors"]
-version = "8.23.0"
+version = "8.24.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

src/bipartite_graph.jl

@@ -63,7 +63,17 @@ end
 function Base.setindex!(m::Matching{U}, v::Union{Integer, U}, i::Integer) where {U}
     if m.inv_match !== nothing
         oldv = m.match[i]
-        isa(oldv, Int) && (m.inv_match[oldv] = unassigned)
+        # TODO: maybe default Matching to always have an `inv_match`?
+
+        # To maintain the invariant that `m.inv_match[m.match[i]] == i`, we need
+        # to unassign the matching at `m.inv_match[v]` if it exists.
+        if v isa Int && (iv = m.inv_match[v]) isa Int
+            m.match[iv] = unassigned
+        end
+        if isa(oldv, Int)
+            @assert m.inv_match[oldv] == i
+            m.inv_match[oldv] = unassigned
+        end
         isa(v, Int) && (m.inv_match[v] = i)
     end
     return m.match[i] = v

src/inputoutput.jl

@@ -160,9 +160,9 @@ has_var(ex, x) = x ∈ Set(get_variables(ex))
 # Build control function
 
 """
-    (f_oop, f_ip), dvs, p = generate_control_function(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys); implicit_dae = false, ddvs = if implicit_dae
+    (f_oop, f_ip), dvs, p = generate_control_function(sys::AbstractODESystem, inputs = unbound_inputs(sys); implicit_dae = false, ddvs = if implicit_dae
 
-For a system `sys` that has unbound inputs (as determined by [`unbound_inputs`](@ref)), generate a function with additional input argument `in`
+For a system `sys` with inputs (as determined by [`unbound_inputs`](@ref) or user specified), generate a function with additional input argument `in`
 ```
 f_oop : (u,in,p,t)      -> rhs
 f_ip  : (uout,u,in,p,t) -> nothing
@@ -187,7 +187,7 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu
         error("No unbound inputs were found in system.")
     end
 
-    sys, diff_idxs, alge_idxs = io_preprocessing(sys, inputs, []; simplify, kwargs...)
+    sys, _ = io_preprocessing(sys, inputs, []; simplify, kwargs...)
 
     dvs = states(sys)
     ps = parameters(sys)
@@ -269,8 +269,8 @@ function inputs_to_parameters!(state::TransformationState, io)
         @assert new_v > 0
         new_var_to_diff[new_i] = new_v
     end
-    @set! structure.var_to_diff = new_var_to_diff
-    @set! structure.graph = new_graph
+    @set! structure.var_to_diff = complete(new_var_to_diff)
+    @set! structure.graph = complete(new_graph)
 
     @set! sys.eqs = map(Base.Fix2(substitute, input_to_parameters), equations(sys))
     @set! sys.states = setdiff(states(sys), keys(input_to_parameters))

src/structural_transformation/StructuralTransformations.jl

@@ -21,7 +21,8 @@ using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative, Di
                        ExtraVariablesSystemException,
                        get_postprocess_fbody, vars!,
                        IncrementalCycleTracker, add_edge_checked!, topological_sort,
-                       invalidate_cache!, Substitutions, get_or_construct_tearing_state
+                       invalidate_cache!, Substitutions, get_or_construct_tearing_state,
+                       AliasGraph
 
 using ModelingToolkit.BipartiteGraphs
 import .BipartiteGraphs: invview

src/structural_transformation/pantelides.jl

@@ -74,16 +74,47 @@ end
 
 Perform Pantelides algorithm.
 """
-function pantelides!(state::TransformationState; maxiters = 8000)
+function pantelides!(state::TransformationState, ag::Union{AliasGraph, Nothing} = nothing;
+                     maxiters = 8000)
     @unpack graph, solvable_graph, var_to_diff, eq_to_diff = state.structure
     neqs = nsrcs(graph)
     nvars = nv(var_to_diff)
     vcolor = falses(nvars)
     ecolor = falses(neqs)
     var_eq_matching = Matching(nvars)
     neqs′ = neqs
+    nnonemptyeqs = count(eq -> !isempty(𝑠neighbors(graph, eq)) && eq_to_diff[eq] === nothing,
+                         1:neqs′)
+
+    # Allow matching for the highest differentiated variable that
+    # currently appears in an equation (or implicit equation in a side ag)
+    varwhitelist = falses(nvars)
+    for var in 1:nvars
+        if var_to_diff[var] === nothing && !varwhitelist[var]
+            while isempty(𝑑neighbors(graph, var)) && (ag === nothing || !haskey(ag, var))
+                var′ = invview(var_to_diff)[var]
+                var′ === nothing && break
+                var = var′
+            end
+            if !isempty(𝑑neighbors(graph, var))
+                if ag !== nothing && haskey(ag, var)
+                    # TODO: remove lower diff vars from whitelist
+                    c, a = ag[var]
+                    iszero(c) || (varwhitelist[a] = true)
+                else
+                    varwhitelist[var] = true
+                end
+            end
+        end
+    end
+
+    if nnonemptyeqs > count(varwhitelist)
+        throw(InvalidSystemException("System is structurally singular"))
+    end
+
     for k in 1:neqs′
         eq′ = k
+        eq_to_diff[eq′] === nothing || continue
         isempty(𝑠neighbors(graph, eq′)) && continue
         pathfound = false
         # In practice, `maxiters=8000` should never be reached, otherwise, the
@@ -93,7 +124,6 @@ function pantelides!(state::TransformationState; maxiters = 8000)
             #
             # the derivatives and algebraic variables are zeros in the variable
             # association list
-            varwhitelist = var_to_diff .== nothing
             resize!(vcolor, nvars)
             fill!(vcolor, false)
             resize!(ecolor, neqs)
@@ -103,11 +133,16 @@ function pantelides!(state::TransformationState; maxiters = 8000)
             pathfound && break # terminating condition
             for var in eachindex(vcolor)
                 vcolor[var] || continue
-                # introduce a new variable
-                nvars += 1
-                var_diff = var_derivative!(state, var)
-                push!(var_eq_matching, unassigned)
-                @assert length(var_eq_matching) == var_diff
+                if var_to_diff[var] === nothing
+                    # introduce a new variable
+                    nvars += 1
+                    var_diff = var_derivative!(state, var)
+                    push!(var_eq_matching, unassigned)
+                    push!(varwhitelist, false)
+                    @assert length(var_eq_matching) == var_diff
+                end
+                varwhitelist[var] = false
+                varwhitelist[var_to_diff[var]] = true
             end
 
             for eq in eachindex(ecolor)

src/structural_transformation/partial_state_selection.jl

@@ -1,6 +1,7 @@
-function partial_state_selection_graph!(state::TransformationState)
+function partial_state_selection_graph!(state::TransformationState;
+                                        ag::Union{AliasGraph, Nothing} = nothing)
     find_solvables!(state; allow_symbolic = true)
-    var_eq_matching = complete(pantelides!(state))
+    var_eq_matching = complete(pantelides!(state, ag))
     complete!(state.structure)
     partial_state_selection_graph!(state.structure, var_eq_matching)
 end
@@ -86,6 +87,9 @@ function pss_graph_modia!(structure::SystemStructure, var_eq_matching, varlevel,
             filter!(eq -> invview(ict.graph.matching)[eq] === unassigned, to_tear_eqs)
             tearEquations!(ict, solvable_graph.fadjlist, to_tear_eqs, BitSet(to_tear_vars),
                            nothing)
+            for var in to_tear_vars
+                var_eq_matching[var] = unassigned
+            end
             for var in to_tear_vars
                 var_eq_matching[var] = ict.graph.matching[var]
             end
@@ -94,7 +98,10 @@ function pss_graph_modia!(structure::SystemStructure, var_eq_matching, varlevel,
         end
     end
     for var in 1:ndsts(graph)
-        if varlevel[var] !== 0 && var_eq_matching[var] === unassigned
+        dv = var_to_diff[var]
+        # If `var` is not algebraic (not differentiated nor a dummy derivative),
+        # then it's a SelectedState
+        if !(dv === nothing || (varlevel[dv] !== 0 && var_eq_matching[dv] === unassigned))
             var_eq_matching[var] = SelectedState()
         end
     end
@@ -116,12 +123,15 @@ function partial_state_selection_graph!(structure::SystemStructure, var_eq_match
     end
 
     varlevel = map(1:ndsts(graph)) do var
-        level = 0
+        graph_level = level = 0
         while var_to_diff[var] !== nothing
             var = var_to_diff[var]
             level += 1
+            if !isempty(𝑑neighbors(graph, var))
+                graph_level = level
+            end
         end
-        level
+        graph_level
     end
 
     inv_varlevel = map(1:ndsts(graph)) do var
@@ -135,7 +145,7 @@ function partial_state_selection_graph!(structure::SystemStructure, var_eq_match
 
     # TODO: Should pantelides just return this?
     for var in 1:ndsts(graph)
-        if var_to_diff[var] !== nothing
+        if varlevel[var] !== 0
             var_eq_matching[var] = unassigned
         end
     end
@@ -147,11 +157,13 @@ function partial_state_selection_graph!(structure::SystemStructure, var_eq_match
     var_eq_matching
 end
 
-function dummy_derivative_graph!(state::TransformationState, jac = nothing, state_priority = nothing; kwargs...)
+function dummy_derivative_graph!(state::TransformationState, jac = nothing,
+                                 (ag, diff_va) = (nothing, nothing);
+                                 state_priority = nothing, kwargs...)
     state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
-    var_eq_matching = complete(pantelides!(state))
+    var_eq_matching = complete(pantelides!(state, ag))
     complete!(state.structure)
-    dummy_derivative_graph!(state.structure, var_eq_matching, jac, state_priority)
+    dummy_derivative_graph!(state.structure, var_eq_matching, jac, (ag, diff_va), state_priority)
 end
 
 function compute_diff_level(diff_to_x)
@@ -171,7 +183,8 @@ function compute_diff_level(diff_to_x)
     return xlevel, maxlevel
 end
 
-function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, jac, state_priority)
+function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, jac,
+                                 (ag, diff_va), state_priority)
     @unpack eq_to_diff, var_to_diff, graph = structure
     diff_to_eq = invview(eq_to_diff)
     diff_to_var = invview(var_to_diff)
@@ -240,6 +253,18 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
             vars = [diff_to_var[var] for var in vars if diff_to_var[var] !== nothing]
         end
     end
+    if diff_va !== nothing
+        n_dummys = length(dummy_derivatives)
+        needed = count(x -> x isa Int, diff_to_eq) - n_dummys
+        n = 0
+        for v in diff_va
+            c, a = ag[v]
+            n += 1
+            push!(dummy_derivatives, iszero(c) ? v : a)
+            needed == n && break
+            continue
+        end
+    end
 
     dummy_derivatives_set = BitSet(dummy_derivatives)
     # We can eliminate variables that are not a selected state (differential
@@ -249,6 +274,9 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
         dummy_derivatives_set = dummy_derivatives_set
 
         v -> begin
+            if ag !== nothing
+                haskey(ag, v) && return false
+            end
             dv = var_to_diff[v]
             dv === nothing || dv in dummy_derivatives_set
         end
@@ -264,7 +292,8 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
                                        Union{Unassigned, SelectedState};
                                        varfilter = can_eliminate)
     for v in eachindex(var_eq_matching)
-        can_eliminate(v) && continue
+        dv = var_to_diff[v]
+        (dv === nothing || dv in dummy_derivatives_set) && continue
         var_eq_matching[v] = SelectedState()
     end
 

src/structural_transformation/symbolics_tearing.jl

@@ -567,7 +567,7 @@ function tearing_reassemble(state::TearingState, var_eq_matching; simplify = fal
     eq_to_diff = new_eq_to_diff
     diff_to_var = invview(var_to_diff)
 
-    @set! state.structure.graph = graph
+    @set! state.structure.graph = complete(graph)
     @set! state.structure.var_to_diff = var_to_diff
     @set! state.structure.eq_to_diff = eq_to_diff
     @set! state.fullvars = fullvars = fullvars[invvarsperm]
@@ -669,6 +669,6 @@ function dummy_derivative(sys, state = TearingState(sys); simplify = false, kwar
             p
         end
     end
-    var_eq_matching = dummy_derivative_graph!(state, jac, state_priority; kwargs...)
+    var_eq_matching = dummy_derivative_graph!(state, jac; state_priority, kwargs...)
     tearing_reassemble(state, var_eq_matching; simplify = simplify)
 end

src/systems/abstractsystem.jl

@@ -405,7 +405,7 @@ function namespace_assignment(eq::Assignment, sys)
     Assignment(_lhs, _rhs)
 end
 
-function namespace_expr(O, sys, n = nameof(sys)) where {T}
+function namespace_expr(O, sys, n = nameof(sys))
     ivs = independent_variables(sys)
     O = unwrap(O)
     if any(isequal(O), ivs)

src/systems/alias_elimination.jl

@@ -6,7 +6,10 @@ const KEEP = typemin(Int)
 
 function alias_eliminate_graph!(state::TransformationState)
     mm = linear_subsys_adjmat(state)
-    size(mm, 1) == 0 && return AliasGraph(ndsts(state.structure.graph)), mm, BitSet() # No linear subsystems
+    if size(mm, 1) == 0
+        ag = AliasGraph(ndsts(state.structure.graph))
+        return ag, mm, ag, mm, BitSet() # No linear subsystems
+    end
 
     @unpack graph, var_to_diff = state.structure
 
@@ -39,7 +42,7 @@ alias_elimination(sys) = alias_elimination!(TearingState(sys; quick_cancel = tru
 function alias_elimination!(state::TearingState)
     sys = state.sys
     complete!(state.structure)
-    ag, mm, updated_diff_vars = alias_eliminate_graph!(state)
+    ag, mm, complete_ag, complete_mm, updated_diff_vars = alias_eliminate_graph!(state)
     isempty(ag) && return sys
 
     fullvars = state.fullvars
@@ -543,6 +546,7 @@ function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
     #
     nvars = ndsts(graph)
     ag = AliasGraph(nvars)
+    complete_ag = AliasGraph(nvars)
     mm, echelon_mm = simple_aliases!(ag, graph, var_to_diff, mm_orig)
 
     # Step 3: Handle differentiated variables
@@ -702,6 +706,7 @@ function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
             while (iv = diff_to_var[v]) in zero_vars_set
                 v = iv
             end
+            complete_ag[v] = 0
             if diff_to_var[v] === nothing # `v` is reducible
                 dag[v] = 0
             end
@@ -729,6 +734,13 @@ function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
     # Step 4: Merge dag and ag
     removed_aliases = BitSet()
     merged_ag = AliasGraph(nvars)
+    for (v, (c, a)) in dag
+        complete_ag[v] = c => a
+    end
+    for (v, (c, a)) in ag
+        (processed[v] || (!iszero(a) && processed[a])) && continue
+        complete_ag[v] = c => a
+    end
     for (v, (c, a)) in dag
         # D(x) ~ D(y) cannot be removed if x and y are not aliases
         if v != a && !iszero(a)
@@ -789,7 +801,8 @@ function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
         update_graph_neighbors!(graph, ag)
     end
 
-    return ag, mm, updated_diff_vars
+    complete_mm = reduce!(copy(echelon_mm), mm_orig, complete_ag, size(echelon_mm, 1))
+    return ag, mm, complete_ag, complete_mm, updated_diff_vars
 end
 
 function update_graph_neighbors!(graph, ag)