-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgrAlg.nim
187 lines (172 loc) · 8.64 KB
/
grAlg.nim
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
## Various algorithms on abstract directed graphs (digraphs, `dg`). Beyond an
## integral type for node Ids, they only need a span `n` of such & iterators for
## `nodes(dg)` & `dests(dg, node)`. Ids need not be dense on `0..<n`, but algos
## do allocate `n`-proportional space. TODO Maybe make routines be abstract Re:
## nodes containing needed metadata and ref/ptr node Ids?
#TODO Max Flow? Other?
import std/[deques, algorithm, tables], gaPrioQ, cligen/sysUt; export tables
template reverse*(dg, nodes, dests, result): untyped =
## Add reverse arcs of `dg` to `result` (must have `mgetOrPut(Id, seq[Id])`).
for x in nodes(dg):
for y in dests(dg, x):
result.mgetOrPut(y, @[]).add x
const gaMaxPath* {.intdefine.} = 10_000 # Diams even this big are very rare
template trace(pred, b, e, result) =
var p = e # Trace path from `pred[e]` links.
while p != b and result.len < gaMaxPath:
result.add p; p = pred[p]
result.add b
result.reverse
template shortestPathBFS*[I](dg; n: int; b,e: I; dests, did, pred, q): untyped =
## Breadth First Search For Shortest b -> e path; User provided mem buffers.
var result: seq[I] # Initially all nodes have did[i]=false
did.setLen n; zeroMem did[0].addr, n*did[0].sizeof
pred.setLen n # Predecessors in first found path
clear q; addLast q, b # Initialize q with just `b`.
did[b] = true # Mark first node done
var done = false # Loop is unfinished
while not done and q.len > 0: # Standard BFS algorithm
let x = popFirst q # Need FIFO order to get SHORTEST path
for y in dests(dg, x): #..Other orders halt @achievable paths.
if not did[y]:
did[y] = true
pred[y] = x
if y == e: # FOUND: trace reverse path & return
trace pred, b, e, result
done = true; break # Could be `return` in a `proc`.
else:
addLast q, y
result
template shortestPathBFS*[I](dg; n: int; b,e: I; dests): untyped =
## Breadth First Search For Shortest b -> e path; Uses new mem buffers.
var did = newSeqNoInit[bool](n) # Already zeroed in mem reusing entry pt
var pred = newSeqNoInit[I](n)
var q = initDeque[I](32) # Nodes to check
shortestPathBFS(dg, n, b, e, dests, did, pred, q)
template shortestPathPFS*[I](dg; n: int; b,e: I; nodes, dests): untyped=
## Dijkstra Min Cost Path Algorithm for b -> e; Unlike most other algos here,
## `dests` must be compatible with `for (dest, cost: C) in dests(dg, n): ..`.
## As with all Dijkstra, length/costs must be > 0 (but can be `float`).
type C = typeof(for s in nodes(dg): (for (d,c) in dests(dg, s): c))
var result: seq[I]
var cost = newSeq[C](n) # This uses about 12*n space
var pred = newSeqNoInit[I](n) # Dijkstra Min Cost Path
var idx = newSeq[I](n) # map[x] == heap index
if true: # Scope for `proc`
proc iSet(k: I, i: int) = idx[k.int] = I(i + 1) # 0 encodes MISSING
var q: PrioQ[C, I] # Another 8*n space
for i in nodes(dg):
cost[i] = (if i == b: C(0) else: C.high)
idx[i] = i + 1
push q, cost[i], i, iSet
while q.len > 0:
let (xC, x) = q.pop(iSet)
idx[x] = 0 # Mark completed early to not re-do
if xC != C.high: # reachable
for (y, yC) in dests(dg, x):
if idx[y] != 0: # shortest remains undetermined
let alt = cost[x] + yC
if alt < cost[y]:
cost[y] = alt
pred[y] = x
edit q, alt, idx[y] - 1, iSet
trace pred, b, e, result # trace path into result
result
proc udcRoot*(up: var seq[int], x: int): int {.inline.} =
## Find root in Union-Find U)nD)irected C)omponent algo w/path compression.
result = x # Find root defined by parent == self
while up[result] != result:
result = up[result]
var x = x # Compress path afterwards
while up[x] != result: #..by doing up[all x in path] <- result
let up0 = up[x]; up[x] = result; x = up0
proc udcJoin*(up, sz: var seq[int]; x, y: int) {.inline.} =
## U)nD)irectedC)omponent Join (aka Union) in a Union-Find algorithm.
let x = udcRoot(up, x) # Join/union by size
let y = udcRoot(up, y)
if y != x: # x & y are not already joined
if sz[x] < sz[y]: # Attach smaller tree..
up[x] = y #..to root of larger
sz[y] += sz[x] # and update size of larger
else: # Mirror case of above
up[y] = x
sz[x] += sz[y]
template unDirComponIdSz*(dg; n: int; nodes, dests): untyped =
## Evals to an `up` componId array for each node Id in digraph `dg` treated as
## a reciprocated/undirected graph.
var up = newSeq[int](n) # nodeId -> parent id
var sz = newSeq[int](n) # nodeId -> sz
for x in nodes(dg): # Initialize:
up[x] = x; sz[x] = 1 # parents all self & sizes all 1
for x in nodes(dg): # Incorp edges via union-find/join-root
for y in dests(dg, x): udcJoin up, sz, x, y
(up, sz) # Now udcRoot(up, x) == componentId
template unDirCompons*(dg; n: int; nodes, dests): untyped =
## Evals to a `Table[cId, seq[nodeId]]` of connected components
var (up, sz) = unDirComponIdSz(dg, n, nodes, dests)
var t = initTable[int, seq[int]](n) # Compon id => nodeId list
for x in nodes(dg): # Add nodeIds to cId keys
t.mgetOrPut(udcRoot(up, int(x)), newSeqOfCap[int](sz[x])).add x
t
template unDirComponSizes*(dg; n: int; nodes, dests): untyped =
## Evals to a `CountTable[int]` of component sizes.
var (up, _ {.used.}) = unDirComponIdSz(dg, n, nodes, dests)
var t: CountTable[int]
for x in nodes(dg): t.inc udcRoot(up, int(x))
t
template minSpanTree*(dg; n: int; nodes, dests): untyped =
## Evals to a Min Cost Spanning Tree via Kruskal's Algorithm. `dests` here is
## a cost/weighted-one like `shortestPathPFS`.
type I = typeof(for s in nodes(dg): s)
type C = typeof(for s in nodes(dg): (for (d,c) in dests(dg, s): c))
var result, arcs: seq[tuple[cost: C; src, dst: I]]
var up = newSeq[int](n) # nodeId -> parent id
var sz = newSeq[int](n) # nodeId -> sz
var nNode = 0
for src in nodes(dg): # Initialize:
inc nNode
up[src] = src; sz[src] = 1 # parents all self & sizes all 1
for (dst, cost) in dests(dg, src):
arcs.add (cost, src, dst) # accumulate arcs
arcs.sort # Conditionally nsort?
for arc in arcs:
if result.len == nNode: break
if udcRoot(up, arc.src) != udcRoot(up, arc.dst):
udcJoin up, sz, arc.src, arc.dst
result.add arc
result
template transClosure*[I](dg; n: int; r: I; nodes, dests): untyped =
## Evals to a `seq[I]` containing the transitive closure of root `r` on `dg`.
var did = newSeq[bool](n) # Could save space with `result: HashSet[I]` or
var result: seq[I] #..just `seq[T].contains` for really small TCs.
if true: # scope for `proc`
proc visit(x: I) = # Depth First Search (DFS)
result.add x; did[x] = true
for y in dests(dg, x):
if not did[y]: visit(y)
visit(r)
result
template topoSort*[I](dg; n: int; nodes, dests): untyped =
## Eval to a topological sort of the graph, or `len==0 seq` if cyclic.
var did = newSeq[bool](n) # Could save space with `result: HashSet[I]` or
var result: seq[I] #..just `seq[T].contains` for really small TCs.
if true: # scope for `procs`
proc visit(x: I) = # Depth First Search (DFS)
did[x] = true # Mark current node as visited.
for y in dests(dg, x): # Recurse to all kids of this node
if not did[y]: visit(y)
result.add x # Add node to stack storing reversed result
for b in nodes(dg):
if not did[b]: visit(b)
result.reverse
result
template isCyclic*[I](dg; n: int; nodes, dests; tsort: seq[I]): untyped =
## Eval to true if graph w/pre-computed topoSort `tsort` has any cycles.
var ix = newSeqNoInit[I](n) # Position of node in topological sort
var result = false
for j in 0 ..< tsort.len: ix[tsort[j].int] = I(j)
block both:
for x in nodes(dg):
for y in dests(dg, x):
if ix[x] > ix[y]: result = true; break both
result