ref: b3243d75043cf1d70beb88d2a36eaebfe85c2c3f
dir: /maxflow.c/
/*
* Edmonds-Karp algorithm for finding a maximum flow and minimum
* cut in a network. Almost identical to the Ford-Fulkerson
* algorithm, but apparently using breadth-first search to find the
* _shortest_ augmenting path is a good way to guarantee
* termination and ensure the time complexity is not dependent on
* the actual value of the maximum flow. I don't understand why
* that should be, but it's claimed on the Internet that it's been
* proved, and that's good enough for me. I prefer BFS to DFS
* anyway :-)
*/
#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include "maxflow.h"
#include "puzzles.h" /* for snewn/sfree */
int maxflow_with_scratch(void *scratch, int nv, int source, int sink,
int ne, const int *edges, const int *backedges,
const int *capacity, int *flow, int *cut)
{
int *todo = (int *)scratch;
int *prev = todo + nv;
int *firstedge = todo + 2*nv;
int *firstbackedge = todo + 3*nv;
int i, j, head, tail, from, to;
int totalflow;
/*
* Scan the edges array to find the index of the first edge
* from each node.
*/
j = 0;
for (i = 0; i < ne; i++)
while (j <= edges[2*i])
firstedge[j++] = i;
while (j < nv)
firstedge[j++] = ne;
assert(j == nv);
/*
* Scan the backedges array to find the index of the first edge
* _to_ each node.
*/
j = 0;
for (i = 0; i < ne; i++)
while (j <= edges[2*backedges[i]+1])
firstbackedge[j++] = i;
while (j < nv)
firstbackedge[j++] = ne;
assert(j == nv);
/*
* Start the flow off at zero on every edge.
*/
for (i = 0; i < ne; i++)
flow[i] = 0;
totalflow = 0;
/*
* Repeatedly look for an augmenting path, and follow it.
*/
while (1) {
/*
* Set up the prev array.
*/
for (i = 0; i < nv; i++)
prev[i] = -1;
/*
* Initialise the to-do list for BFS.
*/
head = tail = 0;
todo[tail++] = source;
/*
* Now do the BFS loop.
*/
while (head < tail && prev[sink] <= 0) {
from = todo[head++];
/*
* Try all the forward edges out of node `from'. For a
* forward edge to be valid, it must have flow
* currently less than its capacity.
*/
for (i = firstedge[from]; i < ne && edges[2*i] == from; i++) {
to = edges[2*i+1];
if (to == source || prev[to] >= 0)
continue;
if (capacity[i] >= 0 && flow[i] >= capacity[i])
continue;
/*
* This is a valid augmenting edge. Visit node `to'.
*/
prev[to] = 2*i;
todo[tail++] = to;
}
/*
* Try all the backward edges into node `from'. For a
* backward edge to be valid, it must have flow
* currently greater than zero.
*/
for (i = firstbackedge[from];
j = backedges[i], i < ne && edges[2*j+1]==from; i++) {
to = edges[2*j];
if (to == source || prev[to] >= 0)
continue;
if (flow[j] <= 0)
continue;
/*
* This is a valid augmenting edge. Visit node `to'.
*/
prev[to] = 2*j+1;
todo[tail++] = to;
}
}
/*
* If prev[sink] is non-null, we have found an augmenting
* path.
*/
if (prev[sink] >= 0) {
int max;
/*
* Work backwards along the path figuring out the
* maximum flow we can add.
*/
to = sink;
max = -1;
while (to != source) {
int spare;
/*
* Find the edge we're currently moving along.
*/
i = prev[to];
from = edges[i];
assert(from != to);
/*
* Determine the spare capacity of this edge.
*/
if (i & 1)
spare = flow[i / 2]; /* backward edge */
else if (capacity[i / 2] >= 0)
spare = capacity[i / 2] - flow[i / 2]; /* forward edge */
else
spare = -1; /* unlimited forward edge */
assert(spare != 0);
if (max < 0 || (spare >= 0 && spare < max))
max = spare;
to = from;
}
/*
* Fail an assertion if max is still < 0, i.e. there is
* an entirely unlimited path from source to sink. Also
* max should not _be_ zero, because by construction
* this _should_ be an augmenting path.
*/
assert(max > 0);
/*
* Now work backwards along the path again, this time
* actually adjusting the flow.
*/
to = sink;
while (to != source) {
/*
* Find the edge we're currently moving along.
*/
i = prev[to];
from = edges[i];
assert(from != to);
/*
* Adjust the edge.
*/
if (i & 1)
flow[i / 2] -= max; /* backward edge */
else
flow[i / 2] += max; /* forward edge */
to = from;
}
/*
* And adjust the overall flow counter.
*/
totalflow += max;
continue;
}
/*
* If we reach here, we have failed to find an augmenting
* path, which means we're done. Output the `cut' array if
* required, and leave.
*/
if (cut) {
for (i = 0; i < nv; i++) {
if (i == source || prev[i] >= 0)
cut[i] = 0;
else
cut[i] = 1;
}
}
return totalflow;
}
}
int maxflow_scratch_size(int nv)
{
return (nv * 4) * sizeof(int);
}
void maxflow_setup_backedges(int ne, const int *edges, int *backedges)
{
int i, n;
for (i = 0; i < ne; i++)
backedges[i] = i;
/*
* We actually can't use the C qsort() function, because we'd
* need to pass `edges' as a context parameter to its
* comparator function. So instead I'm forced to implement my
* own sorting algorithm internally, which is a pest. I'll use
* heapsort, because I like it.
*/
#define LESS(i,j) ( (edges[2*(i)+1] < edges[2*(j)+1]) || \
(edges[2*(i)+1] == edges[2*(j)+1] && \
edges[2*(i)] < edges[2*(j)]) )
#define PARENT(n) ( ((n)-1)/2 )
#define LCHILD(n) ( 2*(n)+1 )
#define RCHILD(n) ( 2*(n)+2 )
#define SWAP(i,j) do { int swaptmp = (i); (i) = (j); (j) = swaptmp; } while (0)
/*
* Phase 1: build the heap. We want the _largest_ element at
* the top.
*/
n = 0;
while (n < ne) {
n++;
/*
* Swap element n with its parent repeatedly to preserve
* the heap property.
*/
i = n-1;
while (i > 0) {
int p = PARENT(i);
if (LESS(backedges[p], backedges[i])) {
SWAP(backedges[p], backedges[i]);
i = p;
} else
break;
}
}
/*
* Phase 2: repeatedly remove the largest element and stick it
* at the top of the array.
*/
while (n > 0) {
/*
* The largest element is at position 0. Put it at the top,
* and swap the arbitrary element from that position into
* position 0.
*/
n--;
SWAP(backedges[0], backedges[n]);
/*
* Now repeatedly move that arbitrary element down the heap
* by swapping it with the more suitable of its children.
*/
i = 0;
while (1) {
int lc, rc;
lc = LCHILD(i);
rc = RCHILD(i);
if (lc >= n)
break; /* we've hit bottom */
if (rc >= n) {
/*
* Special case: there is only one child to check.
*/
if (LESS(backedges[i], backedges[lc]))
SWAP(backedges[i], backedges[lc]);
/* _Now_ we've hit bottom. */
break;
} else {
/*
* The common case: there are two children and we
* must check them both.
*/
if (LESS(backedges[i], backedges[lc]) ||
LESS(backedges[i], backedges[rc])) {
/*
* Pick the more appropriate child to swap with
* (i.e. the one which would want to be the
* parent if one were above the other - as one
* is about to be).
*/
if (LESS(backedges[lc], backedges[rc])) {
SWAP(backedges[i], backedges[rc]);
i = rc;
} else {
SWAP(backedges[i], backedges[lc]);
i = lc;
}
} else {
/* This element is in the right place; we're done. */
break;
}
}
}
}
#undef LESS
#undef PARENT
#undef LCHILD
#undef RCHILD
#undef SWAP
}
int maxflow(int nv, int source, int sink,
int ne, const int *edges, const int *capacity,
int *flow, int *cut)
{
void *scratch;
int *backedges;
int size;
int ret;
/*
* Allocate the space.
*/
size = ne * sizeof(int) + maxflow_scratch_size(nv);
backedges = smalloc(size);
if (!backedges)
return -1;
scratch = backedges + ne;
/*
* Set up the backedges array.
*/
maxflow_setup_backedges(ne, edges, backedges);
/*
* Call the main function.
*/
ret = maxflow_with_scratch(scratch, nv, source, sink, ne, edges,
backedges, capacity, flow, cut);
/*
* Free the scratch space.
*/
sfree(backedges);
/*
* And we're done.
*/
return ret;
}
#ifdef TESTMODE
#define MAXEDGES 256
#define MAXVERTICES 128
#define ADDEDGE(i,j) do{edges[ne*2] = (i); edges[ne*2+1] = (j); ne++;}while(0)
int compare_edge(const void *av, const void *bv)
{
const int *a = (const int *)av;
const int *b = (const int *)bv;
if (a[0] < b[0])
return -1;
else if (a[0] > b[0])
return +1;
else if (a[1] < b[1])
return -1;
else if (a[1] > b[1])
return +1;
else
return 0;
}
int main(void)
{
int edges[MAXEDGES*2], ne, nv;
int capacity[MAXEDGES], flow[MAXEDGES], cut[MAXVERTICES];
int source, sink, p, q, i, j, ret;
/*
* Use this algorithm to find a maximal complete matching in a
* bipartite graph.
*/
ne = 0;
nv = 0;
source = nv++;
p = nv;
nv += 5;
q = nv;
nv += 5;
sink = nv++;
for (i = 0; i < 5; i++) {
capacity[ne] = 1;
ADDEDGE(source, p+i);
}
for (i = 0; i < 5; i++) {
capacity[ne] = 1;
ADDEDGE(q+i, sink);
}
j = ne;
capacity[ne] = 1; ADDEDGE(p+0,q+0);
capacity[ne] = 1; ADDEDGE(p+1,q+0);
capacity[ne] = 1; ADDEDGE(p+1,q+1);
capacity[ne] = 1; ADDEDGE(p+2,q+1);
capacity[ne] = 1; ADDEDGE(p+2,q+2);
capacity[ne] = 1; ADDEDGE(p+3,q+2);
capacity[ne] = 1; ADDEDGE(p+3,q+3);
capacity[ne] = 1; ADDEDGE(p+4,q+3);
/* capacity[ne] = 1; ADDEDGE(p+2,q+4); */
qsort(edges, ne, 2*sizeof(int), compare_edge);
ret = maxflow(nv, source, sink, ne, edges, capacity, flow, cut);
printf("ret = %d\n", ret);
for (i = 0; i < ne; i++)
printf("flow %d: %d -> %d\n", flow[i], edges[2*i], edges[2*i+1]);
for (i = 0; i < nv; i++)
if (cut[i] == 0)
printf("difficult set includes %d\n", i);
return 0;
}
#endif