Merge branch 'master' into myb/complete_ag

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,10 @@ 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)
+        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

@@ -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,14 +74,37 @@ 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)), 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
+            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)) || (ag !== nothing && haskey(ag, var))
+                varwhitelist[var] = true
+            end
+        end
+    end
+
+    if nnonemptyeqs > count(varwhitelist)
+        throw(InvalidSystemException("System is structurally singular"))
+    end
+
     for k in 1:neqs′
         eq′ = k
         isempty(𝑠neighbors(graph, eq′)) && continue
@@ -93,7 +116,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 +125,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
@@ -119,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
@@ -138,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

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]

src/systems/connectors.jl

@@ -198,25 +198,39 @@ function connection2set!(connectionsets, namespace, ss, isouter)
     end
 end
 
-function generate_connection_set(sys::AbstractSystem)
+function generate_connection_set(sys::AbstractSystem, find = nothing, replace = nothing)
     connectionsets = ConnectionSet[]
-    sys = generate_connection_set!(connectionsets, sys)
+    sys = generate_connection_set!(connectionsets, sys, find, replace)
     sys, merge(connectionsets)
 end
 
-function generate_connection_set!(connectionsets, sys::AbstractSystem, namespace = nothing)
+function generate_connection_set!(connectionsets, sys::AbstractSystem, find, replace,
+                                  namespace = nothing)
     subsys = get_systems(sys)
 
     isouter = generate_isouter(sys)
     eqs′ = get_eqs(sys)
     eqs = Equation[]
 
     cts = [] # connections
+    extra_states = []
     for eq in eqs′
-        if eq.lhs isa Connection
-            push!(cts, get_systems(eq.rhs))
+        lhs = eq.lhs
+        rhs = eq.rhs
+        if find !== nothing && find(rhs, namespace)
+            neweq, extra_state = replace(rhs, namespace)
+            if extra_state isa AbstractArray
+                append!(extra_states, unwrap.(extra_state))
+            elseif extra_state !== nothing
+                push!(extra_states, extra_state)
+            end
+            neweq isa AbstractArray ? append!(eqs, neweq) : push!(eqs, neweq)
         else
-            push!(eqs, eq) # split connections and equations
+            if lhs isa Number || lhs isa Symbolic
+                push!(eqs, eq) # split connections and equations
+            else
+                push!(cts, get_systems(rhs))
+            end
         end
     end
 
@@ -235,7 +249,10 @@ function generate_connection_set!(connectionsets, sys::AbstractSystem, namespace
     end
 
     # pre order traversal
-    @set! sys.systems = map(s -> generate_connection_set!(connectionsets, s,
+    if !isempty(extra_states)
+        @set! sys.states = [get_states(sys); extra_states]
+    end
+    @set! sys.systems = map(s -> generate_connection_set!(connectionsets, s, find, replace,
                                                           renamespace(namespace, nameof(s))),
                             subsys)
     @set! sys.eqs = eqs
@@ -305,8 +322,9 @@ function generate_connection_equations_and_stream_connections(csets::AbstractVec
     eqs, stream_connections
 end
 
-function expand_connections(sys::AbstractSystem; debug = false, tol = 1e-10)
-    sys, csets = generate_connection_set(sys)
+function expand_connections(sys::AbstractSystem, find = nothing, replace = nothing;
+                            debug = false, tol = 1e-10)
+    sys, csets = generate_connection_set(sys, find, replace)
     ceqs, instream_csets = generate_connection_equations_and_stream_connections(csets)
     _sys = expand_instream(instream_csets, sys; debug = debug, tol = tol)
     sys = flatten(sys, true)

src/systems/systemstructure.jl

@@ -350,7 +350,8 @@ function TearingState(sys; quick_cancel = false, check = true)
     eq_to_diff = DiffGraph(nsrcs(graph))
 
     return TearingState(sys, fullvars,
-                        SystemStructure(var_to_diff, eq_to_diff, graph, nothing), Any[])
+                        SystemStructure(complete(var_to_diff), complete(eq_to_diff),
+                                        complete(graph), nothing), Any[])
 end
 
 function linear_subsys_adjmat(state::TransformationState)

test/structural_transformation/utils.jl

@@ -21,7 +21,7 @@ StructuralTransformations.find_solvables!(state)
 sss = state.structure
 @unpack graph, solvable_graph, var_to_diff = sss
 @test graph.fadjlist == [[1, 7], [2, 8], [3, 5, 9], [4, 6, 9], [5, 6]]
-@test graph.badjlist == 9
+@test length(graph.badjlist) == 9
 @test ne(graph) == nnz(incidence_matrix(graph)) == 12
 @test nv(solvable_graph) == 9 + 5
 let N = nothing