Finish tearing refactor

Author: Keno

src/bipartite_graph.jl

@@ -1,7 +1,8 @@
 module BipartiteGraphs
 
 export BipartiteEdge, BipartiteGraph, DiCMOBiGraph, Unassigned, unassigned,
-        Matching, ResidualCMOGraph
+        Matching, ResidualCMOGraph, InducedCondensationGraph, maximal_matching,
+        construct_augmenting_path!
 
 export 𝑠vertices, 𝑑vertices, has_𝑠vertex, has_𝑑vertex, 𝑠neighbors, 𝑑neighbors,
        𝑠edges, 𝑑edges, nsrcs, ndsts, SRC, DST, set_neighbors!, invview,
@@ -18,6 +19,14 @@ struct Unassigned
     global unassigned
     const unassigned = Unassigned.instance
 end
+# Behaves as a scalar
+Base.length(u::Unassigned) = 1
+Base.size(u::Unassigned) = ()
+Base.iterate(u::Unassigned) = (unassigned, nothing)
+Base.iterate(u::Unassigned, state) = nothing
+
+Base.show(io::IO, ::Unassigned) =
+    printstyled(io, "u"; color=:light_black)
 
 struct Matching{U #=> :Unassigned =#, V<:AbstractVector} <: AbstractVector{Union{U, Int}}
     match::V
@@ -30,6 +39,7 @@ function Matching{V}(m::Matching) where {V}
     Matching{V}(convert(VUT, m.match),
         m.inv_match === nothing ? nothing : convert(VUT, m.inv_match))
 end
+Matching(m::Matching) = m
 Matching{U}(v::V) where {U, V<:AbstractVector} = Matching{U, V}(v, nothing)
 Matching{U}(v::V, iv::Union{V, Nothing}) where {U, V<:AbstractVector} = Matching{U, V}(v, iv)
 Matching(v::V) where {U, V<:AbstractVector{Union{U, Int}}} =
@@ -135,10 +145,13 @@ mutable struct BipartiteGraph{I<:Integer, M} <: Graphs.AbstractGraph{I}
     metadata::M
 end
 BipartiteGraph(ne::Integer, fadj::AbstractVector, badj::Union{AbstractVector,Integer}=maximum(maximum, fadj); metadata=nothing) = BipartiteGraph(ne, fadj, badj, metadata)
+BipartiteGraph(fadj::AbstractVector, badj::Union{AbstractVector,Integer}=maximum(maximum, fadj); metadata=nothing) =
+    BipartiteGraph(mapreduce(length, +, fadj; init=0), fadj, badj, metadata)
 
 @noinline require_complete(g::BipartiteGraph) = g.badjlist isa AbstractVector || throw(ArgumentError("The graph has no back edges. Use `complete`."))
 
 function invview(g::BipartiteGraph)
+    require_complete(g)
     BipartiteGraph(g.ne, g.badjlist, g.fadjlist)
 end
 
@@ -215,7 +228,53 @@ ndsts(g::BipartiteGraph) = length(𝑑vertices(g))
 function Graphs.has_edge(g::BipartiteGraph, edge::BipartiteEdge)
     @unpack src, dst = edge
     (src in 𝑠vertices(g) && dst in 𝑑vertices(g)) || return false  # edge out of bounds
-    insorted(𝑠neighbors(src), dst)
+    insorted(dst, 𝑠neighbors(g, src))
+end
+Base.in(edge::BipartiteEdge, g::BipartiteGraph) = Graphs.has_edge(g, edge)
+
+### Maximal matching
+"""
+    construct_augmenting_path!(m::Matching, g::BipartiteGraph, vsrc, dstfilter, vcolor=falses(ndsts(g)), ecolor=falses(nsrcs(g))) -> path_found::Bool
+
+Try to construct an augmenting path in matching and if such a path is found,
+update the matching accordingly.
+"""
+function construct_augmenting_path!(matching::Matching, g::BipartiteGraph, vsrc, dstfilter, dcolor=falses(ndsts(g)), scolor=falses(nsrcs(g)))
+    scolor[vsrc] = true
+
+    # if a `vdst` is unassigned and the edge `vsrc <=> vdst` exists
+    for vdst in 𝑠neighbors(g, vsrc)
+        if dstfilter(vdst) && matching[vdst] === unassigned
+            matching[vdst] = vsrc
+            return true
+        end
+    end
+
+    # for every `vsrc` such that edge `vsrc <=> vdst` exists and `vdst` is uncolored
+    for vdst in 𝑠neighbors(g, vsrc)
+        (dstfilter(vdst) && !dcolor[vdst]) || continue
+        dcolor[vdst] = true
+        if construct_augmenting_path!(matching, g, matching[vdst], dstfilter, dcolor, scolor)
+            matching[vdst] = vsrc
+            return true
+        end
+    end
+    return false
+end
+
+"""
+    maximal_matching(g::BipartiteGraph, [srcfilter], [dstfilter])
+
+For a bipartite graph `g`, construct a maximal matching of destination to source
+vertices, subject to the constraint that vertices for which `srcfilter` or `dstfilter`,
+return `false` may not be matched.
+"""
+function maximal_matching(g::BipartiteGraph, srcfilter=vsrc->true, dstfilter=vdst->true)
+    matching = Matching(ndsts(g))
+    foreach(Iterators.filter(srcfilter, 𝑠vertices(g))) do vsrc
+        construct_augmenting_path!(matching, g, vsrc, dstfilter)
+    end
+    return matching
 end
 
 ###
@@ -508,4 +567,66 @@ function Graphs.neighbors(rcg::ResidualCMOGraph, v::Integer)
             invview(rcg.matching)[vsrc] === unassigned))
 end
 
+# TODO: Fix the function in Graphs to do this instead
+function Graphs.neighborhood(rcg::ResidualCMOGraph, v::Integer)
+    worklist = Int[v]
+    ns = BitSet()
+    while !isempty(worklist)
+        v′ = popfirst!(worklist)
+        for n in neighbors(rcg, v′)
+            if !(n in ns)
+                push!(ns, n)
+                push!(worklist, n)
+            end
+        end
+    end
+    return ns
+end
+
+"""
+    struct InducedCondensationGraph
+
+For some bipartite-graph and an orientation induced on its destination contraction,
+records the condensation DAG of the digraph formed by the orientation. I.e. this
+is a DAG of connected components formed by the destination vertices of some
+underlying bipartite graph.
+
+N.B.: This graph does not store explicit neighbor relations of the sccs.
+Therefor, the edge multiplicity is derived from the underlying bipartite graph,
+i.e. this graph is not strict.
+"""
+struct InducedCondensationGraph{G <: BipartiteGraph} <: AbstractGraph{Vector{Union{Int, Vector{Int}}}}
+    graph::G
+    # Records the members of a strongly connected component. For efficiency,
+    # trivial sccs (with one vertex member) are stored inline. Note: the sccs
+    # here are stored in topological order.
+    sccs::Vector{Union{Int, Vector{Int}}}
+    # Maps the vertices back to the scc of which they are a part
+    scc_assignment::Vector{Int}
+end
+
+function InducedCondensationGraph(g::BipartiteGraph, sccs::Vector{Union{Int, Vector{Int}}})
+    scc_assignment = Vector{Int}(undef, ndsts(g))
+    for (i, c) in enumerate(sccs)
+        for v in c
+            scc_assignment[v] = i
+        end
+    end
+    InducedCondensationGraph(g, sccs, scc_assignment)
+end
+
+Graphs.is_directed(::Type{<:InducedCondensationGraph}) = true
+Graphs.nv(icg::InducedCondensationGraph) = length(icg.sccs)
+Graphs.vertices(icg::InducedCondensationGraph) = icg.sccs
+
+_neighbors(icg::InducedCondensationGraph, cc::Integer) =
+    Iterators.flatten(Iterators.flatten(rcg.graph.fadjlist[vsrc] for vsrc in rcg.graph.badjlist[v]) for v in icg.sccs[cc])
+
+Graphs.outneighbors(rcg::InducedCondensationGraph, v::Integer) =
+    (scc_assignment[n] for n in _neighbors(rcg, v) if scc_assignment[n] > v)
+
+Graphs.inneighbors(rcg::InducedCondensationGraph, v::Integer) =
+    (scc_assignment[n] for n in _neighbors(rcg, v) if scc_assignment[n] < v)
+
+
 end # module

src/structural_transformation/bipartite_tearing/modia_tearing.jl

@@ -42,14 +42,30 @@ function tearEquations!(ict::IncrementalCycleTracker, Gsolvable, es::Vector{Int}
         end
     end
 
-    vSolved = filter(v->G.matching[v] !== unassigned, topological_sort(ict))
-    inv_matching = Union{Missing, Int}[missing for _ = 1:nv(G)]
-    for (v, eq) in pairs(G.matching)
-        eq === unassigned && continue
-        inv_matching[v] = eq
+    return ict
+end
+
+"""
+    tear_graph_modia(sys) -> sys
+
+Tear the bipartite graph in a system. End users are encouraged to call [`structural_simplify`](@ref)
+instead, which calls this function internally.
+"""
+function tear_graph_modia(graph::BipartiteGraph, solvable_graph::BipartiteGraph; varfilter=v->true, eqfilter=eq->true)
+    var_eq_matching = complete(maximal_matching(graph, eqfilter, varfilter))
+    var_sccs::Vector{Union{Vector{Int}, Int}} = find_var_sccs(graph, var_eq_matching)
+
+    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]
+
+        ict = IncrementalCycleTracker(DiCMOBiGraph{true}(graph); dir=:in)
+        tearEquations!(ict, solvable_graph.fadjlist, ieqs, filtered_vars)
+
+        for var in vars
+            var_eq_matching[var] = ict.graph.matching[var]
+        end
     end
-    eSolved = getindex.(Ref(inv_matching), vSolved)
-    vTear = setdiff(vs, vSolved)
-    eResidue = setdiff(es, eSolved)
-    return (eSolved, vSolved, eResidue, vTear)
+
+    return var_eq_matching
 end

src/structural_transformation/codegen.jl

@@ -2,9 +2,9 @@ using LinearAlgebra
 
 const MAX_INLINE_NLSOLVE_SIZE = 8
 
-function torn_system_jacobian_sparsity(sys)
+function torn_system_jacobian_sparsity(sys, var_eq_matching, var_sccs)
     s = structure(sys)
-    @unpack fullvars, graph, partitions = s
+    @unpack fullvars, graph = s
 
     # The sparsity pattern of `nlsolve(f, u, p)` w.r.t `p` is difficult to
     # determine in general. Consider the "simplest" case, a linear system. We
@@ -44,8 +44,9 @@ function torn_system_jacobian_sparsity(sys)
     # dependencies.
     avars2dvars = Dict{Int,Set{Int}}()
     c = 0
-    for partition in partitions
-        @unpack e_residual, v_residual = partition
+    for scc in var_sccs
+        v_residual = scc
+        e_residual = [var_eq_matching[c] for c in v_residual if var_eq_matching[c] !== unassigned]
         # initialization
         for tvar in v_residual
             avars2dvars[tvar] = Set{Int}()
@@ -94,41 +95,6 @@ function torn_system_jacobian_sparsity(sys)
     sparse(I, J, true)
 end
 
-"""
-    partitions_dag(s::SystemStructure)
-
-Return a DAG (sparse matrix) of partitions to use for parallelism.
-"""
-function partitions_dag(s::SystemStructure)
-    @unpack partitions, graph = s
-
-    # `partvars[i]` contains all the states that appear in `partitions[i]`
-    partvars = map(partitions) do partition
-        ipartvars = Set{Int}()
-        for req in partition.e_residual
-            union!(ipartvars, 𝑠neighbors(graph, req))
-        end
-        ipartvars
-    end
-
-    I, J = Int[], Int[]
-    n = length(partitions)
-    for (i, partition) in enumerate(partitions)
-        for j in i+1:n
-            # The only way for a later partition `j` to depend on an earlier
-            # partition `i` is when `partvars[j]` contains one of tearing
-            # variables of partition `i`.
-            if !isdisjoint(partvars[j], partition.v_residual)
-                # j depends on i
-                push!(I, i)
-                push!(J, j)
-            end
-        end
-    end
-
-    sparse(I, J, true, n, n)
-end
-
 """
     exprs = gen_nlsolve(eqs::Vector{Equation}, vars::Vector, u0map::Dict; checkbounds = true)
 
@@ -205,19 +171,20 @@ function build_torn_function(
     end
 
     s = structure(sys)
-    @unpack fullvars, partitions = s
+    @unpack fullvars = s
+    var_eq_matching, var_sccs = algebraic_variables_scc(sys)
 
     states = map(i->s.fullvars[i], diffvars_range(s))
     mass_matrix_diag = ones(length(states))
     torn_expr = []
     defs = defaults(sys)
 
     needs_extending = false
-    for p in partitions
-        @unpack e_residual, v_residual = p
-        torn_eqs = eqs[e_residual]
-        torn_vars = fullvars[v_residual]
-        if length(e_residual) <= max_inlining_size
+    for scc in var_sccs
+        torn_vars = [s.fullvars[var] for var in scc if var_eq_matching[var] !== unassigned]
+        torn_eqs = [eqs[var_eq_matching[var]] for var in scc if var_eq_matching[var] !== unassigned]
+        isempty(torn_eqs) && continue
+        if length(torn_eqs) <= max_inlining_size
             append!(torn_expr, gen_nlsolve(torn_eqs, torn_vars, defs, checkbounds=checkbounds))
         else
             needs_extending = true
@@ -260,15 +227,15 @@ function build_torn_function(
         observedfun = let sys = sys, dict = Dict()
             function generated_observed(obsvar, u, p, t)
                 obs = get!(dict, value(obsvar)) do
-                    build_observed_function(sys, obsvar, checkbounds=checkbounds)
+                    build_observed_function(sys, obsvar, var_eq_matching, var_sccs, checkbounds=checkbounds)
                 end
                 obs(u, p, t)
             end
         end
 
         ODEFunction{true}(
                           @RuntimeGeneratedFunction(expr),
-                          sparsity = torn_system_jacobian_sparsity(sys),
+                          sparsity = torn_system_jacobian_sparsity(sys, var_eq_matching, var_sccs),
                           syms = syms,
                           observed = observedfun,
                           mass_matrix = mass_matrix,
@@ -277,24 +244,24 @@ function build_torn_function(
 end
 
 """
-    find_solve_sequence(partitions, vars)
+    find_solve_sequence(sccs, vars)
 
 given a set of `vars`, find the groups of equations we need to solve for
 to obtain the solution to `vars`
 """
-function find_solve_sequence(partitions, vars)
-    subset = filter(x -> !isdisjoint(x.v_residual, vars), partitions)
+function find_solve_sequence(sccs, vars)
+    subset = filter(i -> !isdisjoint(sccs[i], vars), 1:length(sccs))
     isempty(subset) && return []
-    vars′ = mapreduce(x->x.v_residual, union, subset)
+    vars′ = mapreduce(i->sccs[i], union, subset)
     if vars′ == vars
         return subset
     else
-        return find_solve_sequence(partitions, vars′)
+        return find_solve_sequence(sccs, vars′)
     end
 end
 
 function build_observed_function(
-        sys, ts;
+        sys, ts, var_eq_matching, var_sccs;
         expression=false,
         output_type=Array,
         checkbounds=true
@@ -311,7 +278,7 @@ function build_observed_function(
     dep_vars = collect(setdiff(vars, ivs))
 
     s = structure(sys)
-    @unpack partitions, fullvars, graph = s
+    @unpack fullvars, graph = s
     diffvars = map(i->fullvars[i], diffvars_range(s))
     algvars = map(i->fullvars[i], algvars_range(s))
 
@@ -339,12 +306,12 @@ function build_observed_function(
     end
 
     varidxs = findall(x->x in required_algvars, fullvars)
-    subset = find_solve_sequence(partitions, varidxs)
+    subset = find_solve_sequence(var_sccs, varidxs)
     if !isempty(subset)
         eqs = equations(sys)
 
-        torn_eqs  = map(idxs-> eqs[idxs.e_residual], subset)
-        torn_vars = map(idxs->fullvars[idxs.v_residual], subset)
+        torn_eqs  = map(i->map(v->eqs[var_eq_matching[v]], var_sccs[i]), subset)
+        torn_vars = map(i->map(v->fullvars[v], var_sccs[i]), subset)
         u0map = defaults(sys)
         solves = gen_nlsolve.(torn_eqs, torn_vars, (u0map,); checkbounds=checkbounds)
     else

src/structural_transformation/pantelides.jl

@@ -101,7 +101,7 @@ function pantelides!(graph, var_to_diff; maxiters = 8000)
             fill!(vcolor, false)
             resize!(ecolor, neqs)
             fill!(ecolor, false)
-            pathfound = find_augmenting_path(graph, eq′, var_eq_matching, varwhitelist, vcolor, ecolor)
+            pathfound = construct_augmenting_path!(var_eq_matching, graph, eq′, v->varwhitelist[v], vcolor, ecolor)
             pathfound && break # terminating condition
             for var in eachindex(vcolor); vcolor[var] || continue
                 # introduce a new variable