summaryrefslogtreecommitdiff
path: root/apps/plugins/puzzles/src/maxflow.c
diff options
context:
space:
mode:
Diffstat (limited to 'apps/plugins/puzzles/src/maxflow.c')
-rw-r--r--apps/plugins/puzzles/src/maxflow.c461
1 files changed, 461 insertions, 0 deletions
diff --git a/apps/plugins/puzzles/src/maxflow.c b/apps/plugins/puzzles/src/maxflow.c
new file mode 100644
index 0000000..028946b
--- /dev/null
+++ b/apps/plugins/puzzles/src/maxflow.c
@@ -0,0 +1,461 @@
+/*
+ * 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
generated by cgit v1.1 at 2026-04-25 13:42:42 +0000