Merge pull request #1328 from Keno/kf/graphsscc

Replace open coded tarjan algorithm by Graphs.jl version

Author: YingboMa

src/bipartite_graph.jl

@@ -1,6 +1,6 @@
 module BipartiteGraphs
 
-export BipartiteEdge, BipartiteGraph
+export BipartiteEdge, BipartiteGraph, DiCMOBiGraph, Unassigned, unassigned
 
 export 𝑠vertices, 𝑑vertices, has_𝑠vertex, has_𝑑vertex, 𝑠neighbors, 𝑑neighbors,
        𝑠edges, 𝑑edges, nsrcs, ndsts, SRC, DST
@@ -11,6 +11,12 @@ using SparseArrays
 using Graphs
 using Setfield
 
+### Matching
+struct Unassigned
+    global unassigned
+    const unassigned = Unassigned.instance
+end
+
 ###
 ### Edges & Vertex
 ###
@@ -269,4 +275,56 @@ function Graphs.incidence_matrix(g::BipartiteGraph, val=true)
     S = sparse(I, J, val, nsrcs(g), ndsts(g))
 end
 
+
+"""
+    struct DiCMOBiGraph
+
+This data structure implements a "directed, contracted, matching-oriented" view of an
+original (undirected) bipartite graph. In particular, it performs two largely
+orthogonal functions.
+
+1. It pairs an undirected bipartite graph with a matching of destination vertex.
+
+    This matching is used to induce an orientation on the otherwise undirected graph:
+    Matched edges pass from destination to source, all other edges pass in the opposite
+    direction.
+
+2. It exposes the graph view obtained by contracting the destination vertices into
+   the source edges.
+
+The result of this operation is an induced, directed graph on the source vertices.
+The resulting graph has a few desirable properties. In particular, this graph
+is acyclic if and only if the induced directed graph on the original bipartite
+graph is acyclic.
+"""
+struct DiCMOBiGraph{I, G<:BipartiteGraph{I}, M} <: Graphs.AbstractGraph{I}
+    graph::G
+    matching::M
+end
+Graphs.is_directed(::Type{<:DiCMOBiGraph}) = true
+Graphs.nv(g::DiCMOBiGraph) = nsrcs(g.graph)
+Graphs.vertices(g::DiCMOBiGraph) = 1:nsrcs(g.graph)
+
+struct CMOOutNeighbors{V}
+    g::DiCMOBiGraph
+    v::V
+end
+Graphs.outneighbors(g::DiCMOBiGraph, v) = CMOOutNeighbors(g, v)
+Base.iterate(c::CMOOutNeighbors) = iterate(c, (c.g.graph.fadjlist[c.v],))
+function Base.iterate(c::CMOOutNeighbors, (l, state...))
+    while true
+        r = iterate(l, state...)
+        r === nothing && return nothing
+        # If this is a matched edge, skip it, it's reversed in the induced
+        # directed graph. Otherwise, if there is no matching for this destination
+        # edge, also skip it, since it got delted in the contraction.
+        vdst = c.g.matching[r[1]]
+        if vdst === c.v || vdst === unassigned
+            state = (r[2],)
+            continue
+        end
+        return vdst, (l, r[2])
+    end
+end
+
 end # module

src/structural_transformation/StructuralTransformations.jl

@@ -1,8 +1,5 @@
 module StructuralTransformations
 
-const UNVISITED = typemin(Int)
-const UNASSIGNED = typemin(Int)
-
 using Setfield: @set!, @set
 using UnPack: @unpack
 

src/structural_transformation/pantelides.jl

@@ -62,7 +62,7 @@ function pantelides_reassemble(sys::ODESystem, eqassoc, assign)
     end
 
     final_vars = unique(filter(x->!(operation(x) isa Differential), fullvars))
-    final_eqs = map(identity, filter(x->value(x.lhs) !== nothing, out_eqs[sort(filter(x->x != UNASSIGNED, assign))]))
+    final_eqs = map(identity, filter(x->value(x.lhs) !== nothing, out_eqs[sort(filter(x->x !== unassigned, assign))]))
 
     @set! sys.eqs = final_eqs
     @set! sys.states = final_vars
@@ -84,7 +84,7 @@ function pantelides!(sys::ODESystem; maxiters = 8000)
     nvars = length(varassoc)
     vcolor = falses(nvars)
     ecolor = falses(neqs)
-    assign = fill(UNASSIGNED, nvars)
+    assign = Union{Unassigned, Int}[unassigned for _ = 1:nvars]
     eqassoc = fill(0, neqs)
     neqs′ = neqs
     D = Differential(iv)
@@ -112,7 +112,7 @@ function pantelides!(sys::ODESystem; maxiters = 8000)
                 # the new variable is the derivative of `var`
                 varassoc[var] = nvars
                 push!(varassoc, 0)
-                push!(assign, UNASSIGNED)
+                push!(assign, unassigned)
             end
 
             for eq in eachindex(ecolor); ecolor[eq] || continue

src/structural_transformation/utils.jl

@@ -12,7 +12,7 @@ function find_augmenting_path(g, eq, assign, varwhitelist, vcolor=falses(ndsts(g
 
     # if a `var` is unassigned and the edge `eq <=> var` exists
     for var in 𝑠neighbors(g, eq)
-        if (varwhitelist === nothing || varwhitelist[var]) && assign[var] == UNASSIGNED
+        if (varwhitelist === nothing || varwhitelist[var]) && assign[var] === unassigned
             assign[var] = eq
             return true
         end
@@ -37,7 +37,7 @@ Find equation-variable bipartite matching. `s.graph` is a bipartite graph.
 """
 matching(s::SystemStructure, varwhitelist=nothing, eqwhitelist=nothing) = matching(s.graph, varwhitelist, eqwhitelist)
 function matching(g::BipartiteGraph, varwhitelist=nothing, eqwhitelist=nothing)
-    assign = fill(UNASSIGNED, ndsts(g))
+    assign = Union{Unassigned, Int}[unassigned for _ = 1:ndsts(g)]
     for eq in 𝑠vertices(g)
         if eqwhitelist !== nothing
             eqwhitelist[eq] || continue
@@ -98,7 +98,7 @@ function check_consistency(sys::AbstractSystem)
             inv_assign = inverse_mapping(assign) # extra equations
             bad_idxs = findall(iszero, @view inv_assign[1:nsrcs(graph)])
         else
-            bad_idxs = findall(isequal(UNASSIGNED), assign)
+            bad_idxs = findall(isequal(unassigned), assign)
         end
         error_reporting(sys, bad_idxs, n_highest_vars, iseqs)
     end
@@ -110,7 +110,7 @@ function check_consistency(sys::AbstractSystem)
 
     unassigned_var = []
     for (vj, eq) in enumerate(extended_assign)
-        if eq === UNASSIGNED
+        if eq === unassigned
             push!(unassigned_var, fullvars[vj])
         end
     end
@@ -149,72 +149,10 @@ gives the undirected bipartite graph a direction. When `assign === nothing`, we
 assume that the ``i``-th variable is assigned to the ``i``-th equation.
 """
 function find_scc(g::BipartiteGraph, assign=nothing)
-    id = 0
-    stack = Int[]
-    components = Vector{Int}[]
-    n = nsrcs(g)
-    onstack = falses(n)
-    lowlink = zeros(Int, n)
-    ids = fill(UNVISITED, n)
-
-    for eq in 𝑠vertices(g)
-        if ids[eq] == UNVISITED
-            id = strongly_connected!(stack, onstack, components, lowlink, ids, g, assign, eq, id)
-        end
-    end
-    return components
-end
-
-"""
-    strongly_connected!(stack, onstack, components, lowlink, ids, g, assign, eq, id)
-
-Use Tarjan's algorithm to find strongly connected components.
-"""
-function strongly_connected!(stack, onstack, components, lowlink, ids, g, assign, eq, id)
-    id += 1
-    lowlink[eq] = ids[eq] = id
-
-    # add `eq` to the stack
-    push!(stack, eq)
-    onstack[eq] = true
-
-    # for `adjeq` in the adjacency list of `eq`
-    for var in 𝑠neighbors(g, eq)
-        if assign === nothing
-            adjeq = var
-        else
-            # assign[var] => the equation that's assigned to var
-            adjeq = assign[var]
-            # skip equations that are not assigned
-            adjeq == UNASSIGNED && continue
-        end
-
-        # if `adjeq` is not yet idsed
-        if ids[adjeq] == UNVISITED # visit unvisited nodes
-            id = strongly_connected!(stack, onstack, components, lowlink, ids, g, assign, adjeq, id)
-        end
-        # at the callback of the DFS
-        if onstack[adjeq]
-            lowlink[eq] = min(lowlink[eq], lowlink[adjeq])
-        end
-    end
-
-    # if we are at a start of a strongly connected component
-    if lowlink[eq] == ids[eq]
-        component = Int[]
-        repeat = true
-        # pop until we are at the start of the strongly connected component
-        while repeat
-            w = pop!(stack)
-            onstack[w] = false
-            lowlink[w] = ids[eq]
-            # put `w` in current component
-            push!(component, w)
-            repeat = w != eq
-        end
-        push!(components, sort!(component))
-    end
-    return id
+    cmog = DiCMOBiGraph(g, assign === nothing ? Base.OneTo(nsrcs(g)) : assign)
+    sccs = Graphs.strongly_connected_components(cmog)
+    foreach(sort!, sccs)
+    return sccs
 end
 
 function sorted_incidence_matrix(sys, val=true; only_algeqs=false, only_algvars=false)
@@ -295,7 +233,7 @@ end
 function inverse_mapping(assign)
     invassign = zeros(Int, length(assign))
     for (i, eq) in enumerate(assign)
-        eq <= 0 && continue
+        eq === unassigned && continue
         invassign[eq] = i
     end
     return invassign

src/systems/systemstructure.jl

@@ -80,7 +80,7 @@ Base.@kwdef struct SystemStructure
     algeqs::BitVector
     graph::BipartiteGraph{Int,Vector{Vector{Int}},Int,Nothing}
     solvable_graph::BipartiteGraph{Int,Vector{Vector{Int}},Int,Nothing}
-    assign::Vector{Int}
+    assign::Vector{Union{Int, Unassigned}}
     inv_assign::Vector{Int}
     scc::Vector{Vector{Int}}
     partitions::Vector{SystemPartition}

test/structural_transformation/index_reduction.jl

@@ -34,7 +34,7 @@ pendulum = ODESystem(eqs, t, [x, y, w, z, T], [L, g], name=:pendulum)
 pendulum = initialize_system_structure(pendulum)
 sss = structure(pendulum)
 @unpack graph, fullvars, varassoc = sss
-@test StructuralTransformations.matching(sss, varassoc .== 0) == map(x -> x == 0 ? StructuralTransformations.UNASSIGNED : x, [1, 2, 3, 4, 0, 0, 0, 0, 0])
+@test StructuralTransformations.matching(sss, varassoc .== 0) == map(x -> x == 0 ? StructuralTransformations.unassigned : x, [1, 2, 3, 4, 0, 0, 0, 0, 0])
 
 sys, assign, eqassoc = StructuralTransformations.pantelides!(pendulum)
 sss = structure(sys)