+++ /dev/null
-/*
- * 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
+++ /dev/null
-/*
- * 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 :-)
- */
-
-#ifndef MAXFLOW_MAXFLOW_H
-#define MAXFLOW_MAXFLOW_H
-
-/*
- * The actual algorithm.
- *
- * Inputs:
- *
- * - `scratch' is previously allocated scratch space of a size
- * previously determined by calling `maxflow_scratch_size'.
- *
- * - `nv' is the number of vertices. Vertices are assumed to be
- * numbered from 0 to nv-1.
- *
- * - `source' and `sink' are the distinguished source and sink
- * vertices.
- *
- * - `ne' is the number of edges in the graph.
- *
- * - `edges' is an array of 2*ne integers, giving a (source, dest)
- * pair for each network edge. Edge pairs are expected to be
- * sorted in lexicographic order.
- *
- * - `backedges' is an array of `ne' integers, each a distinct
- * index into `edges'. The edges in `edges', if permuted as
- * specified by this array, should end up sorted in the _other_
- * lexicographic order, i.e. dest taking priority over source.
- *
- * - `capacity' is an array of `ne' integers, giving a maximum
- * flow capacity for each edge. A negative value is taken to
- * indicate unlimited capacity on that edge, but note that there
- * may not be any unlimited-capacity _path_ from source to sink
- * or an assertion will be failed.
- *
- * Output:
- *
- * - `flow' must be non-NULL. It is an array of `ne' integers,
- * each giving the final flow along each edge.
- *
- * - `cut' may be NULL. If non-NULL, it is an array of `nv'
- * integers, which will be set to zero or one on output, in such
- * a way that:
- * + the set of zero vertices includes the source
- * + the set of one vertices includes the sink
- * + the maximum flow capacity between the zero and one vertex
- * sets is achieved (i.e. all edges from a zero vertex to a
- * one vertex are at full capacity, while all edges from a
- * one vertex to a zero vertex have no flow at all).
- *
- * - the returned value from the function is the total flow
- * achieved.
- */
-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);
-
-/*
- * The above function expects its `scratch' and `backedges'
- * parameters to have already been set up. This allows you to set
- * them up once and use them in multiple invocates of the
- * algorithm. Now I provide functions to actually do the setting
- * up.
- */
-int maxflow_scratch_size(int nv);
-void maxflow_setup_backedges(int ne, const int *edges, int *backedges);
-
-/*
- * Simplified version of the above function. All parameters are the
- * same, except that `scratch' and `backedges' are constructed
- * internally. This is the simplest way to call the algorithm as a
- * one-off; however, if you need to call it multiple times on the
- * same network, it is probably better to call the above version
- * directly so that you only construct `scratch' and `backedges'
- * once.
- *
- * Additional return value is now -1, meaning that scratch space
- * could not be allocated.
- */
-int maxflow(int nv, int source, int sink,
- int ne, const int *edges, const int *capacity,
- int *flow, int *cut);
-
-#endif /* MAXFLOW_MAXFLOW_H */
~ian / sgt-puzzles.git / commitdiff