Skip to content

Commit af69c71

Browse files
authored
Merge pull request #1338 from Keno/kf/retear
Refactor tearing code
2 parents 52bdc81 + e08e269 commit af69c71

File tree

9 files changed

+304
-370
lines changed

9 files changed

+304
-370
lines changed

src/ModelingToolkit.jl

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@ using SpecialFunctions, NaNMath
1818
using RuntimeGeneratedFunctions
1919
using Base.Threads
2020
using DiffEqCallbacks
21+
using Graphs
2122
import MacroTools: splitdef, combinedef, postwalk, striplines
2223
import Libdl
2324
using DocStringExtensions
@@ -115,6 +116,12 @@ include("parameters.jl")
115116
include("utils.jl")
116117
include("domains.jl")
117118

119+
# Code that should eventually go elsewhere, but is here for fow
120+
include("compat/bareiss.jl")
121+
if !isdefined(Graphs, :IncrementalCycleTracker)
122+
include("compat/incremental_cycles.jl")
123+
end
124+
118125
include("systems/abstractsystem.jl")
119126

120127
include("systems/diffeqs/odesystem.jl")

src/bipartite_graph.jl

Lines changed: 46 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -280,50 +280,77 @@ end
280280
struct DiCMOBiGraph
281281
282282
This data structure implements a "directed, contracted, matching-oriented" view of an
283-
original (undirected) bipartite graph. In particular, it performs two largely
284-
orthogonal functions.
283+
original (undirected) bipartite graph. It has two modes, depending on the `Transposed`
284+
flag, which switches the direction of the induced matching.
285285
286-
1. It pairs an undirected bipartite graph with a matching of destination vertex.
286+
Essentially the graph adapter performs two largely orthogonal functions
287+
[`Transposed == true` differences are indicated in square brackets]:
288+
289+
1. It pairs an undirected bipartite graph with a matching of the destination vertex.
287290
288291
This matching is used to induce an orientation on the otherwise undirected graph:
289-
Matched edges pass from destination to source, all other edges pass in the opposite
290-
direction.
292+
Matched edges pass from destination to source [source to desination], all other edges
293+
pass in the opposite direction.
291294
292-
2. It exposes the graph view obtained by contracting the destination vertices into
293-
the source edges.
295+
2. It exposes the graph view obtained by contracting the destination [source] vertices
296+
along the matched edges.
294297
295-
The result of this operation is an induced, directed graph on the source vertices.
298+
The result of this operation is an induced, directed graph on the source [destination] vertices.
296299
The resulting graph has a few desirable properties. In particular, this graph
297300
is acyclic if and only if the induced directed graph on the original bipartite
298301
graph is acyclic.
302+
299303
"""
300-
struct DiCMOBiGraph{I, G<:BipartiteGraph{I}, M} <: Graphs.AbstractGraph{I}
304+
struct DiCMOBiGraph{Transposed, I, G<:BipartiteGraph{I}, M} <: Graphs.AbstractGraph{I}
301305
graph::G
302306
matching::M
307+
DiCMOBiGraph{Transposed}(g::G, m::M) where {Transposed, I, G<:BipartiteGraph{I}, M} =
308+
new{Transposed, I, G, M}(g, m)
303309
end
310+
DiCMOBiGraph{Transposed}(g::BipartiteGraph) where {Transposed} = DiCMOBiGraph{Transposed}(g, Union{Unassigned, Int}[unassigned for i = 1:ndsts(g)])
304311
Graphs.is_directed(::Type{<:DiCMOBiGraph}) = true
305-
Graphs.nv(g::DiCMOBiGraph) = nsrcs(g.graph)
306-
Graphs.vertices(g::DiCMOBiGraph) = 1:nsrcs(g.graph)
312+
Graphs.nv(g::DiCMOBiGraph{Transposed}) where {Transposed} = Transposed ? ndsts(g.graph) : nsrcs(g.graph)
313+
Graphs.vertices(g::DiCMOBiGraph{Transposed}) where {Transposed} = Transposed ? 𝑑vertices(g.graph) : 𝑠vertices(g.graph)
307314

308-
struct CMOOutNeighbors{V}
309-
g::DiCMOBiGraph
315+
struct CMONeighbors{Transposed, V}
316+
g::DiCMOBiGraph{Transposed}
310317
v::V
318+
CMONeighbors{Transposed}(g::DiCMOBiGraph{Transposed}, v::V) where {Transposed, V} =
319+
new{Transposed, V}(g, v)
311320
end
312-
Graphs.outneighbors(g::DiCMOBiGraph, v) = CMOOutNeighbors(g, v)
313-
Base.iterate(c::CMOOutNeighbors) = iterate(c, (c.g.graph.fadjlist[c.v],))
314-
function Base.iterate(c::CMOOutNeighbors, (l, state...))
321+
Graphs.outneighbors(g::DiCMOBiGraph{false}, v) = CMONeighbors{false}(g, v)
322+
Base.iterate(c::CMONeighbors{false}) = iterate(c, (c.g.graph.fadjlist[c.v],))
323+
function Base.iterate(c::CMONeighbors{false}, (l, state...))
315324
while true
316325
r = iterate(l, state...)
317326
r === nothing && return nothing
318327
# If this is a matched edge, skip it, it's reversed in the induced
319328
# directed graph. Otherwise, if there is no matching for this destination
320329
# edge, also skip it, since it got delted in the contraction.
321-
vdst = c.g.matching[r[1]]
322-
if vdst === c.v || vdst === unassigned
330+
vsrc = c.g.matching[r[1]]
331+
if vsrc === c.v || vsrc === unassigned
332+
state = (r[2],)
333+
continue
334+
end
335+
return vsrc, (l, r[2])
336+
end
337+
end
338+
339+
Graphs.inneighbors(g::DiCMOBiGraph{true}, v) = CMONeighbors{true}(g, v)
340+
function Base.iterate(c::CMONeighbors{true})
341+
vsrc = c.g.matching[c.v]
342+
vsrc === unassigned && return nothing
343+
iterate(c, (c.g.graph.fadjlist[vsrc],))
344+
end
345+
function Base.iterate(c::CMONeighbors{true}, (l, state...))
346+
while true
347+
r = iterate(l, state...)
348+
r === nothing && return nothing
349+
if r[1] === c.v
323350
state = (r[2],)
324351
continue
325352
end
326-
return vdst, (l, r[2])
353+
return r[1], (l, r[2])
327354
end
328355
end
329356

File renamed without changes.

src/compat/incremental_cycles.jl

Lines changed: 227 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,227 @@
1+
using Base.Iterators: repeated
2+
3+
# Abstract Interface
4+
5+
"""
6+
abstract type IncrementalCycleTracker
7+
8+
The supertype for incremental cycle detection problems. The abstract type
9+
constructor IncrementalCycleTracker(G) may be used to automatically select
10+
a specific incremental cycle detection algorithm. See [`add_edge_checked!`](@ref)
11+
for a usage example.
12+
"""
13+
abstract type IncrementalCycleTracker{I} <: AbstractGraph{I} end
14+
15+
function (::Type{IncrementalCycleTracker})(s::AbstractGraph{I}; in_out_reverse=nothing) where {I}
16+
# TODO: Once we have more algorithms, the poly-algorithm decision goes here.
17+
# For now, we only have Algorithm N.
18+
return DenseGraphICT_BFGT_N{something(in_out_reverse, false)}(s)
19+
end
20+
21+
# Cycle Detection Interface
22+
"""
23+
add_edge_checked!([f!,], ict::IncrementalCycleTracker, v, w)
24+
25+
Using the incremental cycle tracker, ict, check whether adding the edge `v=>w`.
26+
Would introduce a cycle in the underlying graph. If so, return false and leave
27+
the ict intact. If not, update the underlying graph and return true.
28+
29+
# Optional `f!` Argument
30+
31+
By default the `add_edge!` function is used to update the underlying graph.
32+
However, for more complicated graphs, users may wish to manually specify the
33+
graph update operation. This may be accomplished by passing the optional `f!`
34+
callback arhgument. This callback is called on the underlying graph when no
35+
cycle is detected and is required to modify the underlying graph in order to
36+
effectuate the proposed edge addition.
37+
38+
# Batched edge additions
39+
40+
Optionally, either `v` or `w` (depending on the `in_out_reverse` flag) may be a
41+
collection of vertices representing a batched addition of vertices sharing a
42+
common source or target more efficiently than individual updates.
43+
44+
## Example
45+
46+
```jldoctest
47+
julia> G = SimpleDiGraph(3)
48+
49+
julia> ict = IncrementalCycleTracker(G)
50+
BFGT_N cycle tracker on {3, 0} directed simple Int64 graph
51+
52+
julia> add_edge_checked!(ict, 1, 2)
53+
true
54+
55+
julia> collect(edges(G))
56+
1-element Vector{Graphs.SimpleGraphs.SimpleEdge{Int64}}:
57+
Edge 1 => 2
58+
59+
julia> add_edge_checked!(ict, 2, 3)
60+
true
61+
62+
julia> collect(edges(G))
63+
2-element Vector{Graphs.SimpleGraphs.SimpleEdge{Int64}}:
64+
Edge 1 => 2
65+
Edge 2 => 3
66+
67+
julia> add_edge_checked!(ict, 3, 1) # Would add a cycle
68+
false
69+
70+
julia> collect(edges(G))
71+
2-element Vector{Graphs.SimpleGraphs.SimpleEdge{Int64}}:
72+
Edge 1 => 2
73+
Edge 2 => 3
74+
```
75+
"""
76+
function add_edge_checked! end
77+
78+
to_edges(v::Integer, w::Integer) = (v=>w,)
79+
to_edges(v::Integer, ws) = zip(repeated(v), ws)
80+
to_edges(vs, w::Integer) = zip(vs, repeated(w))
81+
82+
add_edge_checked!(ict::IncrementalCycleTracker, vs, ws) = add_edge_checked!(ict, vs, ws) do g
83+
foreach(((v, w),)->add_edge!(g, v, w), to_edges(vs, ws))
84+
end
85+
86+
# Utilities
87+
"""
88+
struct TransactionalVector
89+
90+
A vector with one checkpoint that may be reverted to by calling `revert!`. The setpoint itself
91+
is set by calling `commit!`.
92+
"""
93+
struct TransactionalVector{T} <: AbstractVector{T}
94+
v::Vector{T}
95+
log::Vector{Pair{Int, T}}
96+
TransactionalVector(v::Vector{T}) where {T} =
97+
new{T}(v, Vector{Pair{Int, T}}())
98+
end
99+
100+
function commit!(v::TransactionalVector)
101+
empty!(v.log)
102+
return nothing
103+
end
104+
105+
function revert!(vec::TransactionalVector)
106+
for (idx, val) in reverse(vec.log)
107+
vec.v[idx] = val
108+
end
109+
return nothing
110+
end
111+
112+
function Base.setindex!(vec::TransactionalVector, val, idx)
113+
oldval = vec.v[idx]
114+
vec.v[idx] = val
115+
push!(vec.log, idx=>oldval)
116+
return nothing
117+
end
118+
Base.getindex(vec::TransactionalVector, idx) = vec.v[idx]
119+
Base.size(vec) = size(vec.v)
120+
121+
# Specific Algorithms
122+
123+
const bibliography = """
124+
## References
125+
126+
[BFGT15] Michael A. Bender, Jeremy T. Fineman, Seth Gilbert, and Robert E. Tarjan. 2015
127+
A New Approach to Incremental Cycle Detection and Related Problems.
128+
ACM Trans. Algorithms 12, 2, Article 14 (December 2015), 22 pages.
129+
DOI: http://dx.doi.org/10.1145/2756553
130+
"""
131+
132+
## Bender, Algorithm N
133+
134+
"""
135+
struct DenseGraphICT_BFGT_N
136+
137+
Implements the "Naive" (Algorithm N) Bender-Fineman-Gilbert-Tarjan one-way line search incremental cycle detector
138+
for dense graphs from [BFGT15] (Section 3).
139+
140+
$bibliography
141+
"""
142+
struct DenseGraphICT_BFGT_N{InOutReverse, I, G<:AbstractGraph{I}} <: IncrementalCycleTracker{I}
143+
graph::G
144+
levels::TransactionalVector{Int}
145+
DenseGraphICT_BFGT_N{InOutReverse}(g::G) where {InOutReverse, I, G<:AbstractGraph{I}} =
146+
new{InOutReverse, I, G}(g, TransactionalVector(fill(0, nv(g))))
147+
end
148+
function Base.show(io::IO, ict::DenseGraphICT_BFGT_N)
149+
print(io, "BFGT_N cycle tracker on ")
150+
show(io, ict.graph)
151+
end
152+
153+
function topological_sort(ict::DenseGraphICT_BFGT_N{InOutReverse}) where {InOutReverse}
154+
# The ICT levels are a weak topological ordering, so a sort of the levels
155+
# will give a topological sort of the vertices.
156+
perm = sortperm(ict.levels)
157+
InOutReverse && (perm = reverse(perm))
158+
return perm
159+
end
160+
161+
# Even when both `v` and `w` are integer, we know that `v` would come first, so
162+
# we prefer to check for `v` as the cycle vertex in this case.
163+
add_edge_checked!(f!, ict::DenseGraphICT_BFGT_N{false}, v::Integer, ws) =
164+
_check_cycle_add!(f!, ict, to_edges(v, ws), v)
165+
add_edge_checked!(f!, ict::DenseGraphICT_BFGT_N{true}, vs, w::Integer) =
166+
_check_cycle_add!(f!, ict, to_edges(vs, w), w)
167+
168+
### [BFGT15] Algorithm N
169+
#
170+
# Implementation Notes
171+
#
172+
# This is Algorithm N from [BFGT15] (Section 3), plus limited patching support and
173+
# a number of standard tricks. Namely:
174+
#
175+
# 1. Batching is supported as long as there is only a single source or destination
176+
# vertex. General batching is left as an open problem. The reason that the
177+
# single source/dest batching is easy to add is that we know that either the
178+
# source or the destination vertex is guaranteed to be a part of any cycle
179+
# that we may have added. Thus we're guaranteed to encounter one of the two
180+
# verticies in our cycle validation and the rest of the algorithm goes through
181+
# as usual.
182+
# 2. We opportunistically traverse each edge when we see it and only add it
183+
# to the worklist if we know that traversal will recurse further.
184+
# 3. We add some early out checks to detect we're about to do redundant work.
185+
function _check_cycle_add!(f!, ict::DenseGraphICT_BFGT_N{InOutReverse}, edges, v) where {InOutReverse}
186+
g = ict.graph
187+
worklist = Pair{Int, Int}[]
188+
# TODO: In the case where there's a single target vertex, we could saturate
189+
# the level first before we assign it to the tracked vector to save some
190+
# log space.
191+
for (v, w) in edges
192+
InOutReverse && ((v, w) = (w, v))
193+
if ict.levels[v] < ict.levels[w]
194+
continue
195+
end
196+
v == w && return false
197+
ict.levels[w] = ict.levels[v] + 1
198+
push!(worklist, v=>w)
199+
end
200+
while !isempty(worklist)
201+
(x, y) = popfirst!(worklist)
202+
xlevel = ict.levels[x]
203+
ylevel = ict.levels[y]
204+
if xlevel >= ylevel
205+
# The xlevel may have been incremented further since we added this
206+
# edge to the worklist.
207+
ict.levels[y] = ylevel = xlevel + 1
208+
elseif ylevel > xlevel + 1
209+
# Some edge traversal scheduled for later already incremented this
210+
# level past where we would have been. Delay processing until then.
211+
continue
212+
end
213+
for z in (InOutReverse ? inneighbors(g, y) : outneighbors(g, y))
214+
if z == v
215+
revert!(ict.levels)
216+
return false
217+
end
218+
if ylevel >= ict.levels[z]
219+
ict.levels[z] = ylevel + 1
220+
push!(worklist, y=>z)
221+
end
222+
end
223+
end
224+
commit!(ict.levels)
225+
f!(g)
226+
return true
227+
end

src/structural_transformation/StructuralTransformations.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,8 @@ using ModelingToolkit: ODESystem, AbstractSystem,var_from_nested_derivative, Dif
1818
get_structure, defaults, InvalidSystemException,
1919
ExtraEquationsSystemException,
2020
ExtraVariablesSystemException,
21-
get_postprocess_fbody, vars!
21+
get_postprocess_fbody, vars!,
22+
IncrementalCycleTracker, add_edge_checked!, topological_sort
2223

2324
using ModelingToolkit.BipartiteGraphs
2425
using Graphs

0 commit comments

Comments
 (0)