Merge branch 'SciML:master' into func-affect

Author: dodoplus

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "8.16.0"
+version = "8.17.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -99,4 +99,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AmplNLWriter", "BenchmarkTools", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "Random", "ReferenceTests", "SafeTestsets", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials"]
+test = ["AmplNLWriter", "BenchmarkTools", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "Random", "ReferenceTests", "SafeTestsets", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials"]

src/bipartite_graph.jl

@@ -6,7 +6,7 @@ export BipartiteEdge, BipartiteGraph, DiCMOBiGraph, Unassigned, unassigned,
 
 export 𝑠vertices, 𝑑vertices, has_𝑠vertex, has_𝑑vertex, 𝑠neighbors, 𝑑neighbors,
        𝑠edges, 𝑑edges, nsrcs, ndsts, SRC, DST, set_neighbors!, invview,
-       complete
+       complete, delete_srcs!, delete_dsts!
 
 using DocStringExtensions
 using UnPack
@@ -49,6 +49,10 @@ end
 function Matching(m::Int)
     Matching{Unassigned}(Union{Int, Unassigned}[unassigned for _ in 1:m], nothing)
 end
+function Matching{U}(m::Int) where {U}
+    Matching{Union{Unassigned, U}}(Union{Int, Unassigned, U}[unassigned for _ in 1:m],
+                                   nothing)
+end
 
 Base.size(m::Matching) = Base.size(m.match)
 Base.getindex(m::Matching, i::Integer) = m.match[i]
@@ -65,9 +69,9 @@ function Base.setindex!(m::Matching{U}, v::Union{Integer, U}, i::Integer) where
     return m.match[i] = v
 end
 
-function Base.push!(m::Matching{U}, v::Union{Integer, U}) where {U}
+function Base.push!(m::Matching, v)
     push!(m.match, v)
-    if v !== unassigned && m.inv_match !== nothing
+    if v isa Integer && m.inv_match !== nothing
         m.inv_match[v] = length(m.match)
     end
 end
@@ -346,8 +350,8 @@ vertices, subject to the constraint that vertices for which `srcfilter` or `dstf
 return `false` may not be matched.
 """
 function maximal_matching(g::BipartiteGraph, srcfilter = vsrc -> true,
-                          dstfilter = vdst -> true)
-    matching = Matching(ndsts(g))
+                          dstfilter = vdst -> true, ::Type{U} = Unassigned) where {U}
+    matching = Matching{U}(ndsts(g))
     foreach(Iterators.filter(srcfilter, 𝑠vertices(g))) do vsrc
         construct_augmenting_path!(matching, g, vsrc, dstfilter)
     end
@@ -420,11 +424,15 @@ function Graphs.add_vertex!(g::BipartiteGraph{T}, type::VertType) where {T}
     return true  # vertex successfully added
 end
 
-function set_neighbors!(g::BipartiteGraph, i::Integer, new_neighbors::AbstractVector)
+function set_neighbors!(g::BipartiteGraph, i::Integer, new_neighbors)
     old_neighbors = g.fadjlist[i]
     old_nneighbors = length(old_neighbors)
     new_nneighbors = length(new_neighbors)
-    g.fadjlist[i] = new_neighbors
+    if iszero(new_nneighbors) # this handles Tuple as well
+        empty!(g.fadjlist[i])
+    else
+        g.fadjlist[i] = new_neighbors
+    end
     g.ne += new_nneighbors - old_nneighbors
     if isa(g.badjlist, AbstractVector)
         for n in old_neighbors
@@ -440,6 +448,14 @@ function set_neighbors!(g::BipartiteGraph, i::Integer, new_neighbors::AbstractVe
     end
 end
 
+function delete_srcs!(g::BipartiteGraph, srcs)
+    for s in srcs
+        set_neighbors!(g, s, ())
+    end
+    g
+end
+delete_dsts!(g::BipartiteGraph, srcs) = delete_srcs!(invview(g), srcs)
+
 ###
 ### Edges iteration
 ###

src/inputoutput.jl

@@ -185,8 +185,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,
-                                                 check_bound = false, kwargs...)
+    sys, diff_idxs, alge_idxs = io_preprocessing(sys, inputs, []; simplify, kwargs...)
 
     dvs = states(sys)
     ps = parameters(sys)

src/structural_transformation/bipartite_tearing/modia_tearing.jl

@@ -9,46 +9,90 @@ function try_assign_eq!(ict::IncrementalCycleTracker, vj::Integer, eq::Integer)
     end
 end
 
-function tearEquations!(ict::IncrementalCycleTracker, Gsolvable, es::Vector{Int},
-                        vs::Vector{Int})
+function try_assign_eq!(ict::IncrementalCycleTracker, vars, v_active, eq::Integer,
+                        condition::F = _ -> true) where {F}
     G = ict.graph
-    vActive = BitSet(vs)
+    for vj in vars
+        (vj in v_active && G.matching[vj] === unassigned && condition(vj)) || continue
+        try_assign_eq!(ict, vj, eq) && return true
+    end
+    return false
+end
 
+function tearEquations!(ict::IncrementalCycleTracker, Gsolvable, es::Vector{Int},
+                        v_active::BitSet, isder′::F) where {F}
+    check_der = isder′ !== nothing
+    if check_der
+        has_der = Ref(false)
+        isder = let has_der = has_der, isder′ = isder′
+            v -> begin
+                r = isder′(v)
+                has_der[] |= r
+                r
+            end
+        end
+    end
     for eq in es  # iterate only over equations that are not in eSolvedFixed
-        for vj in Gsolvable[eq]
-            if G.matching[vj] === unassigned && (vj in vActive)
-                r = try_assign_eq!(ict, vj, eq)
-                r && break
+        vs = Gsolvable[eq]
+        if check_der
+            # if there're differentiated variables, then only consider them
+            try_assign_eq!(ict, vs, v_active, eq, isder)
+            if has_der[]
+                has_der[] = false
+                continue
             end
         end
+        try_assign_eq!(ict, vs, v_active, eq)
     end
 
     return ict
 end
 
-function tear_graph_block_modia!(var_eq_matching, graph, solvable_graph, eqs, vars)
+function tear_graph_block_modia!(var_eq_matching, graph, solvable_graph, eqs, vars,
+                                 isder::F) where {F}
     ict = IncrementalCycleTracker(DiCMOBiGraph{true}(graph); dir = :in)
-    tearEquations!(ict, solvable_graph.fadjlist, eqs, vars)
+    tearEquations!(ict, solvable_graph.fadjlist, eqs, vars, isder)
     for var in vars
         var_eq_matching[var] = ict.graph.matching[var]
     end
     return nothing
 end
 
-function tear_graph_modia(structure::SystemStructure; varfilter = v -> true,
-                          eqfilter = eq -> true)
+function tear_graph_modia(structure::SystemStructure, isder::F = nothing,
+                          ::Type{U} = Unassigned;
+                          varfilter::F2 = v -> true,
+                          eqfilter::F3 = eq -> true) where {F, U, F2, F3}
+    # It would be possible here to simply iterate over all variables and attempt to
+    # use tearEquations! to produce a matching that greedily selects the minimal
+    # number of torn variables. However, we can do this process faster if we first
+    # compute the strongly connected components. In the absence of cycles and
+    # non-solvability, a maximal matching on the original graph will give us an
+    # optimal assignment. However, even with cycles, we can use the maximal matching
+    # to give us a good starting point for a good matching and then proceed to
+    # reverse edges in each scc to improve the solution. Note that it is possible
+    # to have optimal solutions that cannot be found by this process. We will not
+    # find them here [TODO: It would be good to have an explicit example of this.]
+
     @unpack graph, solvable_graph = structure
-    var_eq_matching = complete(maximal_matching(graph, eqfilter, varfilter))
+    var_eq_matching = complete(maximal_matching(graph, eqfilter, varfilter, U))
     var_sccs::Vector{Union{Vector{Int}, Int}} = find_var_sccs(graph, var_eq_matching)
 
+    ieqs = Int[]
+    filtered_vars = BitSet()
     for vars in var_sccs
-        filtered_vars = filter(varfilter, vars)
-        ieqs = Int[var_eq_matching[v]
-                   for v in filtered_vars if var_eq_matching[v] !== unassigned]
         for var in vars
+            if varfilter(var)
+                push!(filtered_vars, var)
+                if var_eq_matching[var] !== unassigned
+                    push!(ieqs, var_eq_matching[var])
+                end
+            end
             var_eq_matching[var] = unassigned
         end
-        tear_graph_block_modia!(var_eq_matching, graph, solvable_graph, ieqs, filtered_vars)
+        tear_graph_block_modia!(var_eq_matching, graph, solvable_graph, ieqs, filtered_vars,
+                                isder)
+        empty!(ieqs)
+        empty!(filtered_vars)
     end
 
     return var_eq_matching

src/structural_transformation/codegen.jl

@@ -6,7 +6,6 @@ const MAX_INLINE_NLSOLVE_SIZE = 8
 
 function torn_system_with_nlsolve_jacobian_sparsity(state, var_eq_matching, var_sccs,
                                                     nlsolve_scc_idxs, eqs_idxs, states_idxs)
-    fullvars = state.fullvars
     graph = state.structure.graph
 
     # The sparsity pattern of `nlsolve(f, u, p)` w.r.t `p` is difficult to
@@ -71,7 +70,6 @@ function torn_system_with_nlsolve_jacobian_sparsity(state, var_eq_matching, var_
 
     var2idx = Dict{Int, Int}(v => i for (i, v) in enumerate(states_idxs))
     eqs2idx = Dict{Int, Int}(v => i for (i, v) in enumerate(eqs_idxs))
-    nlsolve_vars_set = BitSet(nlsolve_vars)
 
     I = Int[]
     J = Int[]

src/structural_transformation/partial_state_selection.jl

@@ -84,7 +84,8 @@ function pss_graph_modia!(structure::SystemStructure, var_eq_matching, varlevel,
             end
             filter!(var -> ict.graph.matching[var] === unassigned, to_tear_vars)
             filter!(eq -> invview(ict.graph.matching)[eq] === unassigned, to_tear_eqs)
-            tearEquations!(ict, solvable_graph.fadjlist, to_tear_eqs, to_tear_vars)
+            tearEquations!(ict, solvable_graph.fadjlist, to_tear_eqs, BitSet(to_tear_vars),
+                           nothing)
             for var in to_tear_vars
                 var_eq_matching[var] = ict.graph.matching[var]
             end
@@ -146,57 +147,51 @@ function partial_state_selection_graph!(structure::SystemStructure, var_eq_match
     var_eq_matching
 end
 
-function dummy_derivative_graph!(state::TransformationState, jac = nothing)
+function dummy_derivative_graph!(state::TransformationState, jac = nothing; kwargs...)
+    state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
     var_eq_matching = complete(pantelides!(state))
     complete!(state.structure)
     dummy_derivative_graph!(state.structure, var_eq_matching, jac)
 end
 
-function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, jac)
-    @unpack eq_to_diff, var_to_diff, graph = structure
-    diff_to_eq = invview(eq_to_diff)
-    diff_to_var = invview(var_to_diff)
-    invgraph = invview(graph)
-
-    neqs = nsrcs(graph)
-    eqlevel = zeros(Int, neqs)
+function compute_diff_level(diff_to_x)
+    nxs = length(diff_to_x)
+    xlevel = zeros(Int, nxs)
     maxlevel = 0
-    for i in 1:neqs
+    for i in 1:nxs
         level = 0
-        eq = i
-        while diff_to_eq[eq] !== nothing
-            eq = diff_to_eq[eq]
+        x = i
+        while diff_to_x[x] !== nothing
+            x = diff_to_x[x]
             level += 1
         end
         maxlevel = max(maxlevel, level)
-        eqlevel[i] = level
+        xlevel[i] = level
     end
+    return xlevel, maxlevel
+end
 
-    nvars = ndsts(graph)
-    varlevel = zeros(Int, nvars)
-    for i in 1:nvars
-        level = 0
-        var = i
-        while diff_to_var[var] !== nothing
-            var = diff_to_var[var]
-            level += 1
-        end
-        maxlevel = max(maxlevel, level)
-        varlevel[i] = level
-    end
+function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, jac)
+    @unpack eq_to_diff, var_to_diff, graph = structure
+    diff_to_eq = invview(eq_to_diff)
+    diff_to_var = invview(var_to_diff)
+    invgraph = invview(graph)
+
+    eqlevel, _ = compute_diff_level(diff_to_eq)
 
     var_sccs = find_var_sccs(graph, var_eq_matching)
-    eqcolor = falses(neqs)
+    eqcolor = falses(nsrcs(graph))
     dummy_derivatives = Int[]
     col_order = Int[]
+    nvars = ndsts(graph)
     for vars in var_sccs
         eqs = [var_eq_matching[var] for var in vars if var_eq_matching[var] !== unassigned]
         isempty(eqs) && continue
-        maxlevel = maximum(map(x -> eqlevel[x], eqs))
+        maxlevel = maximum(Base.Fix1(getindex, eqlevel), eqs)
         iszero(maxlevel) && continue
 
         rank_matching = Matching(nvars)
-        for level in maxlevel:-1:1
+        for _ in maxlevel:-1:1
             eqs = filter(eq -> diff_to_eq[eq] !== nothing, eqs)
             nrows = length(eqs)
             iszero(nrows) && break
@@ -220,8 +215,10 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
             else
                 rank = 0
                 for var in vars
+                    # We need `invgraph` here because we are matching from
+                    # variables to equations.
                     pathfound = construct_augmenting_path!(rank_matching, invgraph, var,
-                                                           eq -> eq in eqs_set, eqcolor)
+                                                           Base.Fix2(in, eqs_set), eqcolor)
                     pathfound || continue
                     push!(dummy_derivatives, var)
                     rank += 1
@@ -239,5 +236,32 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
         end
     end
 
-    dummy_derivatives
+    dummy_derivatives_set = BitSet(dummy_derivatives)
+    # We can eliminate variables that are not a selected state (differential
+    # variables). Selected states are differentiated variables that are not
+    # dummy derivatives.
+    can_eliminate = let var_to_diff = var_to_diff,
+        dummy_derivatives_set = dummy_derivatives_set
+
+        v -> begin
+            dv = var_to_diff[v]
+            dv === nothing || dv in dummy_derivatives_set
+        end
+    end
+
+    # We don't want tearing to give us `y_t ~ D(y)`, so we skip equations with
+    # actually differentiated variables.
+    isdiffed = let diff_to_var = diff_to_var, dummy_derivatives_set = dummy_derivatives_set
+        v -> diff_to_var[v] !== nothing && !(v in dummy_derivatives_set)
+    end
+
+    var_eq_matching = tear_graph_modia(structure, isdiffed,
+                                       Union{Unassigned, SelectedState};
+                                       varfilter = can_eliminate)
+    for v in eachindex(var_eq_matching)
+        can_eliminate(v) && continue
+        var_eq_matching[v] = SelectedState()
+    end
+
+    return var_eq_matching
 end