Statistics
| Branch: | Revision:

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

History | View | Annotate | Download (15.6 KB)

1
# -*- coding: utf-8 -*-
2
"""
3
Highest-label preflow-push algorithm for maximum flow problems.
4
"""
5

    
6
__author__ = """ysitu <ysitu@users.noreply.github.com>"""
7
# Copyright (C) 2014 ysitu <ysitu@users.noreply.github.com>
8
# All rights reserved.
9
# BSD license.
10

    
11
from collections import deque
12
from itertools import islice
13
import networkx as nx
14
#from networkx.algorithms.flow.utils import *
15
from ...utils import arbitrary_element
16
from .utils import build_residual_network
17
from .utils import CurrentEdge
18
from .utils import detect_unboundedness
19
from .utils import GlobalRelabelThreshold
20
from .utils import Level
21

    
22
__all__ = ['preflow_push']
23

    
24

    
25
def preflow_push_impl(G, s, t, capacity, residual, global_relabel_freq,
26
                      value_only):
27
    """Implementation of the highest-label preflow-push algorithm.
28
    """
29
    if s not in G:
30
        raise nx.NetworkXError('node %s not in graph' % str(s))
31
    if t not in G:
32
        raise nx.NetworkXError('node %s not in graph' % str(t))
33
    if s == t:
34
        raise nx.NetworkXError('source and sink are the same node')
35

    
36
    if global_relabel_freq is None:
37
        global_relabel_freq = 0
38
    if global_relabel_freq < 0:
39
        raise nx.NetworkXError('global_relabel_freq must be nonnegative.')
40

    
41
    if residual is None:
42
        R = build_residual_network(G, capacity)
43
    else:
44
        R = residual
45

    
46
    detect_unboundedness(R, s, t)
47

    
48
    R_nodes = R.nodes
49
    R_pred = R.pred
50
    R_succ = R.succ
51

    
52
    # Initialize/reset the residual network.
53
    for u in R:
54
        R_nodes[u]['excess'] = 0
55
        for e in R_succ[u].values():
56
            e['flow'] = 0
57

    
58
    def reverse_bfs(src):
59
        """Perform a reverse breadth-first search from src in the residual
60
        network.
61
        """
62
        heights = {src: 0}
63
        q = deque([(src, 0)])
64
        while q:
65
            u, height = q.popleft()
66
            height += 1
67
            for v, attr in R_pred[u].items():
68
                if v not in heights and attr['flow'] < attr['capacity']:
69
                    heights[v] = height
70
                    q.append((v, height))
71
        return heights
72

    
73
    # Initialize heights of the nodes.
74
    heights = reverse_bfs(t)
75

    
76
    if s not in heights:
77
        # t is not reachable from s in the residual network. The maximum flow
78
        # must be zero.
79
        R.graph['flow_value'] = 0
80
        return R
81

    
82
    n = len(R)
83
    # max_height represents the height of the highest level below level n with
84
    # at least one active node.
85
    max_height = max(heights[u] for u in heights if u != s)
86
    heights[s] = n
87

    
88
    grt = GlobalRelabelThreshold(n, R.size(), global_relabel_freq)
89

    
90
    # Initialize heights and 'current edge' data structures of the nodes.
91
    for u in R:
92
        R_nodes[u]['height'] = heights[u] if u in heights else n + 1
93
        R_nodes[u]['curr_edge'] = CurrentEdge(R_succ[u])
94

    
95
    def push(u, v, flow):
96
        """Push flow units of flow from u to v.
97
        """
98
        R_succ[u][v]['flow'] += flow
99
        R_succ[v][u]['flow'] -= flow
100
        R_nodes[u]['excess'] -= flow
101
        R_nodes[v]['excess'] += flow
102

    
103
    # The maximum flow must be nonzero now. Initialize the preflow by
104
    # saturating all edges emanating from s.
105
    for u, attr in R_succ[s].items():
106
        flow = attr['capacity']
107
        if flow > 0:
108
            push(s, u, flow)
109

    
110
    # Partition nodes into levels.
111
    levels = [Level() for i in range(2 * n)]
112
    for u in R:
113
        if u != s and u != t:
114
            level = levels[R_nodes[u]['height']]
115
            if R_nodes[u]['excess'] > 0:
116
                level.active.add(u)
117
            else:
118
                level.inactive.add(u)
119

    
120
    def activate(v):
121
        """Move a node from the inactive set to the active set of its level.
122
        """
123
        if v != s and v != t:
124
            level = levels[R_nodes[v]['height']]
125
            if v in level.inactive:
126
                level.inactive.remove(v)
127
                level.active.add(v)
128

    
129
    def relabel(u):
130
        """Relabel a node to create an admissible edge.
131
        """
132
        grt.add_work(len(R_succ[u]))
133
        return min(R_nodes[v]['height'] for v, attr in R_succ[u].items()
134
                   if attr['flow'] < attr['capacity']) + 1
135

    
136
    def discharge(u, is_phase1):
137
        """Discharge a node until it becomes inactive or, during phase 1 (see
138
        below), its height reaches at least n. The node is known to have the
139
        largest height among active nodes.
140
        """
141
        height = R_nodes[u]['height']
142
        curr_edge = R_nodes[u]['curr_edge']
143
        # next_height represents the next height to examine after discharging
144
        # the current node. During phase 1, it is capped to below n.
145
        next_height = height
146
        levels[height].active.remove(u)
147
        while True:
148
            v, attr = curr_edge.get()
149
            if (height == R_nodes[v]['height'] + 1 and
150
                    attr['flow'] < attr['capacity']):
151
                flow = min(R_nodes[u]['excess'],
152
                           attr['capacity'] - attr['flow'])
153
                push(u, v, flow)
154
                activate(v)
155
                if R_nodes[u]['excess'] == 0:
156
                    # The node has become inactive.
157
                    levels[height].inactive.add(u)
158
                    break
159
            try:
160
                curr_edge.move_to_next()
161
            except StopIteration:
162
                # We have run off the end of the adjacency list, and there can
163
                # be no more admissible edges. Relabel the node to create one.
164
                height = relabel(u)
165
                if is_phase1 and height >= n - 1:
166
                    # Although the node is still active, with a height at least
167
                    # n - 1, it is now known to be on the s side of the minimum
168
                    # s-t cut. Stop processing it until phase 2.
169
                    levels[height].active.add(u)
170
                    break
171
                # The first relabel operation after global relabeling may not
172
                # increase the height of the node since the 'current edge' data
173
                # structure is not rewound. Use height instead of (height - 1)
174
                # in case other active nodes at the same level are missed.
175
                next_height = height
176
        R_nodes[u]['height'] = height
177
        return next_height
178

    
179
    def gap_heuristic(height):
180
        """Apply the gap heuristic.
181
        """
182
        # Move all nodes at levels (height + 1) to max_height to level n + 1.
183
        for level in islice(levels, height + 1, max_height + 1):
184
            for u in level.active:
185
                R_nodes[u]['height'] = n + 1
186
            for u in level.inactive:
187
                R_nodes[u]['height'] = n + 1
188
            levels[n + 1].active.update(level.active)
189
            level.active.clear()
190
            levels[n + 1].inactive.update(level.inactive)
191
            level.inactive.clear()
192

    
193
    def global_relabel(from_sink):
194
        """Apply the global relabeling heuristic.
195
        """
196
        src = t if from_sink else s
197
        heights = reverse_bfs(src)
198
        if not from_sink:
199
            # s must be reachable from t. Remove t explicitly.
200
            del heights[t]
201
        max_height = max(heights.values())
202
        if from_sink:
203
            # Also mark nodes from which t is unreachable for relabeling. This
204
            # serves the same purpose as the gap heuristic.
205
            for u in R:
206
                if u not in heights and R_nodes[u]['height'] < n:
207
                    heights[u] = n + 1
208
        else:
209
            # Shift the computed heights because the height of s is n.
210
            for u in heights:
211
                heights[u] += n
212
            max_height += n
213
        del heights[src]
214
        for u, new_height in heights.items():
215
            old_height = R_nodes[u]['height']
216
            if new_height != old_height:
217
                if u in levels[old_height].active:
218
                    levels[old_height].active.remove(u)
219
                    levels[new_height].active.add(u)
220
                else:
221
                    levels[old_height].inactive.remove(u)
222
                    levels[new_height].inactive.add(u)
223
                R_nodes[u]['height'] = new_height
224
        return max_height
225

    
226
    # Phase 1: Find the maximum preflow by pushing as much flow as possible to
227
    # t.
228

    
229
    height = max_height
230
    while height > 0:
231
        # Discharge active nodes in the current level.
232
        while True:
233
            level = levels[height]
234
            if not level.active:
235
                # All active nodes in the current level have been discharged.
236
                # Move to the next lower level.
237
                height -= 1
238
                break
239
            # Record the old height and level for the gap heuristic.
240
            old_height = height
241
            old_level = level
242
            u = arbitrary_element(level.active)
243
            height = discharge(u, True)
244
            if grt.is_reached():
245
                # Global relabeling heuristic: Recompute the exact heights of
246
                # all nodes.
247
                height = global_relabel(True)
248
                max_height = height
249
                grt.clear_work()
250
            elif not old_level.active and not old_level.inactive:
251
                # Gap heuristic: If the level at old_height is empty (a 'gap'),
252
                # a minimum cut has been identified. All nodes with heights
253
                # above old_height can have their heights set to n + 1 and not
254
                # be further processed before a maximum preflow is found.
255
                gap_heuristic(old_height)
256
                height = old_height - 1
257
                max_height = height
258
            else:
259
                # Update the height of the highest level with at least one
260
                # active node.
261
                max_height = max(max_height, height)
262

    
263
    # A maximum preflow has been found. The excess at t is the maximum flow
264
    # value.
265
    if value_only:
266
        R.graph['flow_value'] = R_nodes[t]['excess']
267
        return R
268

    
269
    # Phase 2: Convert the maximum preflow into a maximum flow by returning the
270
    # excess to s.
271

    
272
    # Relabel all nodes so that they have accurate heights.
273
    height = global_relabel(False)
274
    grt.clear_work()
275

    
276
    # Continue to discharge the active nodes.
277
    while height > n:
278
        # Discharge active nodes in the current level.
279
        while True:
280
            level = levels[height]
281
            if not level.active:
282
                # All active nodes in the current level have been discharged.
283
                # Move to the next lower level.
284
                height -= 1
285
                break
286
            u = arbitrary_element(level.active)
287
            height = discharge(u, False)
288
            if grt.is_reached():
289
                # Global relabeling heuristic.
290
                height = global_relabel(False)
291
                grt.clear_work()
292

    
293
    R.graph['flow_value'] = R_nodes[t]['excess']
294
    return R
295

    
296

    
297
def preflow_push(G, s, t, capacity='capacity', residual=None,
298
                 global_relabel_freq=1, value_only=False):
299
    r"""Find a maximum single-commodity flow using the highest-label
300
    preflow-push algorithm.
301

302
    This function returns the residual network resulting after computing
303
    the maximum flow. See below for details about the conventions
304
    NetworkX uses for defining residual networks.
305

306
    This algorithm has a running time of $O(n^2 \sqrt{m})$ for $n$ nodes and
307
    $m$ edges.
308

309

310
    Parameters
311
    ----------
312
    G : NetworkX graph
313
        Edges of the graph are expected to have an attribute called
314
        'capacity'. If this attribute is not present, the edge is
315
        considered to have infinite capacity.
316

317
    s : node
318
        Source node for the flow.
319

320
    t : node
321
        Sink node for the flow.
322

323
    capacity : string
324
        Edges of the graph G are expected to have an attribute capacity
325
        that indicates how much flow the edge can support. If this
326
        attribute is not present, the edge is considered to have
327
        infinite capacity. Default value: 'capacity'.
328

329
    residual : NetworkX graph
330
        Residual network on which the algorithm is to be executed. If None, a
331
        new residual network is created. Default value: None.
332

333
    global_relabel_freq : integer, float
334
        Relative frequency of applying the global relabeling heuristic to speed
335
        up the algorithm. If it is None, the heuristic is disabled. Default
336
        value: 1.
337

338
    value_only : bool
339
        If False, compute a maximum flow; otherwise, compute a maximum preflow
340
        which is enough for computing the maximum flow value. Default value:
341
        False.
342

343
    Returns
344
    -------
345
    R : NetworkX DiGraph
346
        Residual network after computing the maximum flow.
347

348
    Raises
349
    ------
350
    NetworkXError
351
        The algorithm does not support MultiGraph and MultiDiGraph. If
352
        the input graph is an instance of one of these two classes, a
353
        NetworkXError is raised.
354

355
    NetworkXUnbounded
356
        If the graph has a path of infinite capacity, the value of a
357
        feasible flow on the graph is unbounded above and the function
358
        raises a NetworkXUnbounded.
359

360
    See also
361
    --------
362
    :meth:`maximum_flow`
363
    :meth:`minimum_cut`
364
    :meth:`edmonds_karp`
365
    :meth:`shortest_augmenting_path`
366

367
    Notes
368
    -----
369
    The residual network :samp:`R` from an input graph :samp:`G` has the
370
    same nodes as :samp:`G`. :samp:`R` is a DiGraph that contains a pair
371
    of edges :samp:`(u, v)` and :samp:`(v, u)` iff :samp:`(u, v)` is not a
372
    self-loop, and at least one of :samp:`(u, v)` and :samp:`(v, u)` exists
373
    in :samp:`G`. For each node :samp:`u` in :samp:`R`,
374
    :samp:`R.nodes[u]['excess']` represents the difference between flow into
375
    :samp:`u` and flow out of :samp:`u`.
376

377
    For each edge :samp:`(u, v)` in :samp:`R`, :samp:`R[u][v]['capacity']`
378
    is equal to the capacity of :samp:`(u, v)` in :samp:`G` if it exists
379
    in :samp:`G` or zero otherwise. If the capacity is infinite,
380
    :samp:`R[u][v]['capacity']` will have a high arbitrary finite value
381
    that does not affect the solution of the problem. This value is stored in
382
    :samp:`R.graph['inf']`. For each edge :samp:`(u, v)` in :samp:`R`,
383
    :samp:`R[u][v]['flow']` represents the flow function of :samp:`(u, v)` and
384
    satisfies :samp:`R[u][v]['flow'] == -R[v][u]['flow']`.
385

386
    The flow value, defined as the total flow into :samp:`t`, the sink, is
387
    stored in :samp:`R.graph['flow_value']`. Reachability to :samp:`t` using
388
    only edges :samp:`(u, v)` such that
389
    :samp:`R[u][v]['flow'] < R[u][v]['capacity']` induces a minimum
390
    :samp:`s`-:samp:`t` cut.
391

392
    Examples
393
    --------
394
    >>> import networkx as nx
395
    >>> from networkx.algorithms.flow import preflow_push
396

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

401
    >>> G = nx.DiGraph()
402
    >>> G.add_edge('x','a', capacity=3.0)
403
    >>> G.add_edge('x','b', capacity=1.0)
404
    >>> G.add_edge('a','c', capacity=3.0)
405
    >>> G.add_edge('b','c', capacity=5.0)
406
    >>> G.add_edge('b','d', capacity=4.0)
407
    >>> G.add_edge('d','e', capacity=2.0)
408
    >>> G.add_edge('c','y', capacity=2.0)
409
    >>> G.add_edge('e','y', capacity=3.0)
410
    >>> R = preflow_push(G, 'x', 'y')
411
    >>> flow_value = nx.maximum_flow_value(G, 'x', 'y')
412
    >>> flow_value == R.graph['flow_value']
413
    True
414
    >>> # preflow_push also stores the maximum flow value
415
    >>> # in the excess attribute of the sink node t
416
    >>> flow_value == R.nodes['y']['excess']
417
    True
418
    >>> # For some problems, you might only want to compute a
419
    >>> # maximum preflow.
420
    >>> R = preflow_push(G, 'x', 'y', value_only=True)
421
    >>> flow_value == R.graph['flow_value']
422
    True
423
    >>> flow_value == R.nodes['y']['excess']
424
    True
425

426
    """
427
    R = preflow_push_impl(G, s, t, capacity, residual, global_relabel_freq,
428
                          value_only)
429
    R.graph['algorithm'] = 'preflow_push'
430
    return R