| 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 |