Statistics
| Branch: | Revision:

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

History | View | Annotate | Download (14.2 KB)

1
# -*- coding: utf-8 -*-
2
"""
3
Capacity scaling minimum cost flow algorithm.
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
__all__ = ['capacity_scaling']
12

    
13
from itertools import chain
14
from math import log
15
import networkx as nx
16
from ...utils import BinaryHeap
17
from ...utils import generate_unique_node
18
from ...utils import not_implemented_for
19
from ...utils import arbitrary_element
20

    
21

    
22
def _detect_unboundedness(R):
23
    """Detect infinite-capacity negative cycles.
24
    """
25
    s = generate_unique_node()
26
    G = nx.DiGraph()
27
    G.add_nodes_from(R)
28

    
29
    # Value simulating infinity.
30
    inf = R.graph['inf']
31
    # True infinity.
32
    f_inf = float('inf')
33
    for u in R:
34
        for v, e in R[u].items():
35
            # Compute the minimum weight of infinite-capacity (u, v) edges.
36
            w = f_inf
37
            for k, e in e.items():
38
                if e['capacity'] == inf:
39
                    w = min(w, e['weight'])
40
            if w != f_inf:
41
                G.add_edge(u, v, weight=w)
42

    
43
    if nx.negative_edge_cycle(G):
44
        raise nx.NetworkXUnbounded(
45
            'Negative cost cycle of infinite capacity found. '
46
            'Min cost flow may be unbounded below.')
47

    
48

    
49
@not_implemented_for('undirected')
50
def _build_residual_network(G, demand, capacity, weight):
51
    """Build a residual network and initialize a zero flow.
52
    """
53
    if sum(G.nodes[u].get(demand, 0) for u in G) != 0:
54
        raise nx.NetworkXUnfeasible("Sum of the demands should be 0.")
55

    
56
    R = nx.MultiDiGraph()
57
    R.add_nodes_from((u, {'excess': -G.nodes[u].get(demand, 0),
58
                          'potential': 0}) for u in G)
59

    
60
    inf = float('inf')
61
    # Detect selfloops with infinite capacities and negative weights.
62
    for u, v, e in nx.selfloop_edges(G, data=True):
63
        if e.get(weight, 0) < 0 and e.get(capacity, inf) == inf:
64
            raise nx.NetworkXUnbounded(
65
                'Negative cost cycle of infinite capacity found. '
66
                'Min cost flow may be unbounded below.')
67

    
68
    # Extract edges with positive capacities. Self loops excluded.
69
    if G.is_multigraph():
70
        edge_list = [(u, v, k, e)
71
                     for u, v, k, e in G.edges(data=True, keys=True)
72
                     if u != v and e.get(capacity, inf) > 0]
73
    else:
74
        edge_list = [(u, v, 0, e) for u, v, e in G.edges(data=True)
75
                     if u != v and e.get(capacity, inf) > 0]
76
    # Simulate infinity with the larger of the sum of absolute node imbalances
77
    # the sum of finite edge capacities or any positive value if both sums are
78
    # zero. This allows the infinite-capacity edges to be distinguished for
79
    # unboundedness detection and directly participate in residual capacity
80
    # calculation.
81
    inf = max(sum(abs(R.nodes[u]['excess']) for u in R),
82
              2 * sum(e[capacity] for u, v, k, e in edge_list
83
                      if capacity in e and e[capacity] != inf)) or 1
84
    for u, v, k, e in edge_list:
85
        r = min(e.get(capacity, inf), inf)
86
        w = e.get(weight, 0)
87
        # Add both (u, v) and (v, u) into the residual network marked with the
88
        # original key. (key[1] == True) indicates the (u, v) is in the
89
        # original network.
90
        R.add_edge(u, v, key=(k, True), capacity=r, weight=w, flow=0)
91
        R.add_edge(v, u, key=(k, False), capacity=0, weight=-w, flow=0)
92

    
93
    # Record the value simulating infinity.
94
    R.graph['inf'] = inf
95

    
96
    _detect_unboundedness(R)
97

    
98
    return R
99

    
100

    
101
def _build_flow_dict(G, R, capacity, weight):
102
    """Build a flow dictionary from a residual network.
103
    """
104
    inf = float('inf')
105
    flow_dict = {}
106
    if G.is_multigraph():
107
        for u in G:
108
            flow_dict[u] = {}
109
            for v, es in G[u].items():
110
                flow_dict[u][v] = dict(
111
                    # Always saturate negative selfloops.
112
                    (k, (0 if (u != v or e.get(capacity, inf) <= 0 or
113
                               e.get(weight, 0) >= 0) else e[capacity]))
114
                    for k, e in es.items())
115
            for v, es in R[u].items():
116
                if v in flow_dict[u]:
117
                    flow_dict[u][v].update((k[0], e['flow'])
118
                                           for k, e in es.items()
119
                                           if e['flow'] > 0)
120
    else:
121
        for u in G:
122
            flow_dict[u] = dict(
123
                # Always saturate negative selfloops.
124
                (v, (0 if (u != v or e.get(capacity, inf) <= 0 or
125
                           e.get(weight, 0) >= 0) else e[capacity]))
126
                for v, e in G[u].items())
127
            flow_dict[u].update((v, e['flow']) for v, es in R[u].items()
128
                                for e in es.values() if e['flow'] > 0)
129
    return flow_dict
130

    
131

    
132
def capacity_scaling(G, demand='demand', capacity='capacity', weight='weight',
133
                     heap=BinaryHeap):
134
    r"""Find a minimum cost flow satisfying all demands in digraph G.
135

136
    This is a capacity scaling successive shortest augmenting path algorithm.
137

138
    G is a digraph with edge costs and capacities and in which nodes
139
    have demand, i.e., they want to send or receive some amount of
140
    flow. A negative demand means that the node wants to send flow, a
141
    positive demand means that the node want to receive flow. A flow on
142
    the digraph G satisfies all demand if the net flow into each node
143
    is equal to the demand of that node.
144

145
    Parameters
146
    ----------
147
    G : NetworkX graph
148
        DiGraph or MultiDiGraph on which a minimum cost flow satisfying all
149
        demands is to be found.
150

151
    demand : string
152
        Nodes of the graph G are expected to have an attribute demand
153
        that indicates how much flow a node wants to send (negative
154
        demand) or receive (positive demand). Note that the sum of the
155
        demands should be 0 otherwise the problem in not feasible. If
156
        this attribute is not present, a node is considered to have 0
157
        demand. Default value: 'demand'.
158

159
    capacity : string
160
        Edges of the graph G are expected to have an attribute capacity
161
        that indicates how much flow the edge can support. If this
162
        attribute is not present, the edge is considered to have
163
        infinite capacity. Default value: 'capacity'.
164

165
    weight : string
166
        Edges of the graph G are expected to have an attribute weight
167
        that indicates the cost incurred by sending one unit of flow on
168
        that edge. If not present, the weight is considered to be 0.
169
        Default value: 'weight'.
170

171
    heap : class
172
        Type of heap to be used in the algorithm. It should be a subclass of
173
        :class:`MinHeap` or implement a compatible interface.
174

175
        If a stock heap implementation is to be used, :class:`BinaryHeap` is
176
        recommended over :class:`PairingHeap` for Python implementations without
177
        optimized attribute accesses (e.g., CPython) despite a slower
178
        asymptotic running time. For Python implementations with optimized
179
        attribute accesses (e.g., PyPy), :class:`PairingHeap` provides better
180
        performance. Default value: :class:`BinaryHeap`.
181

182
    Returns
183
    -------
184
    flowCost : integer
185
        Cost of a minimum cost flow satisfying all demands.
186

187
    flowDict : dictionary
188
        If G is a digraph, a dict-of-dicts keyed by nodes such that
189
        flowDict[u][v] is the flow on edge (u, v).
190
        If G is a MultiDiGraph, a dict-of-dicts-of-dicts keyed by nodes
191
        so that flowDict[u][v][key] is the flow on edge (u, v, key).
192

193
    Raises
194
    ------
195
    NetworkXError
196
        This exception is raised if the input graph is not directed,
197
        not connected.
198

199
    NetworkXUnfeasible
200
        This exception is raised in the following situations:
201

202
            * The sum of the demands is not zero. Then, there is no
203
              flow satisfying all demands.
204
            * There is no flow satisfying all demand.
205

206
    NetworkXUnbounded
207
        This exception is raised if the digraph G has a cycle of
208
        negative cost and infinite capacity. Then, the cost of a flow
209
        satisfying all demands is unbounded below.
210

211
    Notes
212
    -----
213
    This algorithm does not work if edge weights are floating-point numbers.
214

215
    See also
216
    --------
217
    :meth:`network_simplex`
218

219
    Examples
220
    --------
221
    A simple example of a min cost flow problem.
222

223
    >>> import networkx as nx
224
    >>> G = nx.DiGraph()
225
    >>> G.add_node('a', demand = -5)
226
    >>> G.add_node('d', demand = 5)
227
    >>> G.add_edge('a', 'b', weight = 3, capacity = 4)
228
    >>> G.add_edge('a', 'c', weight = 6, capacity = 10)
229
    >>> G.add_edge('b', 'd', weight = 1, capacity = 9)
230
    >>> G.add_edge('c', 'd', weight = 2, capacity = 5)
231
    >>> flowCost, flowDict = nx.capacity_scaling(G)
232
    >>> flowCost
233
    24
234
    >>> flowDict # doctest: +SKIP
235
    {'a': {'c': 1, 'b': 4}, 'c': {'d': 1}, 'b': {'d': 4}, 'd': {}}
236

237
    It is possible to change the name of the attributes used for the
238
    algorithm.
239

240
    >>> G = nx.DiGraph()
241
    >>> G.add_node('p', spam = -4)
242
    >>> G.add_node('q', spam = 2)
243
    >>> G.add_node('a', spam = -2)
244
    >>> G.add_node('d', spam = -1)
245
    >>> G.add_node('t', spam = 2)
246
    >>> G.add_node('w', spam = 3)
247
    >>> G.add_edge('p', 'q', cost = 7, vacancies = 5)
248
    >>> G.add_edge('p', 'a', cost = 1, vacancies = 4)
249
    >>> G.add_edge('q', 'd', cost = 2, vacancies = 3)
250
    >>> G.add_edge('t', 'q', cost = 1, vacancies = 2)
251
    >>> G.add_edge('a', 't', cost = 2, vacancies = 4)
252
    >>> G.add_edge('d', 'w', cost = 3, vacancies = 4)
253
    >>> G.add_edge('t', 'w', cost = 4, vacancies = 1)
254
    >>> flowCost, flowDict = nx.capacity_scaling(G, demand = 'spam',
255
    ...                                          capacity = 'vacancies',
256
    ...                                          weight = 'cost')
257
    >>> flowCost
258
    37
259
    >>> flowDict  # doctest: +SKIP
260
    {'a': {'t': 4}, 'd': {'w': 2}, 'q': {'d': 1}, 'p': {'q': 2, 'a': 2}, 't': {'q': 1, 'w': 1}, 'w': {}}
261
    """
262
    R = _build_residual_network(G, demand, capacity, weight)
263

    
264
    inf = float('inf')
265
    # Account cost of negative selfloops.
266
    flow_cost = sum(
267
        0 if e.get(capacity, inf) <= 0 or e.get(weight, 0) >= 0
268
        else e[capacity] * e[weight]
269
        for u, v, e in nx.selfloop_edges(G, data=True))
270

    
271
    # Determine the maxmimum edge capacity.
272
    wmax = max(chain([-inf],
273
                     (e['capacity'] for u, v, e in R.edges(data=True))))
274
    if wmax == -inf:
275
        # Residual network has no edges.
276
        return flow_cost, _build_flow_dict(G, R, capacity, weight)
277

    
278
    R_nodes = R.nodes
279
    R_succ = R.succ
280

    
281
    delta = 2 ** int(log(wmax, 2))
282
    while delta >= 1:
283
        # Saturate Δ-residual edges with negative reduced costs to achieve
284
        # Δ-optimality.
285
        for u in R:
286
            p_u = R_nodes[u]['potential']
287
            for v, es in R_succ[u].items():
288
                for k, e in es.items():
289
                    flow = e['capacity'] - e['flow']
290
                    if e['weight'] - p_u + R_nodes[v]['potential'] < 0:
291
                        flow = e['capacity'] - e['flow']
292
                        if flow >= delta:
293
                            e['flow'] += flow
294
                            R_succ[v][u][(k[0], not k[1])]['flow'] -= flow
295
                            R_nodes[u]['excess'] -= flow
296
                            R_nodes[v]['excess'] += flow
297
        # Determine the Δ-active nodes.
298
        S = set()
299
        T = set()
300
        S_add = S.add
301
        S_remove = S.remove
302
        T_add = T.add
303
        T_remove = T.remove
304
        for u in R:
305
            excess = R_nodes[u]['excess']
306
            if excess >= delta:
307
                S_add(u)
308
            elif excess <= -delta:
309
                T_add(u)
310
        # Repeatedly augment flow from S to T along shortest paths until
311
        # Δ-feasibility is achieved.
312
        while S and T:
313
            s = arbitrary_element(S)
314
            t = None
315
            # Search for a shortest path in terms of reduce costs from s to
316
            # any t in T in the Δ-residual network.
317
            d = {}
318
            pred = {s: None}
319
            h = heap()
320
            h_insert = h.insert
321
            h_get = h.get
322
            h_insert(s, 0)
323
            while h:
324
                u, d_u = h.pop()
325
                d[u] = d_u
326
                if u in T:
327
                    # Path found.
328
                    t = u
329
                    break
330
                p_u = R_nodes[u]['potential']
331
                for v, es in R_succ[u].items():
332
                    if v in d:
333
                        continue
334
                    wmin = inf
335
                    # Find the minimum-weighted (u, v) Δ-residual edge.
336
                    for k, e in es.items():
337
                        if e['capacity'] - e['flow'] >= delta:
338
                            w = e['weight']
339
                            if w < wmin:
340
                                wmin = w
341
                                kmin = k
342
                                emin = e
343
                    if wmin == inf:
344
                        continue
345
                    # Update the distance label of v.
346
                    d_v = d_u + wmin - p_u + R_nodes[v]['potential']
347
                    if h_insert(v, d_v):
348
                        pred[v] = (u, kmin, emin)
349
            if t is not None:
350
                # Augment Δ units of flow from s to t.
351
                while u != s:
352
                    v = u
353
                    u, k, e = pred[v]
354
                    e['flow'] += delta
355
                    R_succ[v][u][(k[0], not k[1])]['flow'] -= delta
356
                # Account node excess and deficit.
357
                R_nodes[s]['excess'] -= delta
358
                R_nodes[t]['excess'] += delta
359
                if R_nodes[s]['excess'] < delta:
360
                    S_remove(s)
361
                if R_nodes[t]['excess'] > -delta:
362
                    T_remove(t)
363
                # Update node potentials.
364
                d_t = d[t]
365
                for u, d_u in d.items():
366
                    R_nodes[u]['potential'] -= d_u - d_t
367
            else:
368
                # Path not found.
369
                S_remove(s)
370
        delta //= 2
371

    
372
    if any(R.nodes[u]['excess'] != 0 for u in R):
373
        raise nx.NetworkXUnfeasible('No flow satisfying all demands.')
374

    
375
    # Calculate the flow cost.
376
    for u in R:
377
        for v, es in R_succ[u].items():
378
            for e in es.values():
379
                flow = e['flow']
380
                if flow > 0:
381
                    flow_cost += flow * e['weight']
382

    
383
    return flow_cost, _build_flow_dict(G, R, capacity, weight)