1 /* -*- mode: C -*- */
2 /*
3 IGraph library.
4 Copyright (C) 2006-2012 Gabor Csardi <csardi.gabor@gmail.com>
5 334 Harvard street, Cambridge, MA 02139 USA
6
7 This program is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
11
12 This program is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
20 02110-1301 USA
21
22 */
23
24 #include "igraph_flow.h"
25 #include "igraph_error.h"
26 #include "igraph_memory.h"
27 #include "igraph_constants.h"
28 #include "igraph_interface.h"
29 #include "igraph_adjlist.h"
30 #include "igraph_conversion.h"
31 #include "igraph_constructors.h"
32 #include "igraph_progress.h"
33 #include "igraph_structural.h"
34 #include "igraph_components.h"
35 #include "igraph_types_internal.h"
36 #include "igraph_math.h"
37 #include "igraph_dqueue.h"
38 #include "igraph_interrupt_internal.h"
39 #include "igraph_topology.h"
40 #include "config.h"
41
42
43 /*
44 * Some general remarks about the functions in this file.
45 *
46 * The following measures can be calculated:
47 * ( 1) s-t maximum flow value, directed graph
48 * ( 2) s-t maximum flow value, undirected graph
49 * ( 3) s-t maximum flow, directed graph
50 * ( 4) s-t maximum flow, undirected graph
51 * ( 5) s-t minimum cut value, directed graph
52 * ( 6) s-t minimum cut value, undirected graph
53 * ( 7) minimum cut value, directed graph
54 * ( 8) minimum cut value, undirected graph
55 * ( 9) s-t minimum cut, directed graph
56 * (10) s-t minimum cut, undirected graph
57 * (11) minimum cut, directed graph
58 * (12) minimum cut, undirected graph
59 * (13) s-t edge connectivity, directed graph
60 * (14) s-t edge connectivity, undirected graph
61 * (15) edge connectivity, directed graph
62 * (16) edge connectivity, undirected graph
63 * (17) s-t vertex connectivity, directed graph
64 * (18) s-t vertex connectivity, undirected graph
65 * (19) vertex connectivity, directed graph
66 * (20) vertex connectivity, undirected graph
67 * (21) s-t number of edge disjoint paths, directed graph
68 * (22) s-t number of edge disjoint paths, undirected graph
69 * (23) s-t number of vertex disjoint paths, directed graph
70 * (24) s-t number of vertex disjoint paths, undirected graph
71 * (25) graph adhesion, directed graph
72 * (26) graph adhesion, undirected graph
73 * (27) graph cohesion, directed graph
74 * (28) graph cohesion, undirected graph
75 *
76 * This is how they are calculated:
77 * ( 1) igraph_maxflow_value, calls igraph_maxflow.
78 * ( 2) igraph_maxflow_value, calls igraph_maxflow, this calls
79 * igraph_i_maxflow_undirected. This transforms the graph into a
80 * directed graph, including two mutual edges instead of every
81 * undirected edge, then igraph_maxflow is called again with the
82 * directed graph.
83 * ( 3) igraph_maxflow, does the push-relabel algorithm, optionally
84 * calculates the cut, the partitions and the flow itself.
85 * ( 4) igraph_maxflow calls igraph_i_maxflow_undirected, this converts
86 * the undirected graph into a directed one, adding two mutual edges
87 * for each undirected edge, then igraph_maxflow is called again,
88 * with the directed graph. After igraph_maxflow returns, we need
89 * to edit the flow (and the cut) to make it sense for the
90 * original graph.
91 * ( 5) igraph_st_mincut_value, we just call igraph_maxflow_value
92 * ( 6) igraph_st_mincut_value, we just call igraph_maxflow_value
93 * ( 7) igraph_mincut_value, we call igraph_maxflow_value (|V|-1)*2
94 * times, from vertex 0 to all other vertices and from all other
95 * vertices to vertex 0
96 * ( 8) We call igraph_i_mincut_value_undirected, that calls
97 * igraph_i_mincut_undirected with partition=partition2=cut=NULL
98 * The Stoer-Wagner algorithm is used.
99 * ( 9) igraph_st_mincut, just calls igraph_maxflow.
100 * (10) igraph_st_mincut, just calls igraph_maxflow.
101 * (11) igraph_mincut, calls igraph_i_mincut_directed, which runs
102 * the maximum flow algorithm 2(|V|-1) times, from vertex zero to
103 * and from all other vertices and stores the smallest cut.
104 * (12) igraph_mincut, igraph_i_mincut_undirected is called,
105 * this is the Stoer-Wagner algorithm
106 * (13) We just call igraph_maxflow_value, back to (1)
107 * (14) We just call igraph_maxflow_value, back to (2)
108 * (15) We just call igraph_mincut_value (possibly after some basic
109 * checks). Back to (7)
110 * (16) We just call igraph_mincut_value (possibly after some basic
111 * checks). Back to (8).
112 * (17) We call igraph_i_st_vertex_connectivity_directed.
113 * That creates a new graph with 2*|V| vertices and smartly chosen
114 * edges, so that the s-t edge connectivity of this graph is the
115 * same as the s-t vertex connectivity of the original graph.
116 * So finally it calls igraph_maxflow_value, go to (1)
117 * (18) We call igraph_i_st_vertex_connectivity_undirected.
118 * We convert the graph to a directed one,
119 * IGRAPH_TO_DIRECTED_MUTUAL method. Then we call
120 * igraph_i_st_vertex_connectivity_directed, see (17).
121 * (19) We call igraph_i_vertex_connectivity_directed.
122 * That calls igraph_st_vertex_connectivity for all pairs of
123 * vertices. Back to (17).
124 * (20) We call igraph_i_vertex_connectivity_undirected.
125 * That converts the graph into a directed one
126 * (IGRAPH_TO_DIRECTED_MUTUAL) and calls the directed version,
127 * igraph_i_vertex_connectivity_directed, see (19).
128 * (21) igraph_edge_disjoint_paths, we just call igraph_maxflow_value, (1).
129 * (22) igraph_edge_disjoint_paths, we just call igraph_maxflow_value, (2).
130 * (23) igraph_vertex_disjoint_paths, if there is a connection between
131 * the two vertices, then we remove that (or all of them if there
132 * are many), as this could mess up vertex connectivity
133 * calculation. The we call
134 * igraph_i_st_vertex_connectivity_directed, see (19).
135 * (24) igraph_vertex_disjoint_paths, if there is a connection between
136 * the two vertices, then we remove that (or all of them if there
137 * are many), as this could mess up vertex connectivity
138 * calculation. The we call
139 * igraph_i_st_vertex_connectivity_undirected, see (20).
140 * (25) We just call igraph_edge_connectivity, see (15).
141 * (26) We just call igraph_edge_connectivity, see (16).
142 * (27) We just call igraph_vertex_connectivity, see (19).
143 * (28) We just call igraph_vertex_connectivity, see (20).
144 */
145
146 /*
147 * This is an internal function that calculates the maximum flow value
148 * on undirected graphs, either for an s-t vertex pair or for the
149 * graph (i.e. all vertex pairs).
150 *
151 * It does it by converting the undirected graph to a corresponding
152 * directed graph, including reciprocal directed edges instead of each
153 * undirected edge.
154 */
155
igraph_i_maxflow_undirected(const igraph_t * graph,igraph_real_t * value,igraph_vector_t * flow,igraph_vector_t * cut,igraph_vector_t * partition,igraph_vector_t * partition2,igraph_integer_t source,igraph_integer_t target,const igraph_vector_t * capacity,igraph_maxflow_stats_t * stats)156 static int igraph_i_maxflow_undirected(const igraph_t *graph,
157 igraph_real_t *value,
158 igraph_vector_t *flow,
159 igraph_vector_t *cut,
160 igraph_vector_t *partition,
161 igraph_vector_t *partition2,
162 igraph_integer_t source,
163 igraph_integer_t target,
164 const igraph_vector_t *capacity,
165 igraph_maxflow_stats_t *stats) {
166 igraph_integer_t no_of_edges = (igraph_integer_t) igraph_ecount(graph);
167 igraph_integer_t no_of_nodes = (igraph_integer_t) igraph_vcount(graph);
168 igraph_vector_t edges;
169 igraph_vector_t newcapacity;
170 igraph_t newgraph;
171 long int i;
172
173 /* We need to convert this to directed by hand, since we need to be
174 sure that the edge ids will be handled properly to build the new
175 capacity vector. */
176
177 IGRAPH_VECTOR_INIT_FINALLY(&edges, 0);
178 IGRAPH_VECTOR_INIT_FINALLY(&newcapacity, no_of_edges * 2);
179 IGRAPH_CHECK(igraph_vector_reserve(&edges, no_of_edges * 4));
180 IGRAPH_CHECK(igraph_get_edgelist(graph, &edges, 0));
181 IGRAPH_CHECK(igraph_vector_resize(&edges, no_of_edges * 4));
182 for (i = 0; i < no_of_edges; i++) {
183 VECTOR(edges)[no_of_edges * 2 + i * 2] = VECTOR(edges)[i * 2 + 1];
184 VECTOR(edges)[no_of_edges * 2 + i * 2 + 1] = VECTOR(edges)[i * 2];
185 VECTOR(newcapacity)[i] = VECTOR(newcapacity)[no_of_edges + i] =
186 capacity ? VECTOR(*capacity)[i] : 1.0;
187 }
188
189 IGRAPH_CHECK(igraph_create(&newgraph, &edges, no_of_nodes, IGRAPH_DIRECTED));
190 IGRAPH_FINALLY(igraph_destroy, &newgraph);
191
192 IGRAPH_CHECK(igraph_maxflow(&newgraph, value, flow, cut, partition,
193 partition2, source, target, &newcapacity, stats));
194
195 if (cut) {
196 long int i, cs = igraph_vector_size(cut);
197 for (i = 0; i < cs; i++) {
198 if (VECTOR(*cut)[i] >= no_of_edges) {
199 VECTOR(*cut)[i] -= no_of_edges;
200 }
201 }
202 }
203
204 /* The flow has one non-zero value for each real-nonreal edge pair,
205 by definition, we convert it to a positive-negative vector. If
206 for an edge the flow is negative that means that it is going
207 from the bigger vertex id to the smaller one. For positive
208 values the direction is the opposite. */
209 if (flow) {
210 long int i;
211 for (i = 0; i < no_of_edges; i++) {
212 VECTOR(*flow)[i] -= VECTOR(*flow)[i + no_of_edges];
213 }
214 IGRAPH_CHECK(igraph_vector_resize(flow, no_of_edges));
215 }
216
217 igraph_destroy(&newgraph);
218 igraph_vector_destroy(&edges);
219 igraph_vector_destroy(&newcapacity);
220 IGRAPH_FINALLY_CLEAN(3);
221
222 return 0;
223 }
224
225 #define FIRST(i) (VECTOR(*first)[(i)])
226 #define LAST(i) (VECTOR(*first)[(i)+1])
227 #define CURRENT(i) (VECTOR(*current)[(i)])
228 #define RESCAP(i) (VECTOR(*rescap)[(i)])
229 #define REV(i) (VECTOR(*rev)[(i)])
230 #define HEAD(i) (VECTOR(*to)[(i)])
231 #define EXCESS(i) (VECTOR(*excess)[(i)])
232 #define DIST(i) (VECTOR(*distance)[(i)])
233 #define DISCHARGE(v) (igraph_i_mf_discharge((v), ¤t, &first, &rescap, \
234 &to, &distance, &excess, \
235 no_of_nodes, source, target, \
236 &buckets, &ibuckets, \
237 &rev, stats, &npushsince, \
238 &nrelabelsince))
239 #define PUSH(v,e,n) (igraph_i_mf_push((v), (e), (n), current, rescap, \
240 excess, target, source, buckets, \
241 ibuckets, distance, rev, stats, \
242 npushsince))
243 #define RELABEL(v) (igraph_i_mf_relabel((v), no_of_nodes, distance, \
244 first, rescap, to, current, \
245 stats, nrelabelsince))
246 #define GAP(b) (igraph_i_mf_gap((b), stats, buckets, ibuckets, \
247 no_of_nodes, distance))
248 #define BFS() (igraph_i_mf_bfs(&bfsq, source, target, no_of_nodes, \
249 &buckets, &ibuckets, &distance, \
250 &first, ¤t, &to, &excess, \
251 &rescap, &rev))
252
igraph_i_mf_gap(long int b,igraph_maxflow_stats_t * stats,igraph_buckets_t * buckets,igraph_dbuckets_t * ibuckets,long int no_of_nodes,igraph_vector_long_t * distance)253 static void igraph_i_mf_gap(long int b, igraph_maxflow_stats_t *stats,
254 igraph_buckets_t *buckets, igraph_dbuckets_t *ibuckets,
255 long int no_of_nodes,
256 igraph_vector_long_t *distance) {
257
258 long int bo;
259 (stats->nogap)++;
260 for (bo = b + 1; bo <= no_of_nodes; bo++) {
261 while (!igraph_dbuckets_empty_bucket(ibuckets, bo)) {
262 long int n = igraph_dbuckets_pop(ibuckets, bo);
263 (stats->nogapnodes)++;
264 DIST(n) = no_of_nodes;
265 }
266 }
267 }
268
igraph_i_mf_relabel(long int v,long int no_of_nodes,igraph_vector_long_t * distance,igraph_vector_long_t * first,igraph_vector_t * rescap,igraph_vector_long_t * to,igraph_vector_long_t * current,igraph_maxflow_stats_t * stats,int * nrelabelsince)269 static void igraph_i_mf_relabel(long int v, long int no_of_nodes,
270 igraph_vector_long_t *distance,
271 igraph_vector_long_t *first,
272 igraph_vector_t *rescap, igraph_vector_long_t *to,
273 igraph_vector_long_t *current,
274 igraph_maxflow_stats_t *stats, int *nrelabelsince) {
275
276 long int min = no_of_nodes;
277 long int k, l, min_edge = 0;
278 (stats->norelabel)++; (*nrelabelsince)++;
279 DIST(v) = no_of_nodes;
280 for (k = FIRST(v), l = LAST(v); k < l; k++) {
281 if (RESCAP(k) > 0 && DIST(HEAD(k)) < min) {
282 min = DIST(HEAD(k));
283 min_edge = k;
284 }
285 }
286 min++;
287 if (min < no_of_nodes) {
288 DIST(v) = min;
289 CURRENT(v) = min_edge;
290 }
291 }
292
igraph_i_mf_push(long int v,long int e,long int n,igraph_vector_long_t * current,igraph_vector_t * rescap,igraph_vector_t * excess,long int target,long int source,igraph_buckets_t * buckets,igraph_dbuckets_t * ibuckets,igraph_vector_long_t * distance,igraph_vector_long_t * rev,igraph_maxflow_stats_t * stats,int * npushsince)293 static void igraph_i_mf_push(long int v, long int e, long int n,
294 igraph_vector_long_t *current,
295 igraph_vector_t *rescap, igraph_vector_t *excess,
296 long int target, long int source,
297 igraph_buckets_t *buckets, igraph_dbuckets_t *ibuckets,
298 igraph_vector_long_t *distance,
299 igraph_vector_long_t *rev, igraph_maxflow_stats_t *stats,
300 int *npushsince) {
301 igraph_real_t delta =
302 RESCAP(e) < EXCESS(v) ? RESCAP(e) : EXCESS(v);
303 (stats->nopush)++; (*npushsince)++;
304 if (EXCESS(n) == 0 && n != target) {
305 igraph_dbuckets_delete(ibuckets, DIST(n), n);
306 igraph_buckets_add(buckets, (long int) DIST(n), n);
307 }
308 RESCAP(e) -= delta;
309 RESCAP(REV(e)) += delta;
310 EXCESS(n) += delta;
311 EXCESS(v) -= delta;
312 }
313
igraph_i_mf_discharge(long int v,igraph_vector_long_t * current,igraph_vector_long_t * first,igraph_vector_t * rescap,igraph_vector_long_t * to,igraph_vector_long_t * distance,igraph_vector_t * excess,long int no_of_nodes,long int source,long int target,igraph_buckets_t * buckets,igraph_dbuckets_t * ibuckets,igraph_vector_long_t * rev,igraph_maxflow_stats_t * stats,int * npushsince,int * nrelabelsince)314 static void igraph_i_mf_discharge(long int v,
315 igraph_vector_long_t *current,
316 igraph_vector_long_t *first,
317 igraph_vector_t *rescap,
318 igraph_vector_long_t *to,
319 igraph_vector_long_t *distance,
320 igraph_vector_t *excess,
321 long int no_of_nodes, long int source,
322 long int target, igraph_buckets_t *buckets,
323 igraph_dbuckets_t *ibuckets,
324 igraph_vector_long_t *rev,
325 igraph_maxflow_stats_t *stats,
326 int *npushsince, int *nrelabelsince) {
327 do {
328 long int i;
329 long int start = (long int) CURRENT(v);
330 long int stop = (long int) LAST(v);
331 for (i = start; i < stop; i++) {
332 if (RESCAP(i) > 0) {
333 long int nei = HEAD(i);
334 if (DIST(v) == DIST(nei) + 1) {
335 PUSH((v), i, nei);
336 if (EXCESS(v) == 0) {
337 break;
338 }
339 }
340 }
341 }
342 if (i == stop) {
343 long int origdist = DIST(v);
344 RELABEL(v);
345 if (igraph_buckets_empty_bucket(buckets, origdist) &&
346 igraph_dbuckets_empty_bucket(ibuckets, origdist)) {
347 GAP(origdist);
348 }
349 if (DIST(v) == no_of_nodes) {
350 break;
351 }
352 } else {
353 CURRENT(v) = i;
354 igraph_dbuckets_add(ibuckets, DIST(v), v);
355 break;
356 }
357 } while (1);
358 }
359
igraph_i_mf_bfs(igraph_dqueue_long_t * bfsq,long int source,long int target,long int no_of_nodes,igraph_buckets_t * buckets,igraph_dbuckets_t * ibuckets,igraph_vector_long_t * distance,igraph_vector_long_t * first,igraph_vector_long_t * current,igraph_vector_long_t * to,igraph_vector_t * excess,igraph_vector_t * rescap,igraph_vector_long_t * rev)360 static void igraph_i_mf_bfs(igraph_dqueue_long_t *bfsq,
361 long int source, long int target,
362 long int no_of_nodes, igraph_buckets_t *buckets,
363 igraph_dbuckets_t *ibuckets,
364 igraph_vector_long_t *distance,
365 igraph_vector_long_t *first, igraph_vector_long_t *current,
366 igraph_vector_long_t *to, igraph_vector_t *excess,
367 igraph_vector_t *rescap, igraph_vector_long_t *rev) {
368
369 long int k, l;
370
371 igraph_buckets_clear(buckets);
372 igraph_dbuckets_clear(ibuckets);
373 igraph_vector_long_fill(distance, no_of_nodes);
374 DIST(target) = 0;
375
376 igraph_dqueue_long_push(bfsq, target);
377 while (!igraph_dqueue_long_empty(bfsq)) {
378 long int node = igraph_dqueue_long_pop(bfsq);
379 long int ndist = DIST(node) + 1;
380 for (k = FIRST(node), l = LAST(node); k < l; k++) {
381 if (RESCAP(REV(k)) > 0) {
382 long int nei = HEAD(k);
383 if (DIST(nei) == no_of_nodes) {
384 DIST(nei) = ndist;
385 CURRENT(nei) = FIRST(nei);
386 if (EXCESS(nei) > 0) {
387 igraph_buckets_add(buckets, ndist, nei);
388 } else {
389 igraph_dbuckets_add(ibuckets, ndist, nei);
390 }
391 igraph_dqueue_long_push(bfsq, nei);
392 }
393 }
394 }
395 }
396 }
397
398 /**
399 * \function igraph_maxflow
400 * Maximum network flow between a pair of vertices
401 *
402 * </para><para>This function implements the Goldberg-Tarjan algorithm for
403 * calculating value of the maximum flow in a directed or undirected
404 * graph. The algorithm was given in Andrew V. Goldberg, Robert
405 * E. Tarjan: A New Approach to the Maximum-Flow Problem, Journal of
406 * the ACM, 35(4), 921-940, 1988. </para>
407 *
408 * <para> The input of the function is a graph, a vector
409 * of real numbers giving the capacity of the edges and two vertices
410 * of the graph, the source and the target. A flow is a function
411 * assigning positive real numbers to the edges and satisfying two
412 * requirements: (1) the flow value is less than the capacity of the
413 * edge and (2) at each vertex except the source and the target, the
414 * incoming flow (i.e. the sum of the flow on the incoming edges) is
415 * the same as the outgoing flow (i.e. the sum of the flow on the
416 * outgoing edges). The value of the flow is the incoming flow at the
417 * target vertex. The maximum flow is the flow with the maximum
418 * value.
419 *
420 * \param graph The input graph, either directed or undirected.
421 * \param value Pointer to a real number, the value of the maximum
422 * will be placed here, unless it is a null pointer.
423 * \param flow If not a null pointer, then it must be a pointer to an
424 * initialized vector. The vector will be resized, and the flow
425 * on each edge will be placed in it, in the order of the edge
426 * ids. For undirected graphs this argument is bit trickier,
427 * since for these the flow direction is not predetermined by
428 * the edge direction. For these graphs the elements of the
429 * \p flow vector can be negative, this means that the flow
430 * goes from the bigger vertex id to the smaller one. Positive
431 * values mean that the flow goes from the smaller vertex id to
432 * the bigger one.
433 * \param cut A null pointer or a pointer to an initialized vector.
434 * If not a null pointer, then the minimum cut corresponding to
435 * the maximum flow is stored here, i.e. all edge ids that are
436 * part of the minimum cut are stored in the vector.
437 * \param partition A null pointer or a pointer to an initialized
438 * vector. If not a null pointer, then the first partition of
439 * the minimum cut that corresponds to the maximum flow will be
440 * placed here. The first partition is always the one that
441 * contains the source vertex.
442 * \param partition2 A null pointer or a pointer to an initialized
443 * vector. If not a null pointer, then the second partition of
444 * the minimum cut that corresponds to the maximum flow will be
445 * placed here. The second partition is always the one that
446 * contains the target vertex.
447 * \param source The id of the source vertex.
448 * \param target The id of the target vertex.
449 * \param capacity Vector containing the capacity of the edges. If NULL, then
450 * every edge is considered to have capacity 1.0.
451 * \param stats Counts of the number of different operations
452 * preformed by the algorithm are stored here.
453 * \return Error code.
454 *
455 * Time complexity: O(|V|^3). In practice it is much faster, but i
456 * cannot prove a better lower bound for the data structure i've
457 * used. In fact, this implementation runs much faster than the
458 * \c hi_pr implementation discussed in
459 * B. V. Cherkassky and A. V. Goldberg: On implementing the
460 * push-relabel method for the maximum flow problem, (Algorithmica,
461 * 19:390--410, 1997) on all the graph classes i've tried.
462 *
463 * \sa \ref igraph_mincut_value(), \ref igraph_edge_connectivity(),
464 * \ref igraph_vertex_connectivity() for
465 * properties based on the maximum flow.
466 *
467 * \example examples/simple/flow.c
468 * \example examples/simple/flow2.c
469 */
470
igraph_maxflow(const igraph_t * graph,igraph_real_t * value,igraph_vector_t * flow,igraph_vector_t * cut,igraph_vector_t * partition,igraph_vector_t * partition2,igraph_integer_t source,igraph_integer_t target,const igraph_vector_t * capacity,igraph_maxflow_stats_t * stats)471 int igraph_maxflow(const igraph_t *graph, igraph_real_t *value,
472 igraph_vector_t *flow, igraph_vector_t *cut,
473 igraph_vector_t *partition, igraph_vector_t *partition2,
474 igraph_integer_t source, igraph_integer_t target,
475 const igraph_vector_t *capacity,
476 igraph_maxflow_stats_t *stats) {
477
478 igraph_integer_t no_of_nodes = (igraph_integer_t) igraph_vcount(graph);
479 igraph_integer_t no_of_orig_edges = (igraph_integer_t) igraph_ecount(graph);
480 igraph_integer_t no_of_edges = 2 * no_of_orig_edges;
481
482 igraph_vector_t rescap, excess;
483 igraph_vector_long_t from, to, rev, distance;
484 igraph_vector_t edges, rank;
485 igraph_vector_long_t current, first;
486 igraph_buckets_t buckets;
487 igraph_dbuckets_t ibuckets;
488
489 igraph_dqueue_long_t bfsq;
490
491 long int i, j, idx;
492 int npushsince = 0, nrelabelsince = 0;
493
494 igraph_maxflow_stats_t local_stats; /* used if the user passed a null pointer for stats */
495
496 if (stats == 0) {
497 stats = &local_stats;
498 }
499
500 if (!igraph_is_directed(graph)) {
501 IGRAPH_CHECK(igraph_i_maxflow_undirected(graph, value, flow, cut,
502 partition, partition2, source,
503 target, capacity, stats));
504 return 0;
505 }
506
507 if (capacity && igraph_vector_size(capacity) != no_of_orig_edges) {
508 IGRAPH_ERROR("Invalid capacity vector", IGRAPH_EINVAL);
509 }
510 if (source < 0 || source >= no_of_nodes || target < 0 || target >= no_of_nodes) {
511 IGRAPH_ERROR("Invalid source or target vertex", IGRAPH_EINVAL);
512 }
513
514 stats->nopush = stats->norelabel = stats->nogap = stats->nogapnodes =
515 stats->nobfs = 0;
516
517 /*
518 * The data structure:
519 * - First of all, we consider every edge twice, first the edge
520 * itself, but also its opposite.
521 * - (from, to) contain all edges (original + opposite), ordered by
522 * the id of the source vertex. During the algorithm we just need
523 * 'to', so from is destroyed soon. We only need it in the
524 * beginning, to create the 'first' pointers.
525 * - 'first' is a pointer vector for 'to', first[i] points to the
526 * first neighbor of vertex i and first[i+1]-1 is the last
527 * neighbor of vertex i. (Unless vertex i is isolate, in which
528 * case first[i]==first[i+1]).
529 * - 'rev' contains a mapping from an edge to its opposite pair
530 * - 'rescap' contains the residual capacities of the edges, this is
531 * initially equal to the capacity of the edges for the original
532 * edges and it is zero for the opposite edges.
533 * - 'excess' contains the excess flow for the vertices. I.e. the flow
534 * that is coming in, but it is not going out.
535 * - 'current' stores the next neighboring vertex to check, for every
536 * vertex, when excess flow is being pushed to neighbors.
537 * - 'distance' stores the distance of the vertices from the source.
538 * - 'rank' and 'edges' are only needed temporarily, for ordering and
539 * storing the edges.
540 * - we use an igraph_buckets_t data structure ('buckets') to find
541 * the vertices with the highest 'distance' values quickly.
542 * This always contains the vertices that have a positive excess
543 * flow.
544 */
545 #undef FIRST
546 #undef LAST
547 #undef CURRENT
548 #undef RESCAP
549 #undef REV
550 #undef HEAD
551 #undef EXCESS
552 #undef DIST
553 #define FIRST(i) (VECTOR(first)[(i)])
554 #define LAST(i) (VECTOR(first)[(i)+1])
555 #define CURRENT(i) (VECTOR(current)[(i)])
556 #define RESCAP(i) (VECTOR(rescap)[(i)])
557 #define REV(i) (VECTOR(rev)[(i)])
558 #define HEAD(i) (VECTOR(to)[(i)])
559 #define EXCESS(i) (VECTOR(excess)[(i)])
560 #define DIST(i) (VECTOR(distance)[(i)])
561
562 igraph_dqueue_long_init(&bfsq, no_of_nodes);
563 IGRAPH_FINALLY(igraph_dqueue_long_destroy, &bfsq);
564 IGRAPH_VECTOR_LONG_INIT_FINALLY(&to, no_of_edges);
565 IGRAPH_VECTOR_LONG_INIT_FINALLY(&rev, no_of_edges);
566 IGRAPH_VECTOR_INIT_FINALLY(&rescap, no_of_edges);
567 IGRAPH_VECTOR_INIT_FINALLY(&excess, no_of_nodes);
568 IGRAPH_VECTOR_LONG_INIT_FINALLY(&distance, no_of_nodes);
569 IGRAPH_VECTOR_LONG_INIT_FINALLY(&first, no_of_nodes + 1);
570
571 IGRAPH_VECTOR_INIT_FINALLY(&rank, no_of_edges);
572 IGRAPH_VECTOR_LONG_INIT_FINALLY(&from, no_of_edges);
573 IGRAPH_VECTOR_INIT_FINALLY(&edges, no_of_edges);
574
575 /* Create the basic data structure */
576 IGRAPH_CHECK(igraph_get_edgelist(graph, &edges, 0));
577 IGRAPH_CHECK(igraph_vector_rank(&edges, &rank, no_of_nodes));
578
579 for (i = 0; i < no_of_edges; i += 2) {
580 long int pos = (long int) VECTOR(rank)[i];
581 long int pos2 = (long int) VECTOR(rank)[i + 1];
582 VECTOR(from)[pos] = VECTOR(edges)[i];
583 VECTOR(to)[pos] = VECTOR(edges)[i + 1];
584 VECTOR(from)[pos2] = VECTOR(edges)[i + 1];
585 VECTOR(to)[pos2] = VECTOR(edges)[i];
586 VECTOR(rev)[pos] = pos2;
587 VECTOR(rev)[pos2] = pos;
588 VECTOR(rescap)[pos] = capacity ? VECTOR(*capacity)[i / 2] : 1.0;
589 VECTOR(rescap)[pos2] = 0.0;
590 }
591
592 /* The first pointers. This is a but trickier, than one would
593 think, because of the possible isolate vertices. */
594
595 idx = -1;
596 for (i = 0; i <= VECTOR(from)[0]; i++) {
597 idx++; VECTOR(first)[idx] = 0;
598 }
599 for (i = 1; i < no_of_edges; i++) {
600 long int n = (long int) (VECTOR(from)[i] -
601 VECTOR(from)[ (long int) VECTOR(first)[idx] ]);
602 for (j = 0; j < n; j++) {
603 idx++; VECTOR(first)[idx] = i;
604 }
605 }
606 idx++;
607 while (idx < no_of_nodes + 1) {
608 VECTOR(first)[idx++] = no_of_edges;
609 }
610
611 igraph_vector_long_destroy(&from);
612 igraph_vector_destroy(&edges);
613 IGRAPH_FINALLY_CLEAN(2);
614
615 if (!flow) {
616 igraph_vector_destroy(&rank);
617 IGRAPH_FINALLY_CLEAN(1);
618 }
619
620 /* And the current pointers, initially the same as the first */
621 IGRAPH_VECTOR_LONG_INIT_FINALLY(¤t, no_of_nodes);
622 for (i = 0; i < no_of_nodes; i++) {
623 VECTOR(current)[i] = VECTOR(first)[i];
624 }
625
626 /* OK, the graph is set up, initialization */
627
628 IGRAPH_CHECK(igraph_buckets_init(&buckets, no_of_nodes + 1, no_of_nodes));
629 IGRAPH_FINALLY(igraph_buckets_destroy, &buckets);
630 IGRAPH_CHECK(igraph_dbuckets_init(&ibuckets, no_of_nodes + 1, no_of_nodes));
631 IGRAPH_FINALLY(igraph_dbuckets_destroy, &ibuckets);
632
633 /* Send as much flow as possible from the source to its neighbors */
634 for (i = FIRST(source), j = LAST(source); i < j; i++) {
635 if (HEAD(i) != source) {
636 igraph_real_t delta = RESCAP(i);
637 RESCAP(i) = 0;
638 RESCAP(REV(i)) += delta;
639 EXCESS(HEAD(i)) += delta;
640 }
641 }
642
643 BFS();
644 (stats->nobfs)++;
645
646 while (!igraph_buckets_empty(&buckets)) {
647 long int vertex = igraph_buckets_popmax(&buckets);
648 DISCHARGE(vertex);
649 if (npushsince > no_of_nodes / 2 && nrelabelsince > no_of_nodes) {
650 (stats->nobfs)++;
651 BFS();
652 npushsince = nrelabelsince = 0;
653 }
654 }
655
656 /* Store the result */
657 if (value) {
658 *value = EXCESS(target);
659 }
660
661 /* If we also need the minimum cut */
662 if (cut || partition || partition2) {
663 /* We need to find all vertices from which the target is reachable
664 in the residual graph. We do a breadth-first search, going
665 backwards. */
666 igraph_dqueue_t Q;
667 igraph_vector_bool_t added;
668 long int marked = 0;
669
670 IGRAPH_CHECK(igraph_vector_bool_init(&added, no_of_nodes));
671 IGRAPH_FINALLY(igraph_vector_bool_destroy, &added);
672
673 IGRAPH_CHECK(igraph_dqueue_init(&Q, 100));
674 IGRAPH_FINALLY(igraph_dqueue_destroy, &Q);
675
676 igraph_dqueue_push(&Q, target);
677 VECTOR(added)[(long int)target] = 1;
678 marked++;
679 while (!igraph_dqueue_empty(&Q)) {
680 long int actnode = (long int) igraph_dqueue_pop(&Q);
681 for (i = FIRST(actnode), j = LAST(actnode); i < j; i++) {
682 long int nei = HEAD(i);
683 if (!VECTOR(added)[nei] && RESCAP(REV(i)) > 0.0) {
684 VECTOR(added)[nei] = 1;
685 marked++;
686 IGRAPH_CHECK(igraph_dqueue_push(&Q, nei));
687 }
688 }
689 }
690 igraph_dqueue_destroy(&Q);
691 IGRAPH_FINALLY_CLEAN(1);
692
693 /* Now we marked each vertex that is on one side of the cut,
694 check the crossing edges */
695
696 if (cut) {
697 igraph_vector_clear(cut);
698 for (i = 0; i < no_of_orig_edges; i++) {
699 long int f = IGRAPH_FROM(graph, i);
700 long int t = IGRAPH_TO(graph, i);
701 if (!VECTOR(added)[f] && VECTOR(added)[t]) {
702 IGRAPH_CHECK(igraph_vector_push_back(cut, i));
703 }
704 }
705 }
706
707 if (partition2) {
708 long int x = 0;
709 IGRAPH_CHECK(igraph_vector_resize(partition2, marked));
710 for (i = 0; i < no_of_nodes; i++) {
711 if (VECTOR(added)[i]) {
712 VECTOR(*partition2)[x++] = i;
713 }
714 }
715 }
716
717 if (partition) {
718 long int x = 0;
719 IGRAPH_CHECK(igraph_vector_resize(partition,
720 no_of_nodes - marked));
721 for (i = 0; i < no_of_nodes; i++) {
722 if (!VECTOR(added)[i]) {
723 VECTOR(*partition)[x++] = i;
724 }
725 }
726 }
727
728 igraph_vector_bool_destroy(&added);
729 IGRAPH_FINALLY_CLEAN(1);
730 }
731
732 if (flow) {
733 /* Initialize the backward distances, with a breadth-first search
734 from the source */
735 igraph_dqueue_t Q;
736 igraph_vector_int_t added;
737 long int j, k, l;
738 igraph_t flow_graph;
739 igraph_vector_t flow_edges;
740 igraph_bool_t dag;
741
742 IGRAPH_CHECK(igraph_vector_int_init(&added, no_of_nodes));
743 IGRAPH_FINALLY(igraph_vector_int_destroy, &added);
744 IGRAPH_CHECK(igraph_dqueue_init(&Q, 100));
745 IGRAPH_FINALLY(igraph_dqueue_destroy, &Q);
746
747 igraph_dqueue_push(&Q, source);
748 igraph_dqueue_push(&Q, 0);
749 VECTOR(added)[(long int)source] = 1;
750 while (!igraph_dqueue_empty(&Q)) {
751 long int actnode = (long int) igraph_dqueue_pop(&Q);
752 long int actdist = (long int) igraph_dqueue_pop(&Q);
753 DIST(actnode) = actdist;
754
755 for (i = FIRST(actnode), j = LAST(actnode); i < j; i++) {
756 long int nei = HEAD(i);
757 if (!VECTOR(added)[nei] && RESCAP(REV(i)) > 0.0) {
758 VECTOR(added)[nei] = 1;
759 IGRAPH_CHECK(igraph_dqueue_push(&Q, nei));
760 IGRAPH_CHECK(igraph_dqueue_push(&Q, actdist + 1));
761 }
762 }
763 } /* !igraph_dqueue_empty(&Q) */
764
765 igraph_vector_int_destroy(&added);
766 igraph_dqueue_destroy(&Q);
767 IGRAPH_FINALLY_CLEAN(2);
768
769 /* Reinitialize the buckets */
770 igraph_buckets_clear(&buckets);
771 for (i = 0; i < no_of_nodes; i++) {
772 if (EXCESS(i) > 0.0 && i != source && i != target) {
773 igraph_buckets_add(&buckets, (long int) DIST(i), i);
774 }
775 }
776
777 /* Now we return the flow to the source */
778 while (!igraph_buckets_empty(&buckets)) {
779 long int vertex = igraph_buckets_popmax(&buckets);
780
781 /* DISCHARGE(vertex) comes here */
782 do {
783 for (i = (long int) CURRENT(vertex), j = LAST(vertex); i < j; i++) {
784 if (RESCAP(i) > 0) {
785 long int nei = HEAD(i);
786
787 if (DIST(vertex) == DIST(nei) + 1) {
788 igraph_real_t delta =
789 RESCAP(i) < EXCESS(vertex) ? RESCAP(i) : EXCESS(vertex);
790 RESCAP(i) -= delta;
791 RESCAP(REV(i)) += delta;
792
793 if (nei != source && EXCESS(nei) == 0.0 &&
794 DIST(nei) != no_of_nodes) {
795 igraph_buckets_add(&buckets, (long int) DIST(nei), nei);
796 }
797
798 EXCESS(nei) += delta;
799 EXCESS(vertex) -= delta;
800
801 if (EXCESS(vertex) == 0) {
802 break;
803 }
804
805 }
806 }
807 }
808
809 if (i == j) {
810
811 /* RELABEL(vertex) comes here */
812 igraph_real_t min;
813 long int min_edge = 0;
814 DIST(vertex) = min = no_of_nodes;
815 for (k = FIRST(vertex), l = LAST(vertex); k < l; k++) {
816 if (RESCAP(k) > 0) {
817 if (DIST(HEAD(k)) < min) {
818 min = DIST(HEAD(k));
819 min_edge = k;
820 }
821 }
822 }
823
824 min++;
825
826 if (min < no_of_nodes) {
827 DIST(vertex) = min;
828 CURRENT(vertex) = min_edge;
829 /* Vertex is still active */
830 igraph_buckets_add(&buckets, (long int) DIST(vertex), vertex);
831 }
832
833 /* TODO: gap heuristics here ??? */
834
835 } else {
836 CURRENT(vertex) = FIRST(vertex);
837 }
838
839 break;
840
841 } while (1);
842 }
843
844 /* We need to eliminate flow cycles now. Before that we check that
845 there is a cycle in the flow graph.
846
847 First we do a couple of DFSes from the source vertex to the
848 target and factor out the paths we find. If there is no more
849 path to the target, then all remaining flow must be in flow
850 cycles, so we don't need it at all.
851
852 Some details. 'stack' contains the whole path of the DFS, both
853 the vertices and the edges, they are alternating in the stack.
854 'current' helps finding the next outgoing edge of a vertex
855 quickly, the next edge of 'v' is FIRST(v)+CURRENT(v). If this
856 is LAST(v), then there are no more edges to try.
857
858 The 'added' vector contains 0 if the vertex was not visited
859 before, 1 if it is currently in 'stack', and 2 if it is not in
860 'stack', but it was visited before. */
861
862 IGRAPH_VECTOR_INIT_FINALLY(&flow_edges, 0);
863 for (i = 0, j = 0; i < no_of_edges; i += 2, j++) {
864 long int pos = (long int) VECTOR(rank)[i];
865 if ((capacity ? VECTOR(*capacity)[j] : 1.0) > RESCAP(pos)) {
866 IGRAPH_CHECK(igraph_vector_push_back(&flow_edges,
867 IGRAPH_FROM(graph, j)));
868 IGRAPH_CHECK(igraph_vector_push_back(&flow_edges,
869 IGRAPH_TO(graph, j)));
870 }
871 }
872 IGRAPH_CHECK(igraph_create(&flow_graph, &flow_edges, no_of_nodes,
873 IGRAPH_DIRECTED));
874 igraph_vector_destroy(&flow_edges);
875 IGRAPH_FINALLY_CLEAN(1);
876 IGRAPH_FINALLY(igraph_destroy, &flow_graph);
877 IGRAPH_CHECK(igraph_is_dag(&flow_graph, &dag));
878 igraph_destroy(&flow_graph);
879 IGRAPH_FINALLY_CLEAN(1);
880
881 if (!dag) {
882 igraph_vector_long_t stack;
883 igraph_vector_t mycap;
884
885 IGRAPH_CHECK(igraph_vector_long_init(&stack, 0));
886 IGRAPH_FINALLY(igraph_vector_long_destroy, &stack);
887 IGRAPH_CHECK(igraph_vector_int_init(&added, no_of_nodes));
888 IGRAPH_FINALLY(igraph_vector_int_destroy, &added);
889 IGRAPH_VECTOR_INIT_FINALLY(&mycap, no_of_edges);
890
891 #define MYCAP(i) (VECTOR(mycap)[(i)])
892
893 for (i = 0; i < no_of_edges; i += 2) {
894 long int pos = (long int) VECTOR(rank)[i];
895 long int pos2 = (long int) VECTOR(rank)[i + 1];
896 MYCAP(pos) = (capacity ? VECTOR(*capacity)[i / 2] : 1.0) - RESCAP(pos);
897 MYCAP(pos2) = 0.0;
898 }
899
900 do {
901 igraph_vector_long_null(¤t);
902 igraph_vector_long_clear(&stack);
903 igraph_vector_int_null(&added);
904
905 IGRAPH_CHECK(igraph_vector_long_push_back(&stack, -1));
906 IGRAPH_CHECK(igraph_vector_long_push_back(&stack, source));
907 VECTOR(added)[(long int)source] = 1;
908 while (!igraph_vector_long_empty(&stack) &&
909 igraph_vector_long_tail(&stack) != target) {
910 long int actnode = igraph_vector_long_tail(&stack);
911 long int edge = FIRST(actnode) + (long int) CURRENT(actnode);
912 long int nei;
913 while (edge < LAST(actnode) && MYCAP(edge) == 0.0) {
914 edge++;
915 }
916 nei = edge < LAST(actnode) ? HEAD(edge) : -1;
917
918 if (edge < LAST(actnode) && !VECTOR(added)[nei]) {
919 /* Go forward along next edge, if the vertex was not
920 visited before */
921 IGRAPH_CHECK(igraph_vector_long_push_back(&stack, edge));
922 IGRAPH_CHECK(igraph_vector_long_push_back(&stack, nei));
923 VECTOR(added)[nei] = 1;
924 CURRENT(actnode) += 1;
925 } else if (edge < LAST(actnode) && VECTOR(added)[nei] == 1) {
926 /* We found a flow cycle, factor it out. Go back in stack
927 until we find 'nei' again, determine the flow along the
928 cycle. */
929 igraph_real_t thisflow = MYCAP(edge);
930 long int idx;
931 for (idx = igraph_vector_long_size(&stack) - 2;
932 idx >= 0 && VECTOR(stack)[idx + 1] != nei; idx -= 2) {
933 long int e = VECTOR(stack)[idx];
934 igraph_real_t rcap = e >= 0 ? MYCAP(e) : MYCAP(edge);
935 if (rcap < thisflow) {
936 thisflow = rcap;
937 }
938 }
939 MYCAP(edge) -= thisflow; RESCAP(edge) += thisflow;
940 for (idx = igraph_vector_long_size(&stack) - 2;
941 idx >= 0 && VECTOR(stack)[idx + 1] != nei; idx -= 2) {
942 long int e = VECTOR(stack)[idx];
943 if (e >= 0) {
944 MYCAP(e) -= thisflow;
945 RESCAP(e) += thisflow;
946 }
947 }
948 CURRENT(actnode) += 1;
949 } else if (edge < LAST(actnode)) { /* && VECTOR(added)[nei]==2 */
950 /* The next edge leads to a vertex that was visited before,
951 but it is currently not in 'stack' */
952 CURRENT(actnode) += 1;
953 } else {
954 /* Go backward, take out the node and the edge that leads to it */
955 igraph_vector_long_pop_back(&stack);
956 igraph_vector_long_pop_back(&stack);
957 VECTOR(added)[actnode] = 2;
958 }
959 }
960
961 /* If non-empty, then it contains a path from source to target
962 in the residual graph. We factor out this path from the flow. */
963 if (!igraph_vector_long_empty(&stack)) {
964 long int pl = igraph_vector_long_size(&stack);
965 igraph_real_t thisflow = EXCESS(target);
966 for (i = 2; i < pl; i += 2) {
967 long int edge = VECTOR(stack)[i];
968 igraph_real_t rcap = MYCAP(edge);
969 if (rcap < thisflow) {
970 thisflow = rcap;
971 }
972 }
973 for (i = 2; i < pl; i += 2) {
974 long int edge = VECTOR(stack)[i];
975 MYCAP(edge) -= thisflow;
976 }
977 }
978
979 } while (!igraph_vector_long_empty(&stack));
980
981 igraph_vector_destroy(&mycap);
982 igraph_vector_int_destroy(&added);
983 igraph_vector_long_destroy(&stack);
984 IGRAPH_FINALLY_CLEAN(3);
985 }
986
987 /* ----------------------------------------------------------- */
988
989 IGRAPH_CHECK(igraph_vector_resize(flow, no_of_orig_edges));
990 for (i = 0, j = 0; i < no_of_edges; i += 2, j++) {
991 long int pos = (long int) VECTOR(rank)[i];
992 VECTOR(*flow)[j] = (capacity ? VECTOR(*capacity)[j] : 1.0) -
993 RESCAP(pos);
994 }
995
996 igraph_vector_destroy(&rank);
997 IGRAPH_FINALLY_CLEAN(1);
998 }
999
1000 igraph_dbuckets_destroy(&ibuckets);
1001 igraph_buckets_destroy(&buckets);
1002 igraph_vector_long_destroy(¤t);
1003 igraph_vector_long_destroy(&first);
1004 igraph_vector_long_destroy(&distance);
1005 igraph_vector_destroy(&excess);
1006 igraph_vector_destroy(&rescap);
1007 igraph_vector_long_destroy(&rev);
1008 igraph_vector_long_destroy(&to);
1009 igraph_dqueue_long_destroy(&bfsq);
1010 IGRAPH_FINALLY_CLEAN(10);
1011
1012 return 0;
1013 }
1014
1015 /**
1016 * \function igraph_maxflow_value
1017 * \brief Maximum flow in a network with the push/relabel algorithm
1018 *
1019 * </para><para>This function implements the Goldberg-Tarjan algorithm for
1020 * calculating value of the maximum flow in a directed or undirected
1021 * graph. The algorithm was given in Andrew V. Goldberg, Robert
1022 * E. Tarjan: A New Approach to the Maximum-Flow Problem, Journal of
1023 * the ACM, 35(4), 921-940, 1988. </para>
1024 *
1025 * <para> The input of the function is a graph, a vector
1026 * of real numbers giving the capacity of the edges and two vertices
1027 * of the graph, the source and the target. A flow is a function
1028 * assigning positive real numbers to the edges and satisfying two
1029 * requirements: (1) the flow value is less than the capacity of the
1030 * edge and (2) at each vertex except the source and the target, the
1031 * incoming flow (i.e. the sum of the flow on the incoming edges) is
1032 * the same as the outgoing flow (i.e. the sum of the flow on the
1033 * outgoing edges). The value of the flow is the incoming flow at the
1034 * target vertex. The maximum flow is the flow with the maximum
1035 * value. </para>
1036 *
1037 * <para> According to a theorem by Ford and Fulkerson
1038 * (L. R. Ford Jr. and D. R. Fulkerson. Maximal flow through a
1039 * network. Canadian J. Math., 8:399-404, 1956.) the maximum flow
1040 * between two vertices is the same as the
1041 * minimum cut between them (also called the minimum s-t cut). So \ref
1042 * igraph_st_mincut_value() gives the same result in all cases as \c
1043 * igraph_maxflow_value().</para>
1044 *
1045 * <para> Note that the value of the maximum flow is the same as the
1046 * minimum cut in the graph.
1047 * \param graph The input graph, either directed or undirected.
1048 * \param value Pointer to a real number, the result will be placed here.
1049 * \param source The id of the source vertex.
1050 * \param target The id of the target vertex.
1051 * \param capacity Vector containing the capacity of the edges. If NULL, then
1052 * every edge is considered to have capacity 1.0.
1053 * \param stats Counts of the number of different operations
1054 * preformed by the algorithm are stored here.
1055 * \return Error code.
1056 *
1057 * Time complexity: O(|V|^3).
1058 *
1059 * \sa \ref igraph_maxflow() to calculate the actual flow.
1060 * \ref igraph_mincut_value(), \ref igraph_edge_connectivity(),
1061 * \ref igraph_vertex_connectivity() for
1062 * properties based on the maximum flow.
1063 */
1064
igraph_maxflow_value(const igraph_t * graph,igraph_real_t * value,igraph_integer_t source,igraph_integer_t target,const igraph_vector_t * capacity,igraph_maxflow_stats_t * stats)1065 int igraph_maxflow_value(const igraph_t *graph, igraph_real_t *value,
1066 igraph_integer_t source, igraph_integer_t target,
1067 const igraph_vector_t *capacity,
1068 igraph_maxflow_stats_t *stats) {
1069
1070 return igraph_maxflow(graph, value, /*flow=*/ 0, /*cut=*/ 0,
1071 /*partition=*/ 0, /*partition1=*/ 0,
1072 source, target, capacity, stats);
1073 }
1074
1075 /**
1076 * \function igraph_st_mincut_value
1077 * \brief The minimum s-t cut in a graph
1078 *
1079 * </para><para> The minimum s-t cut in a weighted (=valued) graph is the
1080 * total minimum edge weight needed to remove from the graph to
1081 * eliminate all paths from a given vertex (\c source) to
1082 * another vertex (\c target). Directed paths are considered in
1083 * directed graphs, and undirected paths in undirected graphs. </para>
1084 *
1085 * <para> The minimum s-t cut between two vertices is known to be same
1086 * as the maximum flow between these two vertices. So this function
1087 * calls \ref igraph_maxflow_value() to do the calculation.
1088 * \param graph The input graph.
1089 * \param value Pointer to a real variable, the result will be stored
1090 * here.
1091 * \param source The id of the source vertex.
1092 * \param target The id of the target vertex.
1093 * \param capacity Pointer to the capacity vector, it should contain
1094 * non-negative numbers and its length should be the same the
1095 * the number of edges in the graph. It can be a null pointer, then
1096 * every edge has unit capacity.
1097 * \return Error code.
1098 *
1099 * Time complexity: O(|V|^3), see also the discussion for \ref
1100 * igraph_maxflow_value(), |V| is the number of vertices.
1101 */
1102
igraph_st_mincut_value(const igraph_t * graph,igraph_real_t * value,igraph_integer_t source,igraph_integer_t target,const igraph_vector_t * capacity)1103 int igraph_st_mincut_value(const igraph_t *graph, igraph_real_t *value,
1104 igraph_integer_t source, igraph_integer_t target,
1105 const igraph_vector_t *capacity) {
1106
1107 if (source == target) {
1108 IGRAPH_ERROR("source and target vertices are the same", IGRAPH_EINVAL);
1109 }
1110
1111 IGRAPH_CHECK(igraph_maxflow_value(graph, value, source, target, capacity, 0));
1112
1113 return 0;
1114 }
1115
1116 /**
1117 * \function igraph_st_mincut
1118 * Minimum cut between a source and a target vertex
1119 *
1120 * Finds the edge set that has the smallest total capacity among all
1121 * edge sets that disconnect the source and target vertices.
1122 *
1123 * </para><para>The calculation is performed using maximum flow
1124 * techniques, by calling \ref igraph_maxflow().
1125 * \param graph The input graph.
1126 * \param value Pointer to a real variable, the value of the cut is
1127 * stored here.
1128 * \param cut Pointer to a real vector, the edge ids that are included
1129 * in the cut are stored here. This argument is ignored if it
1130 * is a null pointer.
1131 * \param partition Pointer to a real vector, the vertex ids of the
1132 * vertices in the first partition of the cut are stored
1133 * here. The first partition is always the one that contains the
1134 * source vertex. This argument is ignored if it is a null pointer.
1135 * \param partition2 Pointer to a real vector, the vertex ids of the
1136 * vertices in the second partition of the cut are stored here.
1137 * The second partition is always the one that contains the
1138 * target vertex. This argument is ignored if it is a null pointer.
1139 * \param source Integer, the id of the source vertex.
1140 * \param target Integer, the id of the target vertex.
1141 * \param capacity Vector containing the capacity of the edges. If a
1142 * null pointer, then every edge is considered to have capacity
1143 * 1.0.
1144 * \return Error code.
1145 *
1146 * \sa \ref igraph_maxflow().
1147 *
1148 * Time complexity: see \ref igraph_maxflow().
1149 */
1150
igraph_st_mincut(const igraph_t * graph,igraph_real_t * value,igraph_vector_t * cut,igraph_vector_t * partition,igraph_vector_t * partition2,igraph_integer_t source,igraph_integer_t target,const igraph_vector_t * capacity)1151 int igraph_st_mincut(const igraph_t *graph, igraph_real_t *value,
1152 igraph_vector_t *cut, igraph_vector_t *partition,
1153 igraph_vector_t *partition2,
1154 igraph_integer_t source, igraph_integer_t target,
1155 const igraph_vector_t *capacity) {
1156
1157 return igraph_maxflow(graph, value, /*flow=*/ 0,
1158 cut, partition, partition2,
1159 source, target, capacity, 0);
1160 }
1161
1162 /* This is a flow-based version, but there is a better one
1163 for undirected graphs */
1164
1165 /* int igraph_i_mincut_value_undirected(const igraph_t *graph, */
1166 /* igraph_real_t *res, */
1167 /* const igraph_vector_t *capacity) { */
1168
1169 /* long int no_of_edges=igraph_ecount(graph); */
1170 /* long int no_of_nodes=igraph_vcount(graph); */
1171 /* igraph_vector_t edges; */
1172 /* igraph_vector_t newcapacity; */
1173 /* igraph_t newgraph; */
1174 /* long int i; */
1175
1176 /* /\* We need to convert this to directed by hand, since we need to be */
1177 /* sure that the edge ids will be handled properly to build the new */
1178 /* capacity vector. *\/ */
1179
1180 /* IGRAPH_VECTOR_INIT_FINALLY(&edges, 0); */
1181 /* IGRAPH_VECTOR_INIT_FINALLY(&newcapacity, no_of_edges*2); */
1182 /* IGRAPH_CHECK(igraph_vector_reserve(&edges, no_of_edges*4)); */
1183 /* IGRAPH_CHECK(igraph_get_edgelist(graph, &edges, 0)); */
1184 /* IGRAPH_CHECK(igraph_vector_resize(&edges, no_of_edges*4)); */
1185 /* for (i=0; i<no_of_edges; i++) { */
1186 /* VECTOR(edges)[no_of_edges*2+i*2] = VECTOR(edges)[i*2+1]; */
1187 /* VECTOR(edges)[no_of_edges*2+i*2+1] = VECTOR(edges)[i*2]; */
1188 /* VECTOR(newcapacity)[i] = VECTOR(newcapacity)[no_of_edges+i] = */
1189 /* capacity ? VECTOR(*capacity)[i] : 1.0 ; */
1190 /* } */
1191
1192 /* IGRAPH_CHECK(igraph_create(&newgraph, &edges, no_of_nodes, IGRAPH_DIRECTED)); */
1193 /* IGRAPH_FINALLY(igraph_destroy, &newgraph); */
1194
1195 /* IGRAPH_CHECK(igraph_mincut_value(&newgraph, res, &newcapacity)); */
1196
1197 /* igraph_destroy(&newgraph); */
1198 /* igraph_vector_destroy(&edges); */
1199 /* igraph_vector_destroy(&newcapacity); */
1200 /* IGRAPH_FINALLY_CLEAN(3); */
1201
1202 /* return 0; */
1203 /* } */
1204
1205 /*
1206 * This is the Stoer-Wagner algorithm, it works for calculating the
1207 * minimum cut for undirected graphs, for the whole graph.
1208 * I.e. this is basically the edge-connectivity of the graph.
1209 * It can also calculate the cut itself, not just the cut value.
1210 */
1211
igraph_i_mincut_undirected(const igraph_t * graph,igraph_real_t * res,igraph_vector_t * partition,igraph_vector_t * partition2,igraph_vector_t * cut,const igraph_vector_t * capacity)1212 static int igraph_i_mincut_undirected(const igraph_t *graph,
1213 igraph_real_t *res,
1214 igraph_vector_t *partition,
1215 igraph_vector_t *partition2,
1216 igraph_vector_t *cut,
1217 const igraph_vector_t *capacity) {
1218
1219 igraph_integer_t no_of_nodes = (igraph_integer_t) igraph_vcount(graph);
1220 igraph_integer_t no_of_edges = (igraph_integer_t) igraph_ecount(graph);
1221
1222 igraph_i_cutheap_t heap;
1223 igraph_real_t mincut = IGRAPH_INFINITY; /* infinity */
1224 long int i;
1225
1226 igraph_adjlist_t adjlist;
1227 igraph_inclist_t inclist;
1228
1229 igraph_vector_t mergehist;
1230 igraph_bool_t calc_cut = partition || partition2 || cut;
1231 long int act_step = 0, mincut_step = 0;
1232
1233 if (capacity && igraph_vector_size(capacity) != no_of_edges) {
1234 IGRAPH_ERROR("Invalid capacity vector size", IGRAPH_EINVAL);
1235 }
1236
1237 /* Check if the graph is connected at all */
1238 {
1239 igraph_vector_t memb, csize;
1240 igraph_integer_t no;
1241 IGRAPH_VECTOR_INIT_FINALLY(&memb, 0);
1242 IGRAPH_VECTOR_INIT_FINALLY(&csize, 0);
1243 IGRAPH_CHECK(igraph_clusters(graph, &memb, &csize, &no,
1244 /*mode=*/ IGRAPH_WEAK));
1245 if (no != 1) {
1246 if (res) {
1247 *res = 0;
1248 }
1249 if (cut) {
1250 igraph_vector_clear(cut);
1251 }
1252 if (partition) {
1253 int j = 0;
1254 IGRAPH_CHECK(igraph_vector_resize(partition,
1255 (long int) VECTOR(csize)[0]));
1256 for (i = 0; i < no_of_nodes; i++) {
1257 if (VECTOR(memb)[i] == 0) {
1258 VECTOR(*partition)[j++] = i;
1259 }
1260 }
1261 }
1262 if (partition2) {
1263 int j = 0;
1264 IGRAPH_CHECK(igraph_vector_resize(partition2, no_of_nodes -
1265 (long int) VECTOR(csize)[0]));
1266 for (i = 0; i < no_of_nodes; i++) {
1267 if (VECTOR(memb)[i] != 0) {
1268 VECTOR(*partition2)[j++] = i;
1269 }
1270 }
1271 }
1272 }
1273 igraph_vector_destroy(&csize);
1274 igraph_vector_destroy(&memb);
1275 IGRAPH_FINALLY_CLEAN(2);
1276
1277 if (no != 1) {
1278 return 0;
1279 }
1280 }
1281
1282 if (calc_cut) {
1283 IGRAPH_VECTOR_INIT_FINALLY(&mergehist, 0);
1284 IGRAPH_CHECK(igraph_vector_reserve(&mergehist, no_of_nodes * 2));
1285 }
1286
1287 IGRAPH_CHECK(igraph_i_cutheap_init(&heap, no_of_nodes));
1288 IGRAPH_FINALLY(igraph_i_cutheap_destroy, &heap);
1289
1290 IGRAPH_CHECK(igraph_inclist_init(graph, &inclist, IGRAPH_OUT));
1291 IGRAPH_FINALLY(igraph_inclist_destroy, &inclist);
1292
1293 IGRAPH_CHECK(igraph_adjlist_init(graph, &adjlist, IGRAPH_OUT));
1294 IGRAPH_FINALLY(igraph_adjlist_destroy, &adjlist);
1295
1296 while (igraph_i_cutheap_size(&heap) >= 2) {
1297
1298 long int last;
1299 igraph_real_t acut;
1300 long int a, n;
1301
1302 igraph_vector_int_t *edges, *edges2;
1303 igraph_vector_int_t *neis, *neis2;
1304
1305 do {
1306 a = igraph_i_cutheap_popmax(&heap);
1307
1308 /* update the weights of the active vertices connected to a */
1309 edges = igraph_inclist_get(&inclist, a);
1310 neis = igraph_adjlist_get(&adjlist, a);
1311 n = igraph_vector_int_size(edges);
1312 for (i = 0; i < n; i++) {
1313 igraph_integer_t edge = (igraph_integer_t) VECTOR(*edges)[i];
1314 igraph_integer_t to = (igraph_integer_t) VECTOR(*neis)[i];
1315 igraph_real_t weight = capacity ? VECTOR(*capacity)[(long int)edge] : 1.0;
1316 igraph_i_cutheap_update(&heap, to, weight);
1317 }
1318
1319 } while (igraph_i_cutheap_active_size(&heap) > 1);
1320
1321 /* Now, there is only one active vertex left,
1322 calculate the cut of the phase */
1323 acut = igraph_i_cutheap_maxvalue(&heap);
1324 last = igraph_i_cutheap_popmax(&heap);
1325
1326 if (acut < mincut) {
1327 mincut = acut;
1328 mincut_step = act_step;
1329 }
1330
1331 if (mincut == 0) {
1332 break;
1333 }
1334
1335 /* And contract the last and the remaining vertex (a and last) */
1336 /* Before actually doing that, make some notes */
1337 act_step++;
1338 if (calc_cut) {
1339 IGRAPH_CHECK(igraph_vector_push_back(&mergehist, a));
1340 IGRAPH_CHECK(igraph_vector_push_back(&mergehist, last));
1341 }
1342 /* First remove the a--last edge if there is one, a is still the
1343 last deactivated vertex */
1344 edges = igraph_inclist_get(&inclist, a);
1345 neis = igraph_adjlist_get(&adjlist, a);
1346 n = igraph_vector_int_size(edges);
1347 for (i = 0; i < n; ) {
1348 if (VECTOR(*neis)[i] == last) {
1349 VECTOR(*neis)[i] = VECTOR(*neis)[n - 1];
1350 VECTOR(*edges)[i] = VECTOR(*edges)[n - 1];
1351 igraph_vector_int_pop_back(neis);
1352 igraph_vector_int_pop_back(edges);
1353 n--;
1354 } else {
1355 i++;
1356 }
1357 }
1358
1359 edges = igraph_inclist_get(&inclist, last);
1360 neis = igraph_adjlist_get(&adjlist, last);
1361 n = igraph_vector_int_size(edges);
1362 for (i = 0; i < n; ) {
1363 if (VECTOR(*neis)[i] == a) {
1364 VECTOR(*neis)[i] = VECTOR(*neis)[n - 1];
1365 VECTOR(*edges)[i] = VECTOR(*edges)[n - 1];
1366 igraph_vector_int_pop_back(neis);
1367 igraph_vector_int_pop_back(edges);
1368 n--;
1369 } else {
1370 i++;
1371 }
1372 }
1373
1374 /* Now rewrite the edge lists of last's neighbors */
1375 neis = igraph_adjlist_get(&adjlist, last);
1376 n = igraph_vector_int_size(neis);
1377 for (i = 0; i < n; i++) {
1378 igraph_integer_t nei = (igraph_integer_t) VECTOR(*neis)[i];
1379 long int n2, j;
1380 neis2 = igraph_adjlist_get(&adjlist, nei);
1381 n2 = igraph_vector_int_size(neis2);
1382 for (j = 0; j < n2; j++) {
1383 if (VECTOR(*neis2)[j] == last) {
1384 VECTOR(*neis2)[j] = a;
1385 }
1386 }
1387 }
1388
1389 /* And append the lists of last to the lists of a */
1390 edges = igraph_inclist_get(&inclist, a);
1391 neis = igraph_adjlist_get(&adjlist, a);
1392 edges2 = igraph_inclist_get(&inclist, last);
1393 neis2 = igraph_adjlist_get(&adjlist, last);
1394 IGRAPH_CHECK(igraph_vector_int_append(edges, edges2));
1395 IGRAPH_CHECK(igraph_vector_int_append(neis, neis2));
1396 igraph_vector_int_clear(edges2); /* TODO: free it */
1397 igraph_vector_int_clear(neis2); /* TODO: free it */
1398
1399 /* Remove the deleted vertex from the heap entirely */
1400 igraph_i_cutheap_reset_undefine(&heap, last);
1401 }
1402
1403 *res = mincut;
1404
1405 igraph_inclist_destroy(&inclist);
1406 igraph_adjlist_destroy(&adjlist);
1407 igraph_i_cutheap_destroy(&heap);
1408 IGRAPH_FINALLY_CLEAN(3);
1409
1410 if (calc_cut) {
1411 long int bignode = (long int) VECTOR(mergehist)[2 * mincut_step + 1];
1412 long int i, idx;
1413 long int size = 1;
1414 char *mark;
1415 mark = igraph_Calloc(no_of_nodes, char);
1416 if (!mark) {
1417 IGRAPH_ERROR("Not enough memory for minimum cut", IGRAPH_ENOMEM);
1418 }
1419 IGRAPH_FINALLY(igraph_free, mark);
1420
1421 /* first count the vertices in the partition */
1422 mark[bignode] = 1;
1423 for (i = mincut_step - 1; i >= 0; i--) {
1424 if ( mark[ (long int) VECTOR(mergehist)[2 * i] ] ) {
1425 size++;
1426 mark [ (long int) VECTOR(mergehist)[2 * i + 1] ] = 1;
1427 }
1428 }
1429
1430 /* now store them, if requested */
1431 if (partition) {
1432 IGRAPH_CHECK(igraph_vector_resize(partition, size));
1433 idx = 0;
1434 VECTOR(*partition)[idx++] = bignode;
1435 for (i = mincut_step - 1; i >= 0; i--) {
1436 if (mark[ (long int) VECTOR(mergehist)[2 * i] ]) {
1437 VECTOR(*partition)[idx++] = VECTOR(mergehist)[2 * i + 1];
1438 }
1439 }
1440 }
1441
1442 /* The other partition too? */
1443 if (partition2) {
1444 IGRAPH_CHECK(igraph_vector_resize(partition2, no_of_nodes - size));
1445 idx = 0;
1446 for (i = 0; i < no_of_nodes; i++) {
1447 if (!mark[i]) {
1448 VECTOR(*partition2)[idx++] = i;
1449 }
1450 }
1451 }
1452
1453 /* The edges in the cut are also requested? */
1454 /* We want as few memory allocated for 'cut' as possible,
1455 so we first collect the edges in mergehist, we don't
1456 need that anymore. Then we copy it to 'cut'; */
1457 if (cut) {
1458 igraph_integer_t from, to;
1459 igraph_vector_clear(&mergehist);
1460 for (i = 0; i < no_of_edges; i++) {
1461 igraph_edge(graph, (igraph_integer_t) i, &from, &to);
1462 if ((mark[(long int)from] && !mark[(long int)to]) ||
1463 (mark[(long int)to] && !mark[(long int)from])) {
1464 IGRAPH_CHECK(igraph_vector_push_back(&mergehist, i));
1465 }
1466 }
1467 igraph_vector_clear(cut);
1468 IGRAPH_CHECK(igraph_vector_append(cut, &mergehist));
1469 }
1470
1471 igraph_free(mark);
1472 igraph_vector_destroy(&mergehist);
1473 IGRAPH_FINALLY_CLEAN(2);
1474 }
1475
1476 return 0;
1477 }
1478
igraph_i_mincut_directed(const igraph_t * graph,igraph_real_t * value,igraph_vector_t * partition,igraph_vector_t * partition2,igraph_vector_t * cut,const igraph_vector_t * capacity)1479 static int igraph_i_mincut_directed(const igraph_t *graph,
1480 igraph_real_t *value,
1481 igraph_vector_t *partition,
1482 igraph_vector_t *partition2,
1483 igraph_vector_t *cut,
1484 const igraph_vector_t *capacity) {
1485 long int i;
1486 long int no_of_nodes = igraph_vcount(graph);
1487 igraph_real_t flow;
1488 igraph_real_t minmaxflow = IGRAPH_INFINITY;
1489 igraph_vector_t mypartition, mypartition2, mycut;
1490 igraph_vector_t *ppartition = 0, *ppartition2 = 0, *pcut = 0;
1491 igraph_vector_t bestpartition, bestpartition2, bestcut;
1492
1493 if (partition) {
1494 IGRAPH_VECTOR_INIT_FINALLY(&bestpartition, 0);
1495 }
1496 if (partition2) {
1497 IGRAPH_VECTOR_INIT_FINALLY(&bestpartition2, 0);
1498 }
1499 if (cut) {
1500 IGRAPH_VECTOR_INIT_FINALLY(&bestcut, 0);
1501 }
1502
1503 if (partition) {
1504 IGRAPH_VECTOR_INIT_FINALLY(&mypartition, 0);
1505 ppartition = &mypartition;
1506 }
1507 if (partition2) {
1508 IGRAPH_VECTOR_INIT_FINALLY(&mypartition2, 0);
1509 ppartition2 = &mypartition2;
1510 }
1511 if (cut) {
1512 IGRAPH_VECTOR_INIT_FINALLY(&mycut, 0);
1513 pcut = &mycut;
1514 }
1515
1516 for (i = 1; i < no_of_nodes; i++) {
1517 IGRAPH_CHECK(igraph_maxflow(graph, /*value=*/ &flow, /*flow=*/ 0,
1518 pcut, ppartition, ppartition2, /*source=*/ 0,
1519 /*target=*/ (igraph_integer_t) i, capacity, 0));
1520 if (flow < minmaxflow) {
1521 minmaxflow = flow;
1522 if (cut) {
1523 IGRAPH_CHECK(igraph_vector_update(&bestcut, &mycut));
1524 }
1525 if (partition) {
1526 IGRAPH_CHECK(igraph_vector_update(&bestpartition, &mypartition));
1527 }
1528 if (partition2) {
1529 IGRAPH_CHECK(igraph_vector_update(&bestpartition2, &mypartition2));
1530 }
1531
1532 if (minmaxflow == 0) {
1533 break;
1534 }
1535 }
1536 IGRAPH_CHECK(igraph_maxflow(graph, /*value=*/ &flow, /*flow=*/ 0,
1537 pcut, ppartition, ppartition2,
1538 /*source=*/ (igraph_integer_t) i,
1539 /*target=*/ 0, capacity, 0));
1540 if (flow < minmaxflow) {
1541 minmaxflow = flow;
1542 if (cut) {
1543 IGRAPH_CHECK(igraph_vector_update(&bestcut, &mycut));
1544 }
1545 if (partition) {
1546 IGRAPH_CHECK(igraph_vector_update(&bestpartition, &mypartition));
1547 }
1548 if (partition2) {
1549 IGRAPH_CHECK(igraph_vector_update(&bestpartition2, &mypartition2));
1550 }
1551
1552 if (minmaxflow == 0) {
1553 break;
1554 }
1555 }
1556 }
1557
1558 if (value) {
1559 *value = minmaxflow;
1560 }
1561
1562 if (cut) {
1563 igraph_vector_destroy(&mycut);
1564 IGRAPH_FINALLY_CLEAN(1);
1565 }
1566 if (partition) {
1567 igraph_vector_destroy(&mypartition);
1568 IGRAPH_FINALLY_CLEAN(1);
1569 }
1570 if (partition2) {
1571 igraph_vector_destroy(&mypartition2);
1572 IGRAPH_FINALLY_CLEAN(1);
1573 }
1574 if (cut) {
1575 IGRAPH_CHECK(igraph_vector_update(cut, &bestcut));
1576 igraph_vector_destroy(&bestcut);
1577 IGRAPH_FINALLY_CLEAN(1);
1578 }
1579 if (partition2) {
1580 IGRAPH_CHECK(igraph_vector_update(partition2, &bestpartition2));
1581 igraph_vector_destroy(&bestpartition2);
1582 IGRAPH_FINALLY_CLEAN(1);
1583 }
1584 if (partition) {
1585 IGRAPH_CHECK(igraph_vector_update(partition, &bestpartition));
1586 igraph_vector_destroy(&bestpartition);
1587 IGRAPH_FINALLY_CLEAN(1);
1588 }
1589
1590 return 0;
1591 }
1592
1593 /**
1594 * \function igraph_mincut
1595 * \brief Calculates the minimum cut in a graph.
1596 *
1597 * This function calculates the minimum cut in a graph.
1598 * The minimum cut is the minimum set of edges which needs to be
1599 * removed to disconnect the graph. The minimum is calculated using
1600 * the weights (\p capacity) of the edges, so the cut with the minimum
1601 * total capacity is calculated.
1602 *
1603 * </para><para> For directed graphs an implementation based on
1604 * calculating 2|V|-2 maximum flows is used.
1605 * For undirected graphs we use the Stoer-Wagner
1606 * algorithm, as described in M. Stoer and F. Wagner: A simple min-cut
1607 * algorithm, Journal of the ACM, 44 585-591, 1997.
1608 *
1609 * </para><para>
1610 * The first implementation of the actual cut calculation for
1611 * undirected graphs was made by Gregory Benison, thanks Greg.
1612 * \param graph The input graph.
1613 * \param value Pointer to a float, the value of the cut will be
1614 * stored here.
1615 * \param partition Pointer to an initialized vector, the ids
1616 * of the vertices in the first partition after separating the
1617 * graph will be stored here. The vector will be resized as
1618 * needed. This argument is ignored if it is a NULL pointer.
1619 * \param partition2 Pointer to an initialized vector the ids
1620 * of the vertices in the second partition will be stored here.
1621 * The vector will be resized as needed. This argument is ignored
1622 * if it is a NULL pointer.
1623 * \param cut Pointer to an initialized vector, the ids of the edges
1624 * in the cut will be stored here. This argument is ignored if it
1625 * is a NULL pointer.
1626 * \param capacity A numeric vector giving the capacities of the
1627 * edges. If a null pointer then all edges have unit capacity.
1628 * \return Error code.
1629 *
1630 * \sa \ref igraph_mincut_value(), a simpler interface for calculating
1631 * the value of the cut only.
1632 *
1633 * Time complexity: for directed graphs it is O(|V|^4), but see the
1634 * remarks at \ref igraph_maxflow(). For undirected graphs it is
1635 * O(|V||E|+|V|^2 log|V|). |V| and |E| are the number of vertices and
1636 * edges respectively.
1637 *
1638 * \example examples/simple/igraph_mincut.c
1639 */
1640
igraph_mincut(const igraph_t * graph,igraph_real_t * value,igraph_vector_t * partition,igraph_vector_t * partition2,igraph_vector_t * cut,const igraph_vector_t * capacity)1641 int igraph_mincut(const igraph_t *graph,
1642 igraph_real_t *value,
1643 igraph_vector_t *partition,
1644 igraph_vector_t *partition2,
1645 igraph_vector_t *cut,
1646 const igraph_vector_t *capacity) {
1647
1648 if (igraph_is_directed(graph)) {
1649 if (partition || partition2 || cut) {
1650 igraph_i_mincut_directed(graph, value, partition, partition2, cut,
1651 capacity);
1652 } else {
1653 return igraph_mincut_value(graph, value, capacity);
1654 }
1655 } else {
1656 IGRAPH_CHECK(igraph_i_mincut_undirected(graph, value, partition,
1657 partition2, cut, capacity));
1658 return IGRAPH_SUCCESS;
1659 }
1660
1661 return 0;
1662 }
1663
1664
igraph_i_mincut_value_undirected(const igraph_t * graph,igraph_real_t * res,const igraph_vector_t * capacity)1665 static int igraph_i_mincut_value_undirected(const igraph_t *graph,
1666 igraph_real_t *res,
1667 const igraph_vector_t *capacity) {
1668 return igraph_i_mincut_undirected(graph, res, 0, 0, 0, capacity);
1669 }
1670
1671 /**
1672 * \function igraph_mincut_value
1673 * \brief The minimum edge cut in a graph
1674 *
1675 * </para><para> The minimum edge cut in a graph is the total minimum
1676 * weight of the edges needed to remove from the graph to make the
1677 * graph \em not strongly connected. (If the original graph is not
1678 * strongly connected then this is zero.) Note that in undirected
1679 * graphs strong connectedness is the same as weak connectedness. </para>
1680 *
1681 * <para> The minimum cut can be calculated with maximum flow
1682 * techniques, although the current implementation does this only for
1683 * directed graphs and a separate non-flow based implementation is
1684 * used for undirected graphs. See Mechthild Stoer and Frank Wagner: A
1685 * simple min-cut algorithm, Journal of the ACM 44 585--591, 1997.
1686 * For directed graphs
1687 * the maximum flow is calculated between a fixed vertex and all the
1688 * other vertices in the graph and this is done in both
1689 * directions. Then the minimum is taken to get the minimum cut.
1690 *
1691 * \param graph The input graph.
1692 * \param res Pointer to a real variable, the result will be stored
1693 * here.
1694 * \param capacity Pointer to the capacity vector, it should contain
1695 * the same number of non-negative numbers as the number of edges in
1696 * the graph. If a null pointer then all edges will have unit capacity.
1697 * \return Error code.
1698 *
1699 * \sa \ref igraph_mincut(), \ref igraph_maxflow_value(), \ref
1700 * igraph_st_mincut_value().
1701 *
1702 * Time complexity: O(log(|V|)*|V|^2) for undirected graphs and
1703 * O(|V|^4) for directed graphs, but see also the discussion at the
1704 * documentation of \ref igraph_maxflow_value().
1705 */
1706
igraph_mincut_value(const igraph_t * graph,igraph_real_t * res,const igraph_vector_t * capacity)1707 int igraph_mincut_value(const igraph_t *graph, igraph_real_t *res,
1708 const igraph_vector_t *capacity) {
1709
1710 long int no_of_nodes = igraph_vcount(graph);
1711 igraph_real_t minmaxflow, flow;
1712 long int i;
1713
1714 minmaxflow = IGRAPH_INFINITY;
1715
1716 if (!igraph_is_directed(graph)) {
1717 IGRAPH_CHECK(igraph_i_mincut_value_undirected(graph, res, capacity));
1718 return 0;
1719 }
1720
1721 for (i = 1; i < no_of_nodes; i++) {
1722 IGRAPH_CHECK(igraph_maxflow_value(graph, &flow, 0, (igraph_integer_t) i,
1723 capacity, 0));
1724 if (flow < minmaxflow) {
1725 minmaxflow = flow;
1726 if (flow == 0) {
1727 break;
1728 }
1729 }
1730 IGRAPH_CHECK(igraph_maxflow_value(graph, &flow, (igraph_integer_t) i, 0,
1731 capacity, 0));
1732 if (flow < minmaxflow) {
1733 minmaxflow = flow;
1734 if (flow == 0) {
1735 break;
1736 }
1737 }
1738 }
1739
1740 if (res) {
1741 *res = minmaxflow;
1742 }
1743
1744 return 0;
1745 }
1746
igraph_i_st_vertex_connectivity_directed(const igraph_t * graph,igraph_integer_t * res,igraph_integer_t source,igraph_integer_t target,igraph_vconn_nei_t neighbors)1747 static int igraph_i_st_vertex_connectivity_directed(const igraph_t *graph,
1748 igraph_integer_t *res,
1749 igraph_integer_t source,
1750 igraph_integer_t target,
1751 igraph_vconn_nei_t neighbors) {
1752
1753 igraph_integer_t no_of_nodes = (igraph_integer_t) igraph_vcount(graph);
1754 igraph_integer_t no_of_edges = (igraph_integer_t) igraph_ecount(graph);
1755 igraph_vector_t edges;
1756 igraph_real_t real_res;
1757 igraph_t newgraph;
1758 long int i;
1759 igraph_bool_t conn1;
1760
1761 if (source < 0 || source >= no_of_nodes || target < 0 || target >= no_of_nodes) {
1762 IGRAPH_ERROR("Invalid source or target vertex", IGRAPH_EINVAL);
1763 }
1764
1765 switch (neighbors) {
1766 case IGRAPH_VCONN_NEI_ERROR:
1767 IGRAPH_CHECK(igraph_are_connected(graph, source, target, &conn1));
1768 if (conn1) {
1769 IGRAPH_ERROR("vertices connected", IGRAPH_EINVAL);
1770 return 0;
1771 }
1772 break;
1773 case IGRAPH_VCONN_NEI_NEGATIVE:
1774 IGRAPH_CHECK(igraph_are_connected(graph, source, target, &conn1));
1775 if (conn1) {
1776 *res = -1;
1777 return 0;
1778 }
1779 break;
1780 case IGRAPH_VCONN_NEI_NUMBER_OF_NODES:
1781 IGRAPH_CHECK(igraph_are_connected(graph, source, target, &conn1));
1782 if (conn1) {
1783 *res = no_of_nodes;
1784 return 0;
1785 }
1786 break;
1787 case IGRAPH_VCONN_NEI_IGNORE:
1788 break;
1789 default:
1790 IGRAPH_ERROR("Unknown `igraph_vconn_nei_t'", IGRAPH_EINVAL);
1791 break;
1792 }
1793
1794 /* Create the new graph */
1795
1796 IGRAPH_VECTOR_INIT_FINALLY(&edges, 0);
1797 IGRAPH_CHECK(igraph_vector_reserve(&edges, 2 * (no_of_edges + no_of_nodes)));
1798 IGRAPH_CHECK(igraph_get_edgelist(graph, &edges, 0));
1799 IGRAPH_CHECK(igraph_vector_resize(&edges, 2 * (no_of_edges + no_of_nodes)));
1800
1801 for (i = 0; i < 2 * no_of_edges; i += 2) {
1802 igraph_integer_t to = (igraph_integer_t) VECTOR(edges)[i + 1];
1803 if (to != source && to != target) {
1804 VECTOR(edges)[i + 1] = no_of_nodes + to;
1805 }
1806 }
1807
1808 for (i = 0; i < no_of_nodes; i++) {
1809 VECTOR(edges)[ 2 * (no_of_edges + i) ] = no_of_nodes + i;
1810 VECTOR(edges)[ 2 * (no_of_edges + i) + 1 ] = i;
1811 }
1812
1813 IGRAPH_CHECK(igraph_create(&newgraph, &edges, 2 * no_of_nodes,
1814 igraph_is_directed(graph)));
1815
1816 igraph_vector_destroy(&edges);
1817 IGRAPH_FINALLY_CLEAN(1);
1818 IGRAPH_FINALLY(igraph_destroy, &newgraph);
1819
1820 /* Do the maximum flow */
1821
1822 no_of_nodes = igraph_vcount(&newgraph);
1823 no_of_edges = igraph_ecount(&newgraph);
1824
1825 IGRAPH_CHECK(igraph_maxflow_value(&newgraph, &real_res,
1826 source, target, 0, 0));
1827 *res = (igraph_integer_t)real_res;
1828
1829 igraph_destroy(&newgraph);
1830 IGRAPH_FINALLY_CLEAN(1);
1831
1832 return 0;
1833 }
1834
igraph_i_st_vertex_connectivity_undirected(const igraph_t * graph,igraph_integer_t * res,igraph_integer_t source,igraph_integer_t target,igraph_vconn_nei_t neighbors)1835 static int igraph_i_st_vertex_connectivity_undirected(const igraph_t *graph,
1836 igraph_integer_t *res,
1837 igraph_integer_t source,
1838 igraph_integer_t target,
1839 igraph_vconn_nei_t neighbors) {
1840
1841 igraph_integer_t no_of_nodes = (igraph_integer_t) igraph_vcount(graph);
1842 igraph_t newgraph;
1843 igraph_bool_t conn;
1844
1845 if (source < 0 || source >= no_of_nodes || target < 0 || target >= no_of_nodes) {
1846 IGRAPH_ERROR("Invalid source or target vertex", IGRAPH_EINVAL);
1847 }
1848
1849 switch (neighbors) {
1850 case IGRAPH_VCONN_NEI_ERROR:
1851 IGRAPH_CHECK(igraph_are_connected(graph, source, target, &conn));
1852 if (conn) {
1853 IGRAPH_ERROR("vertices connected", IGRAPH_EINVAL);
1854 return 0;
1855 }
1856 break;
1857 case IGRAPH_VCONN_NEI_NEGATIVE:
1858 IGRAPH_CHECK(igraph_are_connected(graph, source, target, &conn));
1859 if (conn) {
1860 *res = -1;
1861 return 0;
1862 }
1863 break;
1864 case IGRAPH_VCONN_NEI_NUMBER_OF_NODES:
1865 IGRAPH_CHECK(igraph_are_connected(graph, source, target, &conn));
1866 if (conn) {
1867 *res = no_of_nodes;
1868 return 0;
1869 }
1870 break;
1871 case IGRAPH_VCONN_NEI_IGNORE:
1872 break;
1873 default:
1874 IGRAPH_ERROR("Unknown `igraph_vconn_nei_t'", IGRAPH_EINVAL);
1875 break;
1876 }
1877
1878 IGRAPH_CHECK(igraph_copy(&newgraph, graph));
1879 IGRAPH_FINALLY(igraph_destroy, &newgraph);
1880 IGRAPH_CHECK(igraph_to_directed(&newgraph, IGRAPH_TO_DIRECTED_MUTUAL));
1881
1882 IGRAPH_CHECK(igraph_i_st_vertex_connectivity_directed(&newgraph, res,
1883 source, target,
1884 IGRAPH_VCONN_NEI_IGNORE));
1885
1886 igraph_destroy(&newgraph);
1887 IGRAPH_FINALLY_CLEAN(1);
1888
1889 return 0;
1890 }
1891
1892 /**
1893 * \function igraph_st_vertex_connectivity
1894 * \brief The vertex connectivity of a pair of vertices
1895 *
1896 * </para><para>The vertex connectivity of two vertices (\c source and
1897 * \c target) is the minimum number of vertices that have to be
1898 * deleted to eliminate all paths from \c source to \c
1899 * target. Directed paths are considered in directed graphs.</para>
1900 *
1901 * <para>The vertex connectivity of a pair is the same as the number
1902 * of different (i.e. node-independent) paths from source to
1903 * target.</para>
1904 *
1905 * <para>The current implementation uses maximum flow calculations to
1906 * obtain the result.
1907 * \param graph The input graph.
1908 * \param res Pointer to an integer, the result will be stored here.
1909 * \param source The id of the source vertex.
1910 * \param target The id of the target vertex.
1911 * \param neighbors A constant giving what to do if the two vertices
1912 * are connected. Possible values:
1913 * \c IGRAPH_VCONN_NEI_ERROR, stop with an error message,
1914 * \c IGRAPH_VCONN_NEGATIVE, return -1.
1915 * \c IGRAPH_VCONN_NUMBER_OF_NODES, return the number of nodes.
1916 * \c IGRAPH_VCONN_IGNORE, ignore the fact that the two vertices
1917 * are connected and calculated the number of vertices needed
1918 * to eliminate all paths except for the trivial (direct) paths
1919 * between \c source and \c vertex. TOOD: what about neighbors?
1920 * \return Error code.
1921 *
1922 * Time complexity: O(|V|^3), but see the discussion at \ref
1923 * igraph_maxflow_value().
1924 *
1925 * \sa \ref igraph_vertex_connectivity(),
1926 * \ref igraph_edge_connectivity(),
1927 * \ref igraph_maxflow_value().
1928 */
1929
igraph_st_vertex_connectivity(const igraph_t * graph,igraph_integer_t * res,igraph_integer_t source,igraph_integer_t target,igraph_vconn_nei_t neighbors)1930 int igraph_st_vertex_connectivity(const igraph_t *graph,
1931 igraph_integer_t *res,
1932 igraph_integer_t source,
1933 igraph_integer_t target,
1934 igraph_vconn_nei_t neighbors) {
1935
1936 if (source == target) {
1937 IGRAPH_ERROR("source and target vertices are the same", IGRAPH_EINVAL);
1938 }
1939
1940 if (igraph_is_directed(graph)) {
1941 IGRAPH_CHECK(igraph_i_st_vertex_connectivity_directed(graph, res,
1942 source, target,
1943 neighbors));
1944 } else {
1945 IGRAPH_CHECK(igraph_i_st_vertex_connectivity_undirected(graph, res,
1946 source, target,
1947 neighbors));
1948 }
1949
1950 return 0;
1951 }
1952
igraph_i_vertex_connectivity_directed(const igraph_t * graph,igraph_integer_t * res)1953 static int igraph_i_vertex_connectivity_directed(const igraph_t *graph,
1954 igraph_integer_t *res) {
1955
1956 igraph_integer_t no_of_nodes = (igraph_integer_t) igraph_vcount(graph);
1957 long int i, j;
1958 igraph_integer_t minconn = no_of_nodes - 1, conn;
1959
1960 for (i = 0; i < no_of_nodes; i++) {
1961 for (j = 0; j < no_of_nodes; j++) {
1962 if (i == j) {
1963 continue;
1964 }
1965
1966 IGRAPH_ALLOW_INTERRUPTION();
1967
1968 IGRAPH_CHECK(igraph_st_vertex_connectivity(graph, &conn,
1969 (igraph_integer_t) i,
1970 (igraph_integer_t) j,
1971 IGRAPH_VCONN_NEI_NUMBER_OF_NODES));
1972 if (conn < minconn) {
1973 minconn = conn;
1974 if (conn == 0) {
1975 break;
1976 }
1977 }
1978 }
1979 if (conn == 0) {
1980 break;
1981 }
1982 }
1983
1984 if (res) {
1985 *res = minconn;
1986 }
1987
1988 return 0;
1989 }
1990
igraph_i_vertex_connectivity_undirected(const igraph_t * graph,igraph_integer_t * res)1991 static int igraph_i_vertex_connectivity_undirected(const igraph_t *graph,
1992 igraph_integer_t *res) {
1993 igraph_t newgraph;
1994
1995 IGRAPH_CHECK(igraph_copy(&newgraph, graph));
1996 IGRAPH_FINALLY(igraph_destroy, &newgraph);
1997 IGRAPH_CHECK(igraph_to_directed(&newgraph, IGRAPH_TO_DIRECTED_MUTUAL));
1998
1999 IGRAPH_CHECK(igraph_i_vertex_connectivity_directed(&newgraph, res));
2000
2001 igraph_destroy(&newgraph);
2002 IGRAPH_FINALLY_CLEAN(1);
2003
2004 return 0;
2005 }
2006
2007 /* Use that vertex.connectivity(G) <= edge.connectivity(G) <= min(degree(G)) */
igraph_i_connectivity_checks(const igraph_t * graph,igraph_integer_t * res,igraph_bool_t * found)2008 static int igraph_i_connectivity_checks(const igraph_t *graph,
2009 igraph_integer_t *res,
2010 igraph_bool_t *found) {
2011 igraph_bool_t conn;
2012 *found = 0;
2013
2014 if (igraph_vcount(graph) == 0) {
2015 *res = 0;
2016 *found = 1;
2017 return 0;
2018 }
2019
2020 IGRAPH_CHECK(igraph_is_connected(graph, &conn, IGRAPH_STRONG));
2021 if (!conn) {
2022 *res = 0;
2023 *found = 1;
2024 } else {
2025 igraph_vector_t degree;
2026 IGRAPH_VECTOR_INIT_FINALLY(°ree, 0);
2027 if (!igraph_is_directed(graph)) {
2028 IGRAPH_CHECK(igraph_degree(graph, °ree, igraph_vss_all(),
2029 IGRAPH_OUT, IGRAPH_LOOPS));
2030 if (igraph_vector_min(°ree) == 1) {
2031 *res = 1;
2032 *found = 1;
2033 }
2034 } else {
2035 /* directed, check both in- & out-degree */
2036 IGRAPH_CHECK(igraph_degree(graph, °ree, igraph_vss_all(),
2037 IGRAPH_OUT, IGRAPH_LOOPS));
2038 if (igraph_vector_min(°ree) == 1) {
2039 *res = 1;
2040 *found = 1;
2041 } else {
2042 IGRAPH_CHECK(igraph_degree(graph, °ree, igraph_vss_all(),
2043 IGRAPH_IN, IGRAPH_LOOPS));
2044 if (igraph_vector_min(°ree) == 1) {
2045 *res = 1;
2046 *found = 1;
2047 }
2048 }
2049 }
2050 igraph_vector_destroy(°ree);
2051 IGRAPH_FINALLY_CLEAN(1);
2052 }
2053 return 0;
2054 }
2055
2056 /**
2057 * \function igraph_vertex_connectivity
2058 * The vertex connectivity of a graph
2059 *
2060 * </para><para> The vertex connectivity of a graph is the minimum
2061 * vertex connectivity along each pairs of vertices in the graph.
2062 * </para>
2063 * <para> The vertex connectivity of a graph is the same as group
2064 * cohesion as defined in Douglas R. White and Frank Harary: The
2065 * cohesiveness of blocks in social networks: node connectivity and
2066 * conditional density, Sociological Methodology 31:305--359, 2001.
2067 * \param graph The input graph.
2068 * \param res Pointer to an integer, the result will be stored here.
2069 * \param checks Logical constant. Whether to check that the graph is
2070 * connected and also the degree of the vertices. If the graph is
2071 * not (strongly) connected then the connectivity is obviously zero. Otherwise
2072 * if the minimum degree is one then the vertex connectivity is also
2073 * one. It is a good idea to perform these checks, as they can be
2074 * done quickly compared to the connectivity calculation itself.
2075 * They were suggested by Peter McMahan, thanks Peter.
2076 * \return Error code.
2077 *
2078 * Time complexity: O(|V|^5).
2079 *
2080 * \sa \ref igraph_st_vertex_connectivity(), \ref igraph_maxflow_value(),
2081 * and \ref igraph_edge_connectivity().
2082 */
2083
igraph_vertex_connectivity(const igraph_t * graph,igraph_integer_t * res,igraph_bool_t checks)2084 int igraph_vertex_connectivity(const igraph_t *graph, igraph_integer_t *res,
2085 igraph_bool_t checks) {
2086
2087 igraph_bool_t ret = 0;
2088
2089 if (checks) {
2090 IGRAPH_CHECK(igraph_i_connectivity_checks(graph, res, &ret));
2091 }
2092
2093 /* Are we done yet? */
2094 if (!ret) {
2095 if (igraph_is_directed(graph)) {
2096 IGRAPH_CHECK(igraph_i_vertex_connectivity_directed(graph, res));
2097 } else {
2098 IGRAPH_CHECK(igraph_i_vertex_connectivity_undirected(graph, res));
2099 }
2100 }
2101
2102 return 0;
2103 }
2104
2105 /**
2106 * \function igraph_st_edge_connectivity
2107 * \brief Edge connectivity of a pair of vertices
2108 *
2109 * </para><para> The edge connectivity of two vertices (\c source and
2110 * \c target) in a graph is the minimum number of edges that
2111 * have to be deleted from the graph to eliminate all paths from \c
2112 * source to \c target.</para>
2113 *
2114 * <para>This function uses the maximum flow algorithm to calculate
2115 * the edge connectivity.
2116 * \param graph The input graph, it has to be directed.
2117 * \param res Pointer to an integer, the result will be stored here.
2118 * \param source The id of the source vertex.
2119 * \param target The id of the target vertex.
2120 * \return Error code.
2121 *
2122 * Time complexity: O(|V|^3).
2123 *
2124 * \sa \ref igraph_maxflow_value(), \ref igraph_edge_connectivity(),
2125 * \ref igraph_st_vertex_connectivity(), \ref
2126 * igraph_vertex_connectivity().
2127 */
2128
igraph_st_edge_connectivity(const igraph_t * graph,igraph_integer_t * res,igraph_integer_t source,igraph_integer_t target)2129 int igraph_st_edge_connectivity(const igraph_t *graph, igraph_integer_t *res,
2130 igraph_integer_t source,
2131 igraph_integer_t target) {
2132 igraph_real_t flow;
2133
2134 if (source == target) {
2135 IGRAPH_ERROR("source and target vertices are the same", IGRAPH_EINVAL);
2136 }
2137
2138 IGRAPH_CHECK(igraph_maxflow_value(graph, &flow, source, target, 0, 0));
2139 *res = (igraph_integer_t) flow;
2140
2141 return 0;
2142 }
2143
2144
2145 /**
2146 * \function igraph_edge_connectivity
2147 * \brief The minimum edge connectivity in a graph.
2148 *
2149 * </para><para> This is the minimum of the edge connectivity over all
2150 * pairs of vertices in the graph. </para>
2151 *
2152 * <para>
2153 * The edge connectivity of a graph is the same as group adhesion as
2154 * defined in Douglas R. White and Frank Harary: The cohesiveness of
2155 * blocks in social networks: node connectivity and conditional
2156 * density, Sociological Methodology 31:305--359, 2001.
2157 * \param graph The input graph.
2158 * \param res Pointer to an integer, the result will be stored here.
2159 * \param checks Logical constant. Whether to check that the graph is
2160 * connected and also the degree of the vertices. If the graph is
2161 * not (strongly) connected then the connectivity is obviously zero. Otherwise
2162 * if the minimum degree is one then the edge connectivity is also
2163 * one. It is a good idea to perform these checks, as they can be
2164 * done quickly compared to the connectivity calculation itself.
2165 * They were suggested by Peter McMahan, thanks Peter.
2166 * \return Error code.
2167 *
2168 * Time complexity: O(log(|V|)*|V|^2) for undirected graphs and
2169 * O(|V|^4) for directed graphs, but see also the discussion at the
2170 * documentation of \ref igraph_maxflow_value().
2171 *
2172 * \sa \ref igraph_st_edge_connectivity(), \ref igraph_maxflow_value(),
2173 * \ref igraph_vertex_connectivity().
2174 */
2175
igraph_edge_connectivity(const igraph_t * graph,igraph_integer_t * res,igraph_bool_t checks)2176 int igraph_edge_connectivity(const igraph_t *graph, igraph_integer_t *res,
2177 igraph_bool_t checks) {
2178 igraph_bool_t ret = 0;
2179 igraph_integer_t number_of_nodes = igraph_vcount(graph);
2180
2181 /* igraph_mincut_value returns infinity for the singleton graph,
2182 * which cannot be cast to an integer. We catch this case early
2183 * and postulate the edge-connectivity of this graph to be 0.
2184 * This is consistent with what other software packages return. */
2185 if (number_of_nodes <= 1) {
2186 *res = 0;
2187 return 0;
2188 }
2189
2190 /* Use that vertex.connectivity(G) <= edge.connectivity(G) <= min(degree(G)) */
2191 if (checks) {
2192 IGRAPH_CHECK(igraph_i_connectivity_checks(graph, res, &ret));
2193 }
2194
2195 if (!ret) {
2196 igraph_real_t real_res;
2197 IGRAPH_CHECK(igraph_mincut_value(graph, &real_res, 0));
2198 *res = (igraph_integer_t)real_res;
2199 }
2200
2201 return 0;
2202 }
2203
2204 /**
2205 * \function igraph_edge_disjoint_paths
2206 * \brief The maximum number of edge-disjoint paths between two vertices.
2207 *
2208 * </para><para> A set of paths between two vertices is called
2209 * edge-disjoint if they do not share any edges. The maximum number of
2210 * edge-disjoint paths are calculated by this function using maximum
2211 * flow techniques. Directed paths are considered in directed
2212 * graphs. </para>
2213 *
2214 * <para> Note that the number of disjoint paths is the same as the
2215 * edge connectivity of the two vertices using uniform edge weights.
2216 * \param graph The input graph, can be directed or undirected.
2217 * \param res Pointer to an integer variable, the result will be
2218 * stored here.
2219 * \param source The id of the source vertex.
2220 * \param target The id of the target vertex.
2221 * \return Error code.
2222 *
2223 * Time complexity: O(|V|^3), but see the discussion at \ref
2224 * igraph_maxflow_value().
2225 *
2226 * \sa \ref igraph_vertex_disjoint_paths(), \ref
2227 * igraph_st_edge_connectivity(), \ref igraph_maxflow_value().
2228 */
2229
igraph_edge_disjoint_paths(const igraph_t * graph,igraph_integer_t * res,igraph_integer_t source,igraph_integer_t target)2230 int igraph_edge_disjoint_paths(const igraph_t *graph, igraph_integer_t *res,
2231 igraph_integer_t source,
2232 igraph_integer_t target) {
2233
2234 igraph_real_t flow;
2235
2236 if (source == target) {
2237 IGRAPH_ERROR("Not implemented for source=target", IGRAPH_UNIMPLEMENTED);
2238 }
2239
2240 IGRAPH_CHECK(igraph_maxflow_value(graph, &flow, source, target, 0, 0));
2241
2242 *res = (igraph_integer_t) flow;
2243
2244 return 0;
2245 }
2246
2247 /**
2248 * \function igraph_vertex_disjoint_paths
2249 * \brief Maximum number of vertex-disjoint paths between two vertices.
2250 *
2251 * </para><para> A set of paths between two vertices is called
2252 * vertex-disjoint if they share no vertices. The calculation is
2253 * performed by using maximum flow techniques. </para>
2254 *
2255 * <para> Note that the number of vertex-disjoint paths is the same as
2256 * the vertex connectivity of the two vertices in most cases (if the
2257 * two vertices are not connected by an edge).
2258 * \param graph The input graph.
2259 * \param res Pointer to an integer variable, the result will be
2260 * stored here.
2261 * \param source The id of the source vertex.
2262 * \param target The id of the target vertex.
2263 * \return Error code.
2264 *
2265 * Time complexity: O(|V|^3).
2266 *
2267 * \sa \ref igraph_edge_disjoint_paths(), \ref
2268 * igraph_vertex_connectivity(), \ref igraph_maxflow_value().
2269 */
2270
igraph_vertex_disjoint_paths(const igraph_t * graph,igraph_integer_t * res,igraph_integer_t source,igraph_integer_t target)2271 int igraph_vertex_disjoint_paths(const igraph_t *graph, igraph_integer_t *res,
2272 igraph_integer_t source,
2273 igraph_integer_t target) {
2274
2275 igraph_bool_t conn;
2276
2277 if (source == target) {
2278 IGRAPH_ERROR("The source==target case is not implemented",
2279 IGRAPH_UNIMPLEMENTED);
2280 }
2281
2282 igraph_are_connected(graph, source, target, &conn);
2283 if (conn) {
2284 /* We need to remove every (possibly directed) edge between source
2285 and target and calculate the disjoint paths on the new
2286 graph. Finally we add 1 for the removed connection(s). */
2287 igraph_es_t es;
2288 igraph_vector_t v;
2289 igraph_t newgraph;
2290 IGRAPH_VECTOR_INIT_FINALLY(&v, 2);
2291 VECTOR(v)[0] = source;
2292 VECTOR(v)[1] = target;
2293 IGRAPH_CHECK(igraph_es_multipairs(&es, &v, IGRAPH_DIRECTED));
2294 IGRAPH_FINALLY(igraph_es_destroy, &es);
2295
2296 IGRAPH_CHECK(igraph_copy(&newgraph, graph));
2297 IGRAPH_FINALLY(igraph_destroy, &newgraph);
2298 IGRAPH_CHECK(igraph_delete_edges(&newgraph, es));
2299
2300 if (igraph_is_directed(graph)) {
2301 IGRAPH_CHECK(igraph_i_st_vertex_connectivity_directed(&newgraph, res,
2302 source, target,
2303 IGRAPH_VCONN_NEI_IGNORE));
2304 } else {
2305 IGRAPH_CHECK(igraph_i_st_vertex_connectivity_undirected(&newgraph, res,
2306 source, target,
2307 IGRAPH_VCONN_NEI_IGNORE));
2308 }
2309
2310 if (res) {
2311 *res += 1;
2312 }
2313
2314 IGRAPH_FINALLY_CLEAN(3);
2315 igraph_destroy(&newgraph);
2316 igraph_es_destroy(&es);
2317 igraph_vector_destroy(&v);
2318 }
2319
2320 /* These do nothing if the two vertices are connected,
2321 so it is safe to call them. */
2322
2323 if (igraph_is_directed(graph)) {
2324 IGRAPH_CHECK(igraph_i_st_vertex_connectivity_directed(graph, res,
2325 source, target,
2326 IGRAPH_VCONN_NEI_IGNORE));
2327 } else {
2328 IGRAPH_CHECK(igraph_i_st_vertex_connectivity_undirected(graph, res,
2329 source, target,
2330 IGRAPH_VCONN_NEI_IGNORE));
2331 }
2332
2333 return 0;
2334 }
2335
2336 /**
2337 * \function igraph_adhesion
2338 * \brief Graph adhesion, this is (almost) the same as edge connectivity.
2339 *
2340 * </para><para> This quantity is defined by White and Harary in
2341 * The cohesiveness of blocks in social networks: node connectivity and
2342 * conditional density, (Sociological Methodology 31:305--359, 2001)
2343 * and basically it is the edge connectivity of the graph
2344 * with uniform edge weights.
2345 * \param graph The input graph, either directed or undirected.
2346 * \param res Pointer to an integer, the result will be stored here.
2347 * \param checks Logical constant. Whether to check that the graph is
2348 * connected and also the degree of the vertices. If the graph is
2349 * not (strongly) connected then the adhesion is obviously zero. Otherwise
2350 * if the minimum degree is one then the adhesion is also
2351 * one. It is a good idea to perform these checks, as they can be
2352 * done quickly compared to the edge connectivity calculation itself.
2353 * They were suggested by Peter McMahan, thanks Peter.
2354 * \return Error code.
2355 *
2356 * Time complexity: O(log(|V|)*|V|^2) for undirected graphs and
2357 * O(|V|^4) for directed graphs, but see also the discussion at the
2358 * documentation of \ref igraph_maxflow_value().
2359 *
2360 * \sa \ref igraph_cohesion(), \ref igraph_maxflow_value(), \ref
2361 * igraph_edge_connectivity(), \ref igraph_mincut_value().
2362 */
2363
igraph_adhesion(const igraph_t * graph,igraph_integer_t * res,igraph_bool_t checks)2364 int igraph_adhesion(const igraph_t *graph, igraph_integer_t *res,
2365 igraph_bool_t checks) {
2366 return igraph_edge_connectivity(graph, res, checks);
2367 }
2368
2369 /**
2370 * \function igraph_cohesion
2371 * \brief Graph cohesion, this is the same as vertex connectivity.
2372 *
2373 * </para><para> This quantity was defined by White and Harary in <quote>The
2374 * cohesiveness of blocks in social networks: node connectivity and
2375 * conditional density</quote>, (Sociological Methodology 31:305--359, 2001)
2376 * and it is the same as the vertex connectivity of a
2377 * graph.
2378 * \param graph The input graph.
2379 * \param res Pointer to an integer variable, the result will be
2380 * stored here.
2381 * \param checks Logical constant. Whether to check that the graph is
2382 * connected and also the degree of the vertices. If the graph is
2383 * not (strongly) connected then the cohesion is obviously zero. Otherwise
2384 * if the minimum degree is one then the cohesion is also
2385 * one. It is a good idea to perform these checks, as they can be
2386 * done quickly compared to the vertex connectivity calculation itself.
2387 * They were suggested by Peter McMahan, thanks Peter.
2388 * \return Error code.
2389 *
2390 * Time complexity: O(|V|^4), |V| is the number of vertices. In
2391 * practice it is more like O(|V|^2), see \ref igraph_maxflow_value().
2392 *
2393 * \sa \ref igraph_vertex_connectivity(), \ref igraph_adhesion(),
2394 * \ref igraph_maxflow_value().
2395 */
2396
igraph_cohesion(const igraph_t * graph,igraph_integer_t * res,igraph_bool_t checks)2397 int igraph_cohesion(const igraph_t *graph, igraph_integer_t *res,
2398 igraph_bool_t checks) {
2399
2400 IGRAPH_CHECK(igraph_vertex_connectivity(graph, res, checks));
2401 return 0;
2402 }
2403
2404 /**
2405 * \function igraph_gomory_hu_tree
2406 * \brief Gomory-Hu tree of a graph.
2407 *
2408 * </para><para>
2409 * The Gomory-Hu tree is a concise representation of the value of all the
2410 * maximum flows (or minimum cuts) in a graph. The vertices of the tree
2411 * correspond exactly to the vertices of the original graph in the same order.
2412 * Edges of the Gomory-Hu tree are annotated by flow values. The value of
2413 * the maximum flow (or minimum cut) between an arbitrary (u,v) vertex
2414 * pair in the original graph is then given by the minimum flow value (i.e.
2415 * edge annotation) along the shortest path between u and v in the
2416 * Gomory-Hu tree.
2417 *
2418 * </para><para>This implementation uses Gusfield's algorithm to construct the
2419 * Gomory-Hu tree. See the following paper for more details:
2420 *
2421 * </para><para>
2422 * Gusfield D: Very simple methods for all pairs network flow analysis. SIAM J
2423 * Comput 19(1):143-155, 1990.
2424 *
2425 * \param graph The input graph.
2426 * \param tree Pointer to an uninitialized graph; the result will be
2427 * stored here.
2428 * \param flows Pointer to an uninitialized vector; the flow values
2429 * corresponding to each edge in the Gomory-Hu tree will
2430 * be returned here. You may pass a NULL pointer here if you are
2431 * not interested in the flow values.
2432 * \param capacity Vector containing the capacity of the edges. If NULL, then
2433 * every edge is considered to have capacity 1.0.
2434 * \return Error code.
2435 *
2436 * Time complexity: O(|V|^4) since it performs a max-flow calculation
2437 * between vertex zero and every other vertex and max-flow is
2438 * O(|V|^3).
2439 *
2440 * \sa \ref igraph_maxflow()
2441 */
igraph_gomory_hu_tree(const igraph_t * graph,igraph_t * tree,igraph_vector_t * flows,const igraph_vector_t * capacity)2442 int igraph_gomory_hu_tree(const igraph_t *graph, igraph_t *tree,
2443 igraph_vector_t *flows, const igraph_vector_t *capacity) {
2444
2445 igraph_integer_t no_of_nodes = igraph_vcount(graph);
2446 igraph_integer_t source, target, mid, i, n;
2447 igraph_vector_t neighbors;
2448 igraph_vector_t flow_values;
2449 igraph_vector_t partition;
2450 igraph_vector_t partition2;
2451 igraph_real_t flow_value;
2452
2453 if (igraph_is_directed(graph)) {
2454 IGRAPH_ERROR("Gomory-Hu tree can only be calculated for undirected graphs",
2455 IGRAPH_EINVAL);
2456 }
2457
2458 /* Allocate memory */
2459 IGRAPH_VECTOR_INIT_FINALLY(&neighbors, no_of_nodes);
2460 IGRAPH_VECTOR_INIT_FINALLY(&flow_values, no_of_nodes);
2461 IGRAPH_VECTOR_INIT_FINALLY(&partition, 0);
2462 IGRAPH_VECTOR_INIT_FINALLY(&partition2, 0);
2463
2464 /* Initialize the tree: every edge points to node 0 */
2465 /* Actually, this is done implicitly since both 'neighbors' and 'flow_values' are
2466 * initialized to zero already */
2467
2468 /* For each source vertex except vertex zero... */
2469 for (source = 1; source < no_of_nodes; source++) {
2470 IGRAPH_ALLOW_INTERRUPTION();
2471 IGRAPH_PROGRESS("Gomory-Hu tree", (100.0 * (source - 1)) / (no_of_nodes - 1), 0);
2472
2473 /* Find its current neighbor in the tree */
2474 target = VECTOR(neighbors)[(long int)source];
2475
2476 /* Find the maximum flow between source and target */
2477 IGRAPH_CHECK(igraph_maxflow(graph, &flow_value, 0, 0, &partition, &partition2,
2478 source, target, capacity, 0));
2479
2480 /* Store the maximum flow and determine which side each node is on */
2481 VECTOR(flow_values)[(long int)source] = flow_value;
2482
2483 /* Update the tree */
2484 /* igraph_maxflow() guarantees that the source vertex will be in &partition
2485 * and not in &partition2 */
2486 n = igraph_vector_size(&partition);
2487 for (i = 0; i < n; i++) {
2488 mid = VECTOR(partition)[i];
2489 if (mid > source && VECTOR(neighbors)[(long int)mid] == target) {
2490 VECTOR(neighbors)[(long int)mid] = source;
2491 }
2492 }
2493 }
2494
2495 IGRAPH_PROGRESS("Gomory-Hu tree", 100.0, 0);
2496
2497 /* Re-use the 'partition' vector as an edge list now */
2498 IGRAPH_CHECK(igraph_vector_resize(&partition, 2 * (no_of_nodes - 1)));
2499 for (i = 1, mid = 0; i < no_of_nodes; i++, mid += 2) {
2500 VECTOR(partition)[(long int)mid] = i;
2501 VECTOR(partition)[(long int)mid + 1] = VECTOR(neighbors)[(long int)i];
2502 }
2503
2504 /* Create the tree graph; we use igraph_subgraph_edges here to keep the
2505 * graph and vertex attributes */
2506 IGRAPH_CHECK(igraph_subgraph_edges(graph, tree, igraph_ess_none(), 0));
2507 IGRAPH_CHECK(igraph_add_edges(tree, &partition, 0));
2508
2509 /* Free the allocated memory */
2510 igraph_vector_destroy(&partition2);
2511 igraph_vector_destroy(&partition);
2512 igraph_vector_destroy(&neighbors);
2513 IGRAPH_FINALLY_CLEAN(3);
2514
2515 /* Return the flow values to the caller */
2516 if (flows != 0) {
2517 IGRAPH_CHECK(igraph_vector_update(flows, &flow_values));
2518 if (no_of_nodes > 0) {
2519 igraph_vector_remove(flows, 0);
2520 }
2521 }
2522
2523 /* Free the remaining allocated memory */
2524 igraph_vector_destroy(&flow_values);
2525 IGRAPH_FINALLY_CLEAN(1);
2526
2527 return IGRAPH_SUCCESS;
2528 }
2529