Skip to content

Commit e08e269

Browse files
committed
Address review, vendor Graphs change
1 parent d3b5209 commit e08e269

File tree

6 files changed

+242
-17
lines changed

6 files changed

+242
-17
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")
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

src/structural_transformation/bipartite_tearing/modia_tearing.jl

Lines changed: 6 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,7 @@
1414
################################################
1515

1616
"""
17-
(eSolved, vSolved, eResidue, vTear) = tearEquations!(td, Gsolvable, es, vs; eSolvedFixed=Int[], vSolvedFixed=Int[], vTearFixed=Int[])
17+
(eSolved, vSolved, eResidue, vTear) = tearEquations!(td, Gsolvable, es, vs)
1818
1919
Equations es shall be solved with respect to variables vs. The function returns
2020
the teared equation so that if vTear is given, vSolved can be computed from eSolved
@@ -25,31 +25,23 @@ Variables vs are the union of vSolved and vTear.
2525
2626
Gsolvable defines the variables that can be explicitly solved in every equation without influencing the solution space
2727
(= rank preserving operation).
28-
29-
eSolvedFixed/vSolvedFixed must be a DAG starting at eSolvedFixed/SolvedFixed[1]
3028
"""
31-
function tearEquations!(ict::IncrementalCycleTracker, Gsolvable, es::Vector{Int}, vs::Vector{Int};
32-
eSolvedFixed::Vector{Int}=Int[], vSolvedFixed::Vector{Int}=Int[], vTearFixed::Vector{Int}=Int[])
29+
function tearEquations!(ict::IncrementalCycleTracker, Gsolvable, es::Vector{Int}, vs::Vector{Int})
3330
G = ict.graph
3431
vActive = BitSet(vs)
35-
vMatched = BitSet()
3632

37-
esReduced = setdiff(es, eSolvedFixed)
38-
# println(" es = ", es, ", eSolvedFixed = ", eSolvedFixed, ", esReduced = ", esReduced)
39-
# println(" vs = ", vs, ", vSolvedFixed = ", vSolvedFixed)
40-
for eq in esReduced # iterate only over equations that are not in eSolvedFixed
33+
for eq in es # iterate only over equations that are not in eSolvedFixed
4134
for vj in Gsolvable[eq]
42-
if !(vj in vMatched) && (vj in vActive)
43-
r = add_edge_checked!(ict, Iterators.filter(!=(vj), G.graph.fadjlist[eq]), vj) do G
35+
if G.matching[vj] === unassigned && (vj in vActive)
36+
r = add_edge_checked!(ict, Iterators.filter(!=(vj), 𝑠neighbors(G.graph, eq)), vj) do G
4437
G.matching[vj] = eq
45-
push!(vMatched, vj)
4638
end
4739
r && break
4840
end
4941
end
5042
end
5143

52-
vSolved = filter(in(vMatched), topological_sort(ict))
44+
vSolved = filter(v->G.matching[v] !== unassigned, topological_sort(ict))
5345
inv_matching = Union{Missing, Int}[missing for _ = 1:nv(G)]
5446
for (v, eq) in pairs(G.matching)
5547
eq === unassigned && continue

src/systems/alias_elimination.jl

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,8 +2,6 @@ using SymbolicUtils: Rewriters
22

33
const KEEP = typemin(Int)
44

5-
include("compat/bareiss.jl")
6-
75
function alias_elimination(sys)
86
sys = initialize_system_structure(sys; quick_cancel=true)
97
s = structure(sys)

0 commit comments

Comments
 (0)