Statistics
| Branch: | Revision:

iof-tools / networkxMiCe / networkx-master / networkx / algorithms / flow / boykovkolmogorov.py @ 5cef0f13

History | View | Annotate | Download (13.3 KB)

1
# boykovkolmogorov.py - Boykov Kolmogorov algorithm for maximum flow problems.
2
#
3
# Copyright 2016-2019 NetworkX developers.
4
#
5
# This file is part of NetworkX.
6
#
7
# NetworkX is distributed under a BSD license; see LICENSE.txt for more
8
# information.
9
#
10
# Author: Jordi Torrents <jordi.t21@gmail.com>
11
"""
12
Boykov-Kolmogorov algorithm for maximum flow problems.
13
"""
14
from collections import deque
15
from operator import itemgetter
16

    
17
import networkx as nx
18
from networkx.algorithms.flow.utils import build_residual_network
19

    
20
__all__ = ['boykov_kolmogorov']
21

    
22

    
23
def boykov_kolmogorov(G, s, t, capacity='capacity', residual=None,
24
                      value_only=False, cutoff=None):
25
    r"""Find a maximum single-commodity flow using Boykov-Kolmogorov algorithm.
26

27
    This function returns the residual network resulting after computing
28
    the maximum flow. See below for details about the conventions
29
    NetworkX uses for defining residual networks.
30

31
    This algorithm has worse case complexity $O(n^2 m |C|)$ for $n$ nodes, $m$
32
    edges, and $|C|$ the cost of the minimum cut [1]_. This implementation
33
    uses the marking heuristic defined in [2]_ which improves its running
34
    time in many practical problems.
35

36
    Parameters
37
    ----------
38
    G : NetworkX graph
39
        Edges of the graph are expected to have an attribute called
40
        'capacity'. If this attribute is not present, the edge is
41
        considered to have infinite capacity.
42

43
    s : node
44
        Source node for the flow.
45

46
    t : node
47
        Sink node for the flow.
48

49
    capacity : string
50
        Edges of the graph G are expected to have an attribute capacity
51
        that indicates how much flow the edge can support. If this
52
        attribute is not present, the edge is considered to have
53
        infinite capacity. Default value: 'capacity'.
54

55
    residual : NetworkX graph
56
        Residual network on which the algorithm is to be executed. If None, a
57
        new residual network is created. Default value: None.
58

59
    value_only : bool
60
        If True compute only the value of the maximum flow. This parameter
61
        will be ignored by this algorithm because it is not applicable.
62

63
    cutoff : integer, float
64
        If specified, the algorithm will terminate when the flow value reaches
65
        or exceeds the cutoff. In this case, it may be unable to immediately
66
        determine a minimum cut. Default value: None.
67

68
    Returns
69
    -------
70
    R : NetworkX DiGraph
71
        Residual network after computing the maximum flow.
72

73
    Raises
74
    ------
75
    NetworkXError
76
        The algorithm does not support MultiGraph and MultiDiGraph. If
77
        the input graph is an instance of one of these two classes, a
78
        NetworkXError is raised.
79

80
    NetworkXUnbounded
81
        If the graph has a path of infinite capacity, the value of a
82
        feasible flow on the graph is unbounded above and the function
83
        raises a NetworkXUnbounded.
84

85
    See also
86
    --------
87
    :meth:`maximum_flow`
88
    :meth:`minimum_cut`
89
    :meth:`preflow_push`
90
    :meth:`shortest_augmenting_path`
91

92
    Notes
93
    -----
94
    The residual network :samp:`R` from an input graph :samp:`G` has the
95
    same nodes as :samp:`G`. :samp:`R` is a DiGraph that contains a pair
96
    of edges :samp:`(u, v)` and :samp:`(v, u)` iff :samp:`(u, v)` is not a
97
    self-loop, and at least one of :samp:`(u, v)` and :samp:`(v, u)` exists
98
    in :samp:`G`.
99

100
    For each edge :samp:`(u, v)` in :samp:`R`, :samp:`R[u][v]['capacity']`
101
    is equal to the capacity of :samp:`(u, v)` in :samp:`G` if it exists
102
    in :samp:`G` or zero otherwise. If the capacity is infinite,
103
    :samp:`R[u][v]['capacity']` will have a high arbitrary finite value
104
    that does not affect the solution of the problem. This value is stored in
105
    :samp:`R.graph['inf']`. For each edge :samp:`(u, v)` in :samp:`R`,
106
    :samp:`R[u][v]['flow']` represents the flow function of :samp:`(u, v)` and
107
    satisfies :samp:`R[u][v]['flow'] == -R[v][u]['flow']`.
108

109
    The flow value, defined as the total flow into :samp:`t`, the sink, is
110
    stored in :samp:`R.graph['flow_value']`. If :samp:`cutoff` is not
111
    specified, reachability to :samp:`t` using only edges :samp:`(u, v)` such
112
    that :samp:`R[u][v]['flow'] < R[u][v]['capacity']` induces a minimum
113
    :samp:`s`-:samp:`t` cut.
114

115
    Examples
116
    --------
117
    >>> import networkx as nx
118
    >>> from networkx.algorithms.flow import boykov_kolmogorov
119

120
    The functions that implement flow algorithms and output a residual
121
    network, such as this one, are not imported to the base NetworkX
122
    namespace, so you have to explicitly import them from the flow package.
123

124
    >>> G = nx.DiGraph()
125
    >>> G.add_edge('x','a', capacity=3.0)
126
    >>> G.add_edge('x','b', capacity=1.0)
127
    >>> G.add_edge('a','c', capacity=3.0)
128
    >>> G.add_edge('b','c', capacity=5.0)
129
    >>> G.add_edge('b','d', capacity=4.0)
130
    >>> G.add_edge('d','e', capacity=2.0)
131
    >>> G.add_edge('c','y', capacity=2.0)
132
    >>> G.add_edge('e','y', capacity=3.0)
133
    >>> R = boykov_kolmogorov(G, 'x', 'y')
134
    >>> flow_value = nx.maximum_flow_value(G, 'x', 'y')
135
    >>> flow_value
136
    3.0
137
    >>> flow_value == R.graph['flow_value']
138
    True
139

140
    A nice feature of the Boykov-Kolmogorov algorithm is that a partition
141
    of the nodes that defines a minimum cut can be easily computed based
142
    on the search trees used during the algorithm. These trees are stored
143
    in the graph attribute `trees` of the residual network.
144

145
    >>> source_tree, target_tree = R.graph['trees']
146
    >>> partition = (set(source_tree), set(G) - set(source_tree))
147

148
    Or equivalently:
149

150
    >>> partition = (set(G) - set(target_tree), set(target_tree))
151

152
    References
153
    ----------
154
    .. [1] Boykov, Y., & Kolmogorov, V. (2004). An experimental comparison
155
           of min-cut/max-flow algorithms for energy minimization in vision.
156
           Pattern Analysis and Machine Intelligence, IEEE Transactions on,
157
           26(9), 1124-1137.
158
           http://www.csd.uwo.ca/~yuri/Papers/pami04.pdf
159

160
    .. [2] Vladimir Kolmogorov. Graph-based Algorithms for Multi-camera
161
           Reconstruction Problem. PhD thesis, Cornell University, CS Department,
162
           2003. pp. 109-114.
163
           https://pub.ist.ac.at/~vnk/papers/thesis.pdf
164

165
    """
166
    R = boykov_kolmogorov_impl(G, s, t, capacity, residual, cutoff)
167
    R.graph['algorithm'] = 'boykov_kolmogorov'
168
    return R
169

    
170

    
171
def boykov_kolmogorov_impl(G, s, t, capacity, residual, cutoff):
172
    if s not in G:
173
        raise nx.NetworkXError('node %s not in graph' % str(s))
174
    if t not in G:
175
        raise nx.NetworkXError('node %s not in graph' % str(t))
176
    if s == t:
177
        raise nx.NetworkXError('source and sink are the same node')
178

    
179
    if residual is None:
180
        R = build_residual_network(G, capacity)
181
    else:
182
        R = residual
183

    
184
    # Initialize/reset the residual network.
185
    # This is way too slow
186
    #nx.set_edge_attributes(R, 0, 'flow')
187
    for u in R:
188
        for e in R[u].values():
189
            e['flow'] = 0
190

    
191
    # Use an arbitrary high value as infinite. It is computed
192
    # when building the residual network.
193
    INF = R.graph['inf']
194

    
195
    if cutoff is None:
196
        cutoff = INF
197

    
198
    R_succ = R.succ
199
    R_pred = R.pred
200

    
201
    def grow():
202
        """Bidirectional breadth-first search for the growth stage.
203

204
           Returns a connecting edge, that is and edge that connects
205
           a node from the source search tree with a node from the
206
           target search tree.
207
           The first node in the connecting edge is always from the
208
           source tree and the last node from the target tree.
209
        """
210
        while active:
211
            u = active[0]
212
            if u in source_tree:
213
                this_tree = source_tree
214
                other_tree = target_tree
215
                neighbors = R_succ
216
            else:
217
                this_tree = target_tree
218
                other_tree = source_tree
219
                neighbors = R_pred
220
            for v, attr in neighbors[u].items():
221
                if attr['capacity'] - attr['flow'] > 0:
222
                    if v not in this_tree:
223
                        if v in other_tree:
224
                            return (u, v) if this_tree is source_tree else (v, u)
225
                        this_tree[v] = u
226
                        dist[v] = dist[u] + 1
227
                        timestamp[v] = timestamp[u]
228
                        active.append(v)
229
                    elif v in this_tree and _is_closer(u, v):
230
                        this_tree[v] = u
231
                        dist[v] = dist[u] + 1
232
                        timestamp[v] = timestamp[u]
233
            _ = active.popleft()
234
        return None, None
235

    
236
    def augment(u, v):
237
        """Augmentation stage.
238

239
           Reconstruct path and determine its residual capacity.
240
           We start from a connecting edge, which links a node
241
           from the source tree to a node from the target tree.
242
           The connecting edge is the output of the grow function
243
           and the input of this function.
244
        """
245
        attr = R_succ[u][v]
246
        flow = min(INF, attr['capacity'] - attr['flow'])
247
        path = [u]
248
        # Trace a path from u to s in source_tree.
249
        w = u
250
        while w != s:
251
            n = w
252
            w = source_tree[n]
253
            attr = R_pred[n][w]
254
            flow = min(flow, attr['capacity'] - attr['flow'])
255
            path.append(w)
256
        path.reverse()
257
        # Trace a path from v to t in target_tree.
258
        path.append(v)
259
        w = v
260
        while w != t:
261
            n = w
262
            w = target_tree[n]
263
            attr = R_succ[n][w]
264
            flow = min(flow, attr['capacity'] - attr['flow'])
265
            path.append(w)
266
        # Augment flow along the path and check for saturated edges.
267
        it = iter(path)
268
        u = next(it)
269
        these_orphans = []
270
        for v in it:
271
            R_succ[u][v]['flow'] += flow
272
            R_succ[v][u]['flow'] -= flow
273
            if R_succ[u][v]['flow'] == R_succ[u][v]['capacity']:
274
                if v in source_tree:
275
                    source_tree[v] = None
276
                    these_orphans.append(v)
277
                if u in target_tree:
278
                    target_tree[u] = None
279
                    these_orphans.append(u)
280
            u = v
281
        orphans.extend(sorted(these_orphans, key=dist.get))
282
        return flow
283

    
284
    def adopt():
285
        """Adoption stage.
286

287
           Reconstruct search trees by adopting or discarding orphans.
288
           During augmentation stage some edges got saturated and thus
289
           the source and target search trees broke down to forests, with
290
           orphans as roots of some of its trees. We have to reconstruct
291
           the search trees rooted to source and target before we can grow
292
           them again.
293
        """
294
        while orphans:
295
            u = orphans.popleft()
296
            if u in source_tree:
297
                tree = source_tree
298
                neighbors = R_pred
299
            else:
300
                tree = target_tree
301
                neighbors = R_succ
302
            nbrs = ((n, attr, dist[n]) for n, attr in neighbors[u].items()
303
                    if n in tree)
304
            for v, attr, d in sorted(nbrs, key=itemgetter(2)):
305
                if attr['capacity'] - attr['flow'] > 0:
306
                    if _has_valid_root(v, tree):
307
                        tree[u] = v
308
                        dist[u] = dist[v] + 1
309
                        timestamp[u] = time
310
                        break
311
            else:
312
                nbrs = ((n, attr, dist[n]) for n, attr in neighbors[u].items()
313
                        if n in tree)
314
                for v, attr, d in sorted(nbrs, key=itemgetter(2)):
315
                    if attr['capacity'] - attr['flow'] > 0:
316
                        if v not in active:
317
                            active.append(v)
318
                    if tree[v] == u:
319
                        tree[v] = None
320
                        orphans.appendleft(v)
321
                if u in active:
322
                    active.remove(u)
323
                del tree[u]
324

    
325
    def _has_valid_root(n, tree):
326
        path = []
327
        v = n
328
        while v is not None:
329
            path.append(v)
330
            if v == s or v == t:
331
                base_dist = 0
332
                break
333
            elif timestamp[v] == time:
334
                base_dist = dist[v]
335
                break
336
            v = tree[v]
337
        else:
338
            return False
339
        length = len(path)
340
        for i, u in enumerate(path, 1):
341
            dist[u] = base_dist + length - i
342
            timestamp[u] = time
343
        return True
344

    
345
    def _is_closer(u, v):
346
        return timestamp[v] <= timestamp[u] and dist[v] > dist[u] + 1
347

    
348
    source_tree = {s: None}
349
    target_tree = {t: None}
350
    active = deque([s, t])
351
    orphans = deque()
352
    flow_value = 0
353
    # data structures for the marking heuristic
354
    time = 1
355
    timestamp = {s: time, t: time}
356
    dist = {s: 0, t: 0}
357
    while flow_value < cutoff:
358
        # Growth stage
359
        u, v = grow()
360
        if u is None:
361
            break
362
        time += 1
363
        # Augmentation stage
364
        flow_value += augment(u, v)
365
        # Adoption stage
366
        adopt()
367

    
368
    if flow_value * 2 > INF:
369
        raise nx.NetworkXUnbounded('Infinite capacity path, flow unbounded above.')
370

    
371
    # Add source and target tree in a graph attribute.
372
    # A partition that defines a minimum cut can be directly
373
    # computed from the search trees as explained in the docstrings.
374
    R.graph['trees'] = (source_tree, target_tree)
375
    # Add the standard flow_value graph attribute.
376
    R.graph['flow_value'] = flow_value
377
    return R