86e60e3d |
1 | /* |
2 | * Edmonds-Karp algorithm for finding a maximum flow and minimum |
3 | * cut in a network. Almost identical to the Ford-Fulkerson |
4 | * algorithm, but apparently using breadth-first search to find the |
5 | * _shortest_ augmenting path is a good way to guarantee |
6 | * termination and ensure the time complexity is not dependent on |
7 | * the actual value of the maximum flow. I don't understand why |
8 | * that should be, but it's claimed on the Internet that it's been |
9 | * proved, and that's good enough for me. I prefer BFS to DFS |
10 | * anyway :-) |
11 | */ |
12 | |
13 | #include <assert.h> |
14 | #include <stdlib.h> |
15 | #include <stdio.h> |
16 | |
17 | #include "maxflow.h" |
18 | |
19 | #include "puzzles.h" /* for snewn/sfree */ |
20 | |
21 | int maxflow_with_scratch(void *scratch, int nv, int source, int sink, |
22 | int ne, const int *edges, const int *backedges, |
23 | const int *capacity, int *flow, int *cut) |
24 | { |
25 | int *todo = (int *)scratch; |
26 | int *prev = todo + nv; |
27 | int *firstedge = todo + 2*nv; |
28 | int *firstbackedge = todo + 3*nv; |
29 | int i, j, head, tail, from, to; |
30 | int totalflow; |
31 | |
32 | /* |
33 | * Scan the edges array to find the index of the first edge |
34 | * from each node. |
35 | */ |
36 | j = 0; |
37 | for (i = 0; i < ne; i++) |
38 | while (j <= edges[2*i]) |
39 | firstedge[j++] = i; |
40 | while (j < nv) |
41 | firstedge[j++] = ne; |
42 | assert(j == nv); |
43 | |
44 | /* |
45 | * Scan the backedges array to find the index of the first edge |
46 | * _to_ each node. |
47 | */ |
48 | j = 0; |
49 | for (i = 0; i < ne; i++) |
50 | while (j <= edges[2*backedges[i]+1]) |
51 | firstbackedge[j++] = i; |
52 | while (j < nv) |
53 | firstbackedge[j++] = ne; |
54 | assert(j == nv); |
55 | |
56 | /* |
57 | * Start the flow off at zero on every edge. |
58 | */ |
59 | for (i = 0; i < ne; i++) |
60 | flow[i] = 0; |
61 | totalflow = 0; |
62 | |
63 | /* |
64 | * Repeatedly look for an augmenting path, and follow it. |
65 | */ |
66 | while (1) { |
67 | |
68 | /* |
69 | * Set up the prev array. |
70 | */ |
71 | for (i = 0; i < nv; i++) |
72 | prev[i] = -1; |
73 | |
74 | /* |
75 | * Initialise the to-do list for BFS. |
76 | */ |
77 | head = tail = 0; |
78 | todo[tail++] = source; |
79 | |
80 | /* |
81 | * Now do the BFS loop. |
82 | */ |
83 | while (head < tail && prev[sink] <= 0) { |
84 | from = todo[head++]; |
85 | |
86 | /* |
87 | * Try all the forward edges out of node `from'. For a |
88 | * forward edge to be valid, it must have flow |
89 | * currently less than its capacity. |
90 | */ |
91 | for (i = firstedge[from]; i < ne && edges[2*i] == from; i++) { |
92 | to = edges[2*i+1]; |
93 | if (to == source || prev[to] >= 0) |
94 | continue; |
95 | if (capacity[i] >= 0 && flow[i] >= capacity[i]) |
96 | continue; |
97 | /* |
98 | * This is a valid augmenting edge. Visit node `to'. |
99 | */ |
100 | prev[to] = 2*i; |
101 | todo[tail++] = to; |
102 | } |
103 | |
104 | /* |
105 | * Try all the backward edges into node `from'. For a |
106 | * backward edge to be valid, it must have flow |
107 | * currently greater than zero. |
108 | */ |
109 | for (i = firstbackedge[from]; |
110 | j = backedges[i], i < ne && edges[2*j+1]==from; i++) { |
111 | to = edges[2*j]; |
112 | if (to == source || prev[to] >= 0) |
113 | continue; |
114 | if (flow[j] <= 0) |
115 | continue; |
116 | /* |
117 | * This is a valid augmenting edge. Visit node `to'. |
118 | */ |
119 | prev[to] = 2*j+1; |
120 | todo[tail++] = to; |
121 | } |
122 | } |
123 | |
124 | /* |
125 | * If prev[sink] is non-null, we have found an augmenting |
126 | * path. |
127 | */ |
128 | if (prev[sink] >= 0) { |
129 | int max; |
130 | |
131 | /* |
132 | * Work backwards along the path figuring out the |
133 | * maximum flow we can add. |
134 | */ |
135 | to = sink; |
136 | max = -1; |
137 | while (to != source) { |
138 | int spare; |
139 | |
140 | /* |
141 | * Find the edge we're currently moving along. |
142 | */ |
143 | i = prev[to]; |
144 | from = edges[i]; |
145 | assert(from != to); |
146 | |
147 | /* |
148 | * Determine the spare capacity of this edge. |
149 | */ |
150 | if (i & 1) |
151 | spare = flow[i / 2]; /* backward edge */ |
152 | else if (capacity[i / 2] >= 0) |
153 | spare = capacity[i / 2] - flow[i / 2]; /* forward edge */ |
154 | else |
155 | spare = -1; /* unlimited forward edge */ |
156 | |
157 | assert(spare != 0); |
158 | |
159 | if (max < 0 || (spare >= 0 && spare < max)) |
160 | max = spare; |
161 | |
162 | to = from; |
163 | } |
164 | /* |
165 | * Fail an assertion if max is still < 0, i.e. there is |
166 | * an entirely unlimited path from source to sink. Also |
167 | * max should not _be_ zero, because by construction |
168 | * this _should_ be an augmenting path. |
169 | */ |
170 | assert(max > 0); |
171 | |
172 | /* |
173 | * Now work backwards along the path again, this time |
174 | * actually adjusting the flow. |
175 | */ |
176 | to = sink; |
177 | while (to != source) { |
178 | /* |
179 | * Find the edge we're currently moving along. |
180 | */ |
181 | i = prev[to]; |
182 | from = edges[i]; |
183 | assert(from != to); |
184 | |
185 | /* |
186 | * Adjust the edge. |
187 | */ |
188 | if (i & 1) |
189 | flow[i / 2] -= max; /* backward edge */ |
190 | else |
191 | flow[i / 2] += max; /* forward edge */ |
192 | |
193 | to = from; |
194 | } |
195 | |
196 | /* |
197 | * And adjust the overall flow counter. |
198 | */ |
199 | totalflow += max; |
200 | |
201 | continue; |
202 | } |
203 | |
204 | /* |
205 | * If we reach here, we have failed to find an augmenting |
206 | * path, which means we're done. Output the `cut' array if |
207 | * required, and leave. |
208 | */ |
209 | if (cut) { |
210 | for (i = 0; i < nv; i++) { |
211 | if (i == source || prev[i] >= 0) |
212 | cut[i] = 0; |
213 | else |
214 | cut[i] = 1; |
215 | } |
216 | } |
217 | return totalflow; |
218 | } |
219 | } |
220 | |
221 | int maxflow_scratch_size(int nv) |
222 | { |
223 | return (nv * 4) * sizeof(int); |
224 | } |
225 | |
226 | void maxflow_setup_backedges(int ne, const int *edges, int *backedges) |
227 | { |
228 | int i, n; |
229 | |
230 | for (i = 0; i < ne; i++) |
231 | backedges[i] = i; |
232 | |
233 | /* |
234 | * We actually can't use the C qsort() function, because we'd |
235 | * need to pass `edges' as a context parameter to its |
236 | * comparator function. So instead I'm forced to implement my |
237 | * own sorting algorithm internally, which is a pest. I'll use |
238 | * heapsort, because I like it. |
239 | */ |
240 | |
241 | #define LESS(i,j) ( (edges[2*(i)+1] < edges[2*(j)+1]) || \ |
242 | (edges[2*(i)+1] == edges[2*(j)+1] && \ |
243 | edges[2*(i)] < edges[2*(j)]) ) |
244 | #define PARENT(n) ( ((n)-1)/2 ) |
245 | #define LCHILD(n) ( 2*(n)+1 ) |
246 | #define RCHILD(n) ( 2*(n)+2 ) |
247 | #define SWAP(i,j) do { int swaptmp = (i); (i) = (j); (j) = swaptmp; } while (0) |
248 | |
249 | /* |
250 | * Phase 1: build the heap. We want the _largest_ element at |
251 | * the top. |
252 | */ |
253 | n = 0; |
254 | while (n < ne) { |
255 | n++; |
256 | |
257 | /* |
258 | * Swap element n with its parent repeatedly to preserve |
259 | * the heap property. |
260 | */ |
261 | i = n-1; |
262 | |
263 | while (i > 0) { |
264 | int p = PARENT(i); |
265 | |
266 | if (LESS(backedges[p], backedges[i])) { |
267 | SWAP(backedges[p], backedges[i]); |
268 | i = p; |
269 | } else |
270 | break; |
271 | } |
272 | } |
273 | |
274 | /* |
275 | * Phase 2: repeatedly remove the largest element and stick it |
276 | * at the top of the array. |
277 | */ |
278 | while (n > 0) { |
279 | /* |
280 | * The largest element is at position 0. Put it at the top, |
281 | * and swap the arbitrary element from that position into |
282 | * position 0. |
283 | */ |
284 | n--; |
285 | SWAP(backedges[0], backedges[n]); |
286 | |
287 | /* |
288 | * Now repeatedly move that arbitrary element down the heap |
289 | * by swapping it with the more suitable of its children. |
290 | */ |
291 | i = 0; |
292 | while (1) { |
293 | int lc, rc; |
294 | |
295 | lc = LCHILD(i); |
296 | rc = RCHILD(i); |
297 | |
298 | if (lc >= n) |
299 | break; /* we've hit bottom */ |
300 | |
301 | if (rc >= n) { |
302 | /* |
303 | * Special case: there is only one child to check. |
304 | */ |
305 | if (LESS(backedges[i], backedges[lc])) |
306 | SWAP(backedges[i], backedges[lc]); |
307 | |
308 | /* _Now_ we've hit bottom. */ |
309 | break; |
310 | } else { |
311 | /* |
312 | * The common case: there are two children and we |
313 | * must check them both. |
314 | */ |
315 | if (LESS(backedges[i], backedges[lc]) || |
316 | LESS(backedges[i], backedges[rc])) { |
317 | /* |
318 | * Pick the more appropriate child to swap with |
319 | * (i.e. the one which would want to be the |
320 | * parent if one were above the other - as one |
321 | * is about to be). |
322 | */ |
323 | if (LESS(backedges[lc], backedges[rc])) { |
324 | SWAP(backedges[i], backedges[rc]); |
325 | i = rc; |
326 | } else { |
327 | SWAP(backedges[i], backedges[lc]); |
328 | i = lc; |
329 | } |
330 | } else { |
331 | /* This element is in the right place; we're done. */ |
332 | break; |
333 | } |
334 | } |
335 | } |
336 | } |
337 | |
338 | #undef LESS |
339 | #undef PARENT |
340 | #undef LCHILD |
341 | #undef RCHILD |
342 | #undef SWAP |
343 | |
344 | } |
345 | |
346 | int maxflow(int nv, int source, int sink, |
347 | int ne, const int *edges, const int *capacity, |
348 | int *flow, int *cut) |
349 | { |
350 | void *scratch; |
351 | int *backedges; |
352 | int size; |
353 | int ret; |
354 | |
355 | /* |
356 | * Allocate the space. |
357 | */ |
358 | size = ne * sizeof(int) + maxflow_scratch_size(nv); |
359 | backedges = smalloc(size); |
360 | if (!backedges) |
361 | return -1; |
362 | scratch = backedges + ne; |
363 | |
364 | /* |
365 | * Set up the backedges array. |
366 | */ |
367 | maxflow_setup_backedges(ne, edges, backedges); |
368 | |
369 | /* |
370 | * Call the main function. |
371 | */ |
372 | ret = maxflow_with_scratch(scratch, nv, source, sink, ne, edges, |
373 | backedges, capacity, flow, cut); |
374 | |
375 | /* |
376 | * Free the scratch space. |
377 | */ |
378 | sfree(backedges); |
379 | |
380 | /* |
381 | * And we're done. |
382 | */ |
383 | return ret; |
384 | } |
385 | |
386 | #ifdef TESTMODE |
387 | |
388 | #define MAXEDGES 256 |
389 | #define MAXVERTICES 128 |
390 | #define ADDEDGE(i,j) do{edges[ne*2] = (i); edges[ne*2+1] = (j); ne++;}while(0) |
391 | |
392 | int compare_edge(const void *av, const void *bv) |
393 | { |
394 | const int *a = (const int *)av; |
395 | const int *b = (const int *)bv; |
396 | |
397 | if (a[0] < b[0]) |
398 | return -1; |
399 | else if (a[0] > b[0]) |
400 | return +1; |
401 | else if (a[1] < b[1]) |
402 | return -1; |
403 | else if (a[1] > b[1]) |
404 | return +1; |
405 | else |
406 | return 0; |
407 | } |
408 | |
409 | int main(void) |
410 | { |
411 | int edges[MAXEDGES*2], ne, nv; |
412 | int capacity[MAXEDGES], flow[MAXEDGES], cut[MAXVERTICES]; |
413 | int source, sink, p, q, i, j, ret; |
414 | |
415 | /* |
416 | * Use this algorithm to find a maximal complete matching in a |
417 | * bipartite graph. |
418 | */ |
419 | ne = 0; |
420 | nv = 0; |
421 | source = nv++; |
422 | p = nv; |
423 | nv += 5; |
424 | q = nv; |
425 | nv += 5; |
426 | sink = nv++; |
427 | for (i = 0; i < 5; i++) { |
428 | capacity[ne] = 1; |
429 | ADDEDGE(source, p+i); |
430 | } |
431 | for (i = 0; i < 5; i++) { |
432 | capacity[ne] = 1; |
433 | ADDEDGE(q+i, sink); |
434 | } |
435 | j = ne; |
436 | capacity[ne] = 1; ADDEDGE(p+0,q+0); |
437 | capacity[ne] = 1; ADDEDGE(p+1,q+0); |
438 | capacity[ne] = 1; ADDEDGE(p+1,q+1); |
439 | capacity[ne] = 1; ADDEDGE(p+2,q+1); |
440 | capacity[ne] = 1; ADDEDGE(p+2,q+2); |
441 | capacity[ne] = 1; ADDEDGE(p+3,q+2); |
442 | capacity[ne] = 1; ADDEDGE(p+3,q+3); |
443 | capacity[ne] = 1; ADDEDGE(p+4,q+3); |
444 | /* capacity[ne] = 1; ADDEDGE(p+2,q+4); */ |
445 | qsort(edges, ne, 2*sizeof(int), compare_edge); |
446 | |
447 | ret = maxflow(nv, source, sink, ne, edges, capacity, flow, cut); |
448 | |
449 | printf("ret = %d\n", ret); |
450 | |
451 | for (i = 0; i < ne; i++) |
452 | printf("flow %d: %d -> %d\n", flow[i], edges[2*i], edges[2*i+1]); |
453 | |
454 | for (i = 0; i < nv; i++) |
455 | if (cut[i] == 0) |
456 | printf("difficult set includes %d\n", i); |
457 | |
458 | return 0; |
459 | } |
460 | |
461 | #endif |