1 /* mincut.c */
2
3 /* Written by Andrew Makhorin <mao@gnu.org>, October 2015. */
4
5 #include <limits.h>
6
7 #include "maxflow.h"
8 #include "mincut.h"
9 #include "misc.h"
10
11 /***********************************************************************
12 * NAME
13 *
14 * min_cut - find min cut in undirected capacitated network
15 *
16 * SYNOPSIS
17 *
18 * #include "mincut.h"
19 * int min_cut(int nn, int ne, const int beg[], const int end[],
20 * const cap[], int cut[]);
21 *
22 * DESCRIPTION
23 *
24 * This routine finds min cut in a given undirected network.
25 *
26 * The undirected capacitated network is specified by the parameters
27 * nn, ne, beg, end, and cap. The parameter nn specifies the number of
28 * vertices (nodes), nn >= 2, and the parameter ne specifies the number
29 * of edges, ne >= 0. The network edges are specified by triplets
30 * (beg[k], end[k], cap[k]) for k = 1, ..., ne, where beg[k] < end[k]
31 * are numbers of the first and second nodes of k-th edge, and
32 * cap[k] > 0 is a capacity of k-th edge. Loops and multiple edges are
33 * not allowed.
34 *
35 * Let V be the set of nodes of the network and let W be an arbitrary
36 * non-empty proper subset of V. A cut associated with the subset W is
37 * a subset of all the edges, one node of which belongs to W and other
38 * node belongs to V \ W. The capacity of a cut (W, V \ W) is the sum
39 * of the capacities of all the edges, which belong to the cut. Minimal
40 * cut is a cut, whose capacity is minimal.
41 *
42 * On exit the routine stores flags of nodes v[i], i = 1, ..., nn, to
43 * locations cut[i], where cut[i] = 1 means that v[i] belongs to W and
44 * cut[i] = 0 means that v[i] belongs to V \ W, where W corresponds to
45 * the minimal cut found.
46 *
47 * RETURNS
48 *
49 * The routine returns the capacity of the min cut found. */
50
min_cut(int nn,int ne,const int beg[],const int end[],const cap[],int cut[])51 int min_cut(int nn, int ne, const int beg[/*1+ne*/],
52 const int end[/*1+ne*/], const cap[/*1+ne*/], int cut[/*1+nn*/])
53 { int k;
54 /* sanity checks */
55 xassert(nn >= 2);
56 xassert(ne >= 0);
57 for (k = 1; k <= ne; k++)
58 { xassert(1 <= beg[k] && beg[k] < end[k] && end[k] <= nn);
59 xassert(cap[k] > 0);
60 }
61 /* find min cut */
62 return min_cut_sw(nn, ne, beg, end, cap, cut);
63 }
64
65 /***********************************************************************
66 * NAME
67 *
68 * min_st_cut - find min (s,t)-cut for known max flow
69 *
70 * SYNOPSIS
71 *
72 * #include "mincut.h"
73 *
74 * DESCRIPTION
75 *
76 * This routine finds min (s,t)-cut in a given undirected network that
77 * corresponds to a known max flow from s to t in the network.
78 *
79 * Parameters nn, ne, beg, end, and cap specify the network in the same
80 * way as for the routine min_cut (see above).
81 *
82 * Parameters s and t specify, resp., the number of the source node
83 * and the number of the sink node, s != t, for which the min (s,t)-cut
84 * has to be found.
85 *
86 * Parameter x specifies the known max flow from s to t in the network,
87 * where locations x[1], ..., x[ne] contains elementary flow thru edges
88 * of the network (positive value of x[k] means that the elementary
89 * flow goes from node beg[k] to node end[k], and negative value means
90 * that the flow goes in opposite direction).
91 *
92 * This routine splits the set of nodes V of the network into two
93 * non-empty subsets V(s) and V(t) = V \ V(s), where the source node s
94 * belongs to V(s), the sink node t belongs to V(t), and all edges, one
95 * node of which belongs to V(s) and other one belongs to V(t), are
96 * saturated (that is, x[k] = +cap[k] or x[k] = -cap[k]).
97 *
98 * On exit the routine stores flags of the nodes v[i], i = 1, ..., nn,
99 * to locations cut[i], where cut[i] = 1 means that v[i] belongs to V(s)
100 * and cut[i] = 0 means that v[i] belongs to V(t) = V \ V(s).
101 *
102 * RETURNS
103 *
104 * The routine returns the capacity of min (s,t)-cut, which is the sum
105 * of the capacities of all the edges, which belong to the cut. (Note
106 * that due to theorem by Ford and Fulkerson this value is always equal
107 * to corresponding max flow.)
108 *
109 * ALGORITHM
110 *
111 * To determine the set V(s) the routine simply finds all nodes, which
112 * can be reached from the source node s via non-saturated edges. The
113 * set V(t) is determined as the complement V \ V(s). */
114
min_st_cut(int nn,int ne,const int beg[],const int end[],const int cap[],int s,int t,const int x[],int cut[])115 int min_st_cut(int nn, int ne, const int beg[/*1+ne*/],
116 const int end[/*1+ne*/], const int cap[/*1+ne*/], int s, int t,
117 const int x[/*1+ne*/], int cut[/*1+nn*/])
118 { int i, j, k, p, q, temp, *head1, *next1, *head2, *next2, *list;
119 /* head1[i] points to the first edge with beg[k] = i
120 * next1[k] points to the next edge with the same beg[k]
121 * head2[i] points to the first edge with end[k] = i
122 * next2[k] points to the next edge with the same end[k] */
123 head1 = xalloc(1+nn, sizeof(int));
124 head2 = xalloc(1+nn, sizeof(int));
125 next1 = xalloc(1+ne, sizeof(int));
126 next2 = xalloc(1+ne, sizeof(int));
127 for (i = 1; i <= nn; i++)
128 head1[i] = head2[i] = 0;
129 for (k = 1; k <= ne; k++)
130 { i = beg[k], next1[k] = head1[i], head1[i] = k;
131 j = end[k], next2[k] = head2[j], head2[j] = k;
132 }
133 /* on constructing the set V(s) list[1], ..., list[p-1] contain
134 * nodes, which can be reached from source node and have been
135 * visited, and list[p], ..., list[q] contain nodes, which can be
136 * reached from source node but havn't been visited yet */
137 list = xalloc(1+nn, sizeof(int));
138 for (i = 1; i <= nn; i++)
139 cut[i] = 0;
140 p = q = 1, list[1] = s, cut[s] = 1;
141 while (p <= q)
142 { /* pick next node, which is reachable from the source node and
143 * has not visited yet, and visit it */
144 i = list[p++];
145 /* walk through edges with beg[k] = i */
146 for (k = head1[i]; k != 0; k = next1[k])
147 { j = end[k];
148 xassert(beg[k] == i);
149 /* from v[i] we can reach v[j], if the elementary flow from
150 * v[i] to v[j] is non-saturated */
151 if (cut[j] == 0 && x[k] < +cap[k])
152 list[++q] = j, cut[j] = 1;
153 }
154 /* walk through edges with end[k] = i */
155 for (k = head2[i]; k != 0; k = next2[k])
156 { j = beg[k];
157 xassert(end[k] == i);
158 /* from v[i] we can reach v[j], if the elementary flow from
159 * v[i] to v[j] is non-saturated */
160 if (cut[j] == 0 && x[k] > -cap[k])
161 list[++q] = j, cut[j] = 1;
162 }
163 }
164 /* sink cannot belong to V(s) */
165 xassert(!cut[t]);
166 /* free working arrays */
167 xfree(head1);
168 xfree(head2);
169 xfree(next1);
170 xfree(next2);
171 xfree(list);
172 /* compute capacity of the minimal (s,t)-cut found */
173 temp = 0;
174 for (k = 1; k <= ne; k++)
175 { i = beg[k], j = end[k];
176 if (cut[i] && !cut[j] || !cut[i] && cut[j])
177 temp += cap[k];
178 }
179 /* return to the calling program */
180 return temp;
181 }
182
183 /***********************************************************************
184 * NAME
185 *
186 * min_cut_sw - find min cut with Stoer and Wagner algorithm
187 *
188 * SYNOPSIS
189 *
190 * #include "mincut.h"
191 * int min_cut_sw(int nn, int ne, const int beg[], const int end[],
192 * const cap[], int cut[]);
193 *
194 * DESCRIPTION
195 *
196 * This routine find min cut in a given undirected network with the
197 * algorithm proposed by Stoer and Wagner (see references below).
198 *
199 * Parameters of this routine have the same meaning as for the routine
200 * min_cut (see above).
201 *
202 * RETURNS
203 *
204 * The routine returns the capacity of the min cut found.
205 *
206 * ALGORITHM
207 *
208 * The basic idea of Stoer&Wagner algorithm is the following. Let G be
209 * a capacitated network, and G(s,t) be a network, in which the nodes s
210 * and t are merged into one new node, loops are deleted, but multiple
211 * edges are retained. It is obvious that a minimum cut in G is the
212 * minimum of two quantities: the minimum cut in G(s,t) and a minimum
213 * cut that separates s and t. This allows to find a minimum cut in the
214 * original network solving at most nn max flow problems.
215 *
216 * REFERENCES
217 *
218 * M. Stoer, F. Wagner. A Simple Min Cut Algorithm. Algorithms, ESA'94
219 * LNCS 855 (1994), pp. 141-47.
220 *
221 * J. Cheriyan, R. Ravi. Approximation Algorithms for Network Problems.
222 * Univ. of Waterloo (1998), p. 147. */
223
min_cut_sw(int nn,int ne,const int beg[],const int end[],const cap[],int cut[])224 int min_cut_sw(int nn, int ne, const int beg[/*1+ne*/],
225 const int end[/*1+ne*/], const cap[/*1+ne*/], int cut[/*1+nn*/])
226 { int i, j, k, min_cut, flow, temp, *head1, *next1, *head2, *next2;
227 int I, J, K, S, T, DEG, NV, NE, *HEAD, *NEXT, *NUMB, *BEG, *END,
228 *CAP, *X, *ADJ, *SUM, *CUT;
229 /* head1[i] points to the first edge with beg[k] = i
230 * next1[k] points to the next edge with the same beg[k]
231 * head2[i] points to the first edge with end[k] = i
232 * next2[k] points to the next edge with the same end[k] */
233 head1 = xalloc(1+nn, sizeof(int));
234 head2 = xalloc(1+nn, sizeof(int));
235 next1 = xalloc(1+ne, sizeof(int));
236 next2 = xalloc(1+ne, sizeof(int));
237 for (i = 1; i <= nn; i++)
238 head1[i] = head2[i] = 0;
239 for (k = 1; k <= ne; k++)
240 { i = beg[k], next1[k] = head1[i], head1[i] = k;
241 j = end[k], next2[k] = head2[j], head2[j] = k;
242 }
243 /* an auxiliary network used in the algorithm is resulted from
244 * the original network by merging some nodes into one supernode;
245 * all variables and arrays related to this auxiliary network are
246 * denoted in CAPS */
247 /* HEAD[I] points to the first node of the original network that
248 * belongs to the I-th supernode
249 * NEXT[i] points to the next node of the original network that
250 * belongs to the same supernode as the i-th node
251 * NUMB[i] is a supernode, which the i-th node belongs to */
252 /* initially the auxiliary network is equivalent to the original
253 * network, i.e. each supernode consists of one node */
254 NV = nn;
255 HEAD = xalloc(1+nn, sizeof(int));
256 NEXT = xalloc(1+nn, sizeof(int));
257 NUMB = xalloc(1+nn, sizeof(int));
258 for (i = 1; i <= nn; i++)
259 HEAD[i] = i, NEXT[i] = 0, NUMB[i] = i;
260 /* number of edges in the auxiliary network is never greater than
261 * in the original one */
262 BEG = xalloc(1+ne, sizeof(int));
263 END = xalloc(1+ne, sizeof(int));
264 CAP = xalloc(1+ne, sizeof(int));
265 X = xalloc(1+ne, sizeof(int));
266 /* allocate some auxiliary arrays */
267 ADJ = xalloc(1+nn, sizeof(int));
268 SUM = xalloc(1+nn, sizeof(int));
269 CUT = xalloc(1+nn, sizeof(int));
270 /* currently no min cut is found so far */
271 min_cut = INT_MAX;
272 /* main loop starts here */
273 while (NV > 1)
274 { /* build the set of edges of the auxiliary network */
275 NE = 0;
276 /* multiple edges are not allowed in the max flow algorithm,
277 * so we can replace each multiple edge, which is the result
278 * of merging nodes into supernodes, by a single edge, whose
279 * capacity is the sum of capacities of particular edges;
280 * these summary capacities will be stored in the array SUM */
281 for (I = 1; I <= NV; I++)
282 SUM[I] = 0.0;
283 for (I = 1; I <= NV; I++)
284 { /* DEG is number of single edges, which connects I-th
285 * supernode and some J-th supernode, where I < J */
286 DEG = 0;
287 /* walk thru nodes that belong to I-th supernode */
288 for (i = HEAD[I]; i != 0; i = NEXT[i])
289 { /* i-th node belongs to I-th supernode */
290 /* walk through edges with beg[k] = i */
291 for (k = head1[i]; k != 0; k = next1[k])
292 { j = end[k];
293 /* j-th node belongs to J-th supernode */
294 J = NUMB[j];
295 /* ignore loops and edges with I > J */
296 if (I >= J)
297 continue;
298 /* add an edge that connects I-th and J-th supernodes
299 * (if not added yet) */
300 if (SUM[J] == 0.0)
301 ADJ[++DEG] = J;
302 /* sum up the capacity of the original edge */
303 xassert(cap[k] > 0.0);
304 SUM[J] += cap[k];
305 }
306 /* walk through edges with end[k] = i */
307 for (k = head2[i]; k != 0; k = next2[k])
308 { j = beg[k];
309 /* j-th node belongs to J-th supernode */
310 J = NUMB[j];
311 /* ignore loops and edges with I > J */
312 if (I >= J)
313 continue;
314 /* add an edge that connects I-th and J-th supernodes
315 * (if not added yet) */
316 if (SUM[J] == 0.0)
317 ADJ[++DEG] = J;
318 /* sum up the capacity of the original edge */
319 xassert(cap[k] > 0.0);
320 SUM[J] += cap[k];
321 }
322 }
323 /* add single edges connecting I-th supernode with other
324 * supernodes to the auxiliary network; restore the array
325 * SUM for subsequent use */
326 for (K = 1; K <= DEG; K++)
327 { NE++;
328 xassert(NE <= ne);
329 J = ADJ[K];
330 BEG[NE] = I, END[NE] = J, CAP[NE] = SUM[J];
331 SUM[J] = 0.0;
332 }
333 }
334 /* choose two arbitrary supernodes of the auxiliary network,
335 * one of which is the source and other is the sink */
336 S = 1, T = NV;
337 /* determine max flow from S to T */
338 flow = max_flow(NV, NE, BEG, END, CAP, S, T, X);
339 /* if the min cut that separates supernodes S and T is less
340 * than the currently known, remember it */
341 if (min_cut > flow)
342 { min_cut = flow;
343 /* find min (s,t)-cut in the auxiliary network */
344 temp = min_st_cut(NV, NE, BEG, END, CAP, S, T, X, CUT);
345 /* (Ford and Fulkerson insist on this) */
346 xassert(flow == temp);
347 /* build corresponding min cut in the original network */
348 for (i = 1; i <= nn; i++) cut[i] = CUT[NUMB[i]];
349 /* if the min cut capacity is zero (i.e. the network has
350 * unconnected components), the search can be prematurely
351 * terminated */
352 if (min_cut == 0)
353 break;
354 }
355 /* now merge all nodes of the original network, which belong
356 * to the supernodes S and T, into one new supernode; this is
357 * attained by carrying all nodes from T to S (for the sake of
358 * convenience T should be the last supernode) */
359 xassert(T == NV);
360 /* assign new references to nodes from T */
361 for (i = HEAD[T]; i != 0; i = NEXT[i])
362 NUMB[i] = S;
363 /* find last entry in the node list of S */
364 i = HEAD[S];
365 xassert(i != 0);
366 while (NEXT[i] != 0)
367 i = NEXT[i];
368 /* and attach to it the node list of T */
369 NEXT[i] = HEAD[T];
370 /* decrease number of nodes in the auxiliary network */
371 NV--;
372 }
373 /* free working arrays */
374 xfree(HEAD);
375 xfree(NEXT);
376 xfree(NUMB);
377 xfree(BEG);
378 xfree(END);
379 xfree(CAP);
380 xfree(X);
381 xfree(ADJ);
382 xfree(SUM);
383 xfree(CUT);
384 xfree(head1);
385 xfree(head2);
386 xfree(next1);
387 xfree(next2);
388 /* return to the calling program */
389 return min_cut;
390 }
391
392 /* eof */
393