New upstream release
[sgt-puzzles/ydirson.git] / maxflow.c
blob028946b9bd400b10eef95a76b4bd6c6c6e377ad6
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 :-)
13 #include <assert.h>
14 #include <stdlib.h>
15 #include <stdio.h>
17 #include "maxflow.h"
19 #include "puzzles.h" /* for snewn/sfree */
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)
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;
33 * Scan the edges array to find the index of the first edge
34 * from each node.
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);
45 * Scan the backedges array to find the index of the first edge
46 * _to_ each node.
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);
57 * Start the flow off at zero on every edge.
59 for (i = 0; i < ne; i++)
60 flow[i] = 0;
61 totalflow = 0;
64 * Repeatedly look for an augmenting path, and follow it.
66 while (1) {
69 * Set up the prev array.
71 for (i = 0; i < nv; i++)
72 prev[i] = -1;
75 * Initialise the to-do list for BFS.
77 head = tail = 0;
78 todo[tail++] = source;
81 * Now do the BFS loop.
83 while (head < tail && prev[sink] <= 0) {
84 from = todo[head++];
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.
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;
98 * This is a valid augmenting edge. Visit node `to'.
100 prev[to] = 2*i;
101 todo[tail++] = to;
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.
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;
117 * This is a valid augmenting edge. Visit node `to'.
119 prev[to] = 2*j+1;
120 todo[tail++] = to;
125 * If prev[sink] is non-null, we have found an augmenting
126 * path.
128 if (prev[sink] >= 0) {
129 int max;
132 * Work backwards along the path figuring out the
133 * maximum flow we can add.
135 to = sink;
136 max = -1;
137 while (to != source) {
138 int spare;
141 * Find the edge we're currently moving along.
143 i = prev[to];
144 from = edges[i];
145 assert(from != to);
148 * Determine the spare capacity of this edge.
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 */
157 assert(spare != 0);
159 if (max < 0 || (spare >= 0 && spare < max))
160 max = spare;
162 to = from;
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.
170 assert(max > 0);
173 * Now work backwards along the path again, this time
174 * actually adjusting the flow.
176 to = sink;
177 while (to != source) {
179 * Find the edge we're currently moving along.
181 i = prev[to];
182 from = edges[i];
183 assert(from != to);
186 * Adjust the edge.
188 if (i & 1)
189 flow[i / 2] -= max; /* backward edge */
190 else
191 flow[i / 2] += max; /* forward edge */
193 to = from;
197 * And adjust the overall flow counter.
199 totalflow += max;
201 continue;
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.
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;
217 return totalflow;
221 int maxflow_scratch_size(int nv)
223 return (nv * 4) * sizeof(int);
226 void maxflow_setup_backedges(int ne, const int *edges, int *backedges)
228 int i, n;
230 for (i = 0; i < ne; i++)
231 backedges[i] = i;
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.
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)
250 * Phase 1: build the heap. We want the _largest_ element at
251 * the top.
253 n = 0;
254 while (n < ne) {
255 n++;
258 * Swap element n with its parent repeatedly to preserve
259 * the heap property.
261 i = n-1;
263 while (i > 0) {
264 int p = PARENT(i);
266 if (LESS(backedges[p], backedges[i])) {
267 SWAP(backedges[p], backedges[i]);
268 i = p;
269 } else
270 break;
275 * Phase 2: repeatedly remove the largest element and stick it
276 * at the top of the array.
278 while (n > 0) {
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.
284 n--;
285 SWAP(backedges[0], backedges[n]);
288 * Now repeatedly move that arbitrary element down the heap
289 * by swapping it with the more suitable of its children.
291 i = 0;
292 while (1) {
293 int lc, rc;
295 lc = LCHILD(i);
296 rc = RCHILD(i);
298 if (lc >= n)
299 break; /* we've hit bottom */
301 if (rc >= n) {
303 * Special case: there is only one child to check.
305 if (LESS(backedges[i], backedges[lc]))
306 SWAP(backedges[i], backedges[lc]);
308 /* _Now_ we've hit bottom. */
309 break;
310 } else {
312 * The common case: there are two children and we
313 * must check them both.
315 if (LESS(backedges[i], backedges[lc]) ||
316 LESS(backedges[i], backedges[rc])) {
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).
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;
330 } else {
331 /* This element is in the right place; we're done. */
332 break;
338 #undef LESS
339 #undef PARENT
340 #undef LCHILD
341 #undef RCHILD
342 #undef SWAP
346 int maxflow(int nv, int source, int sink,
347 int ne, const int *edges, const int *capacity,
348 int *flow, int *cut)
350 void *scratch;
351 int *backedges;
352 int size;
353 int ret;
356 * Allocate the space.
358 size = ne * sizeof(int) + maxflow_scratch_size(nv);
359 backedges = smalloc(size);
360 if (!backedges)
361 return -1;
362 scratch = backedges + ne;
365 * Set up the backedges array.
367 maxflow_setup_backedges(ne, edges, backedges);
370 * Call the main function.
372 ret = maxflow_with_scratch(scratch, nv, source, sink, ne, edges,
373 backedges, capacity, flow, cut);
376 * Free the scratch space.
378 sfree(backedges);
381 * And we're done.
383 return ret;
386 #ifdef TESTMODE
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)
392 int compare_edge(const void *av, const void *bv)
394 const int *a = (const int *)av;
395 const int *b = (const int *)bv;
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;
409 int main(void)
411 int edges[MAXEDGES*2], ne, nv;
412 int capacity[MAXEDGES], flow[MAXEDGES], cut[MAXVERTICES];
413 int source, sink, p, q, i, j, ret;
416 * Use this algorithm to find a maximal complete matching in a
417 * bipartite graph.
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);
431 for (i = 0; i < 5; i++) {
432 capacity[ne] = 1;
433 ADDEDGE(q+i, sink);
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);
447 ret = maxflow(nv, source, sink, ne, edges, capacity, flow, cut);
449 printf("ret = %d\n", ret);
451 for (i = 0; i < ne; i++)
452 printf("flow %d: %d -> %d\n", flow[i], edges[2*i], edges[2*i+1]);
454 for (i = 0; i < nv; i++)
455 if (cut[i] == 0)
456 printf("difficult set includes %d\n", i);
458 return 0;
461 #endif