maxflow.c

   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