Replace open coded tarjan algorithm by Graphs.jl version

Besides avoiding redundancy and making used of the more optimized
version in Graphs.jl, I think the extra abstraction also gives an
insight into what exactly the induced orientation of the graph is
that we're using to find strongly connected components.

While we're at it also replace the hardcoded integer sentintel by
a singleton type to align more with how we're doing this elsewhere
in Julia.

Author: Keno

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)