Statistics
| Branch: | Revision:

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

History | View | Annotate | Download (20.7 KB)

1
# -*- coding: utf-8 -*-
2
"""
3
Minimum cost flow algorithms on directed connected graphs.
4
"""
5

    
6
__author__ = """Loïc Séguin-C. <loicseguin@gmail.com>"""
7
# Copyright (C) 2010 Loïc Séguin-C. <loicseguin@gmail.com>
8
# All rights reserved.
9
# BSD license.
10

    
11
__all__ = ['network_simplex']
12

    
13
from itertools import chain, islice, repeat
14
from math import ceil, sqrt
15
import networkx as nx
16
from networkx.utils import not_implemented_for
17

    
18
try:
19
    from itertools import izip as zip
20
except ImportError:
21
    pass
22
try:
23
    range = xrange
24
except NameError:
25
    pass
26

    
27

    
28
@not_implemented_for('undirected')
29
def network_simplex(G, demand='demand', capacity='capacity', weight='weight'):
30
    r"""Find a minimum cost flow satisfying all demands in digraph G.
31

32
    This is a primal network simplex algorithm that uses the leaving
33
    arc rule to prevent cycling.
34

35
    G is a digraph with edge costs and capacities and in which nodes
36
    have demand, i.e., they want to send or receive some amount of
37
    flow. A negative demand means that the node wants to send flow, a
38
    positive demand means that the node want to receive flow. A flow on
39
    the digraph G satisfies all demand if the net flow into each node
40
    is equal to the demand of that node.
41

42
    Parameters
43
    ----------
44
    G : NetworkX graph
45
        DiGraph on which a minimum cost flow satisfying all demands is
46
        to be found.
47

48
    demand : string
49
        Nodes of the graph G are expected to have an attribute demand
50
        that indicates how much flow a node wants to send (negative
51
        demand) or receive (positive demand). Note that the sum of the
52
        demands should be 0 otherwise the problem in not feasible. If
53
        this attribute is not present, a node is considered to have 0
54
        demand. Default value: 'demand'.
55

56
    capacity : string
57
        Edges of the graph G are expected to have an attribute capacity
58
        that indicates how much flow the edge can support. If this
59
        attribute is not present, the edge is considered to have
60
        infinite capacity. Default value: 'capacity'.
61

62
    weight : string
63
        Edges of the graph G are expected to have an attribute weight
64
        that indicates the cost incurred by sending one unit of flow on
65
        that edge. If not present, the weight is considered to be 0.
66
        Default value: 'weight'.
67

68
    Returns
69
    -------
70
    flowCost : integer, float
71
        Cost of a minimum cost flow satisfying all demands.
72

73
    flowDict : dictionary
74
        Dictionary of dictionaries keyed by nodes such that
75
        flowDict[u][v] is the flow edge (u, v).
76

77
    Raises
78
    ------
79
    NetworkXError
80
        This exception is raised if the input graph is not directed,
81
        not connected or is a multigraph.
82

83
    NetworkXUnfeasible
84
        This exception is raised in the following situations:
85

86
            * The sum of the demands is not zero. Then, there is no
87
              flow satisfying all demands.
88
            * There is no flow satisfying all demand.
89

90
    NetworkXUnbounded
91
        This exception is raised if the digraph G has a cycle of
92
        negative cost and infinite capacity. Then, the cost of a flow
93
        satisfying all demands is unbounded below.
94

95
    Notes
96
    -----
97
    This algorithm is not guaranteed to work if edge weights or demands
98
    are floating point numbers (overflows and roundoff errors can
99
    cause problems). As a workaround you can use integer numbers by
100
    multiplying the relevant edge attributes by a convenient
101
    constant factor (eg 100).
102

103
    See also
104
    --------
105
    cost_of_flow, max_flow_min_cost, min_cost_flow, min_cost_flow_cost
106

107
    Examples
108
    --------
109
    A simple example of a min cost flow problem.
110

111
    >>> import networkx as nx
112
    >>> G = nx.DiGraph()
113
    >>> G.add_node('a', demand=-5)
114
    >>> G.add_node('d', demand=5)
115
    >>> G.add_edge('a', 'b', weight=3, capacity=4)
116
    >>> G.add_edge('a', 'c', weight=6, capacity=10)
117
    >>> G.add_edge('b', 'd', weight=1, capacity=9)
118
    >>> G.add_edge('c', 'd', weight=2, capacity=5)
119
    >>> flowCost, flowDict = nx.network_simplex(G)
120
    >>> flowCost
121
    24
122
    >>> flowDict # doctest: +SKIP
123
    {'a': {'c': 1, 'b': 4}, 'c': {'d': 1}, 'b': {'d': 4}, 'd': {}}
124

125
    The mincost flow algorithm can also be used to solve shortest path
126
    problems. To find the shortest path between two nodes u and v,
127
    give all edges an infinite capacity, give node u a demand of -1 and
128
    node v a demand a 1. Then run the network simplex. The value of a
129
    min cost flow will be the distance between u and v and edges
130
    carrying positive flow will indicate the path.
131

132
    >>> G=nx.DiGraph()
133
    >>> G.add_weighted_edges_from([('s', 'u' ,10), ('s' ,'x' ,5),
134
    ...                            ('u', 'v' ,1), ('u' ,'x' ,2),
135
    ...                            ('v', 'y' ,1), ('x' ,'u' ,3),
136
    ...                            ('x', 'v' ,5), ('x' ,'y' ,2),
137
    ...                            ('y', 's' ,7), ('y' ,'v' ,6)])
138
    >>> G.add_node('s', demand = -1)
139
    >>> G.add_node('v', demand = 1)
140
    >>> flowCost, flowDict = nx.network_simplex(G)
141
    >>> flowCost == nx.shortest_path_length(G, 's', 'v', weight='weight')
142
    True
143
    >>> sorted([(u, v) for u in flowDict for v in flowDict[u] if flowDict[u][v] > 0])
144
    [('s', 'x'), ('u', 'v'), ('x', 'u')]
145
    >>> nx.shortest_path(G, 's', 'v', weight = 'weight')
146
    ['s', 'x', 'u', 'v']
147

148
    It is possible to change the name of the attributes used for the
149
    algorithm.
150

151
    >>> G = nx.DiGraph()
152
    >>> G.add_node('p', spam=-4)
153
    >>> G.add_node('q', spam=2)
154
    >>> G.add_node('a', spam=-2)
155
    >>> G.add_node('d', spam=-1)
156
    >>> G.add_node('t', spam=2)
157
    >>> G.add_node('w', spam=3)
158
    >>> G.add_edge('p', 'q', cost=7, vacancies=5)
159
    >>> G.add_edge('p', 'a', cost=1, vacancies=4)
160
    >>> G.add_edge('q', 'd', cost=2, vacancies=3)
161
    >>> G.add_edge('t', 'q', cost=1, vacancies=2)
162
    >>> G.add_edge('a', 't', cost=2, vacancies=4)
163
    >>> G.add_edge('d', 'w', cost=3, vacancies=4)
164
    >>> G.add_edge('t', 'w', cost=4, vacancies=1)
165
    >>> flowCost, flowDict = nx.network_simplex(G, demand='spam',
166
    ...                                         capacity='vacancies',
167
    ...                                         weight='cost')
168
    >>> flowCost
169
    37
170
    >>> flowDict  # doctest: +SKIP
171
    {'a': {'t': 4}, 'd': {'w': 2}, 'q': {'d': 1}, 'p': {'q': 2, 'a': 2}, 't': {'q': 1, 'w': 1}, 'w': {}}
172

173
    References
174
    ----------
175
    .. [1] Z. Kiraly, P. Kovacs.
176
           Efficient implementation of minimum-cost flow algorithms.
177
           Acta Universitatis Sapientiae, Informatica 4(1):67--118. 2012.
178
    .. [2] R. Barr, F. Glover, D. Klingman.
179
           Enhancement of spanning tree labeling procedures for network
180
           optimization.
181
           INFOR 17(1):16--34. 1979.
182
    """
183
    ###########################################################################
184
    # Problem essentials extraction and sanity check
185
    ###########################################################################
186

    
187
    if len(G) == 0:
188
        raise nx.NetworkXError('graph has no nodes')
189

    
190
    # Number all nodes and edges and hereafter reference them using ONLY their
191
    # numbers
192

    
193
    N = list(G)                                # nodes
194
    I = {u: i for i, u in enumerate(N)}        # node indices
195
    D = [G.nodes[u].get(demand, 0) for u in N]  # node demands
196

    
197
    inf = float('inf')
198
    for p, b in zip(N, D):
199
        if abs(b) == inf:
200
            raise nx.NetworkXError('node %r has infinite demand' % (p,))
201

    
202
    multigraph = G.is_multigraph()
203
    S = []  # edge sources
204
    T = []  # edge targets
205
    if multigraph:
206
        K = []  # edge keys
207
    E = {}  # edge indices
208
    U = []  # edge capacities
209
    C = []  # edge weights
210

    
211
    if not multigraph:
212
        edges = G.edges(data=True)
213
    else:
214
        edges = G.edges(data=True, keys=True)
215
    edges = (e for e in edges
216
             if e[0] != e[1] and e[-1].get(capacity, inf) != 0)
217
    for i, e in enumerate(edges):
218
        S.append(I[e[0]])
219
        T.append(I[e[1]])
220
        if multigraph:
221
            K.append(e[2])
222
        E[e[:-1]] = i
223
        U.append(e[-1].get(capacity, inf))
224
        C.append(e[-1].get(weight, 0))
225

    
226
    for e, c in zip(E, C):
227
        if abs(c) == inf:
228
            raise nx.NetworkXError('edge %r has infinite weight' % (e,))
229
    if not multigraph:
230
        edges = nx.selfloop_edges(G, data=True)
231
    else:
232
        edges = nx.selfloop_edges(G, data=True, keys=True)
233
    for e in edges:
234
        if abs(e[-1].get(weight, 0)) == inf:
235
            raise nx.NetworkXError('edge %r has infinite weight' % (e[:-1],))
236

    
237
    ###########################################################################
238
    # Quick infeasibility detection
239
    ###########################################################################
240

    
241
    if sum(D) != 0:
242
        raise nx.NetworkXUnfeasible('total node demand is not zero')
243
    for e, u in zip(E, U):
244
        if u < 0:
245
            raise nx.NetworkXUnfeasible('edge %r has negative capacity' % (e,))
246
    if not multigraph:
247
        edges = nx.selfloop_edges(G, data=True)
248
    else:
249
        edges = nx.selfloop_edges(G, data=True, keys=True)
250
    for e in edges:
251
        if e[-1].get(capacity, inf) < 0:
252
            raise nx.NetworkXUnfeasible(
253
                'edge %r has negative capacity' % (e[:-1],))
254

    
255
    ###########################################################################
256
    # Initialization
257
    ###########################################################################
258

    
259
    # Add a dummy node -1 and connect all existing nodes to it with infinite-
260
    # capacity dummy edges. Node -1 will serve as the root of the
261
    # spanning tree of the network simplex method. The new edges will used to
262
    # trivially satisfy the node demands and create an initial strongly
263
    # feasible spanning tree.
264
    n = len(N)  # number of nodes
265
    for p, d in enumerate(D):
266
        if d > 0:  # Must be greater-than here. Zero-demand nodes must have
267
                   # edges pointing towards the root to ensure strong
268
                   # feasibility.
269
            S.append(-1)
270
            T.append(p)
271
        else:
272
            S.append(p)
273
            T.append(-1)
274
    faux_inf = 3 * max(chain([sum(u for u in U if u < inf),
275
                              sum(abs(c) for c in C)],
276
                             (abs(d) for d in D))) or 1
277
    C.extend(repeat(faux_inf, n))
278
    U.extend(repeat(faux_inf, n))
279

    
280
    # Construct the initial spanning tree.
281
    e = len(E)                                           # number of edges
282
    x = list(chain(repeat(0, e), (abs(d) for d in D)))   # edge flows
283
    pi = [faux_inf if d <= 0 else -faux_inf for d in D]  # node potentials
284
    parent = list(chain(repeat(-1, n), [None]))  # parent nodes
285
    edge = list(range(e, e + n))                 # edges to parents
286
    size = list(chain(repeat(1, n), [n + 1]))    # subtree sizes
287
    next = list(chain(range(1, n), [-1, 0]))     # next nodes in depth-first thread
288
    prev = list(range(-1, n))                    # previous nodes in depth-first thread
289
    last = list(chain(range(n), [n - 1]))        # last descendants in depth-first thread
290

    
291
    ###########################################################################
292
    # Pivot loop
293
    ###########################################################################
294

    
295
    def reduced_cost(i):
296
        """Returns the reduced cost of an edge i.
297
        """
298
        c = C[i] - pi[S[i]] + pi[T[i]]
299
        return c if x[i] == 0 else -c
300

    
301
    def find_entering_edges():
302
        """Yield entering edges until none can be found.
303
        """
304
        if e == 0:
305
            return
306

    
307
        # Entering edges are found by combining Dantzig's rule and Bland's
308
        # rule. The edges are cyclically grouped into blocks of size B. Within
309
        # each block, Dantzig's rule is applied to find an entering edge. The
310
        # blocks to search is determined following Bland's rule.
311
        B = int(ceil(sqrt(e)))  # pivot block size
312
        M = (e + B - 1) // B    # number of blocks needed to cover all edges
313
        m = 0                   # number of consecutive blocks without eligible
314
        # entering edges
315
        f = 0                   # first edge in block
316
        while m < M:
317
            # Determine the next block of edges.
318
            l = f + B
319
            if l <= e:
320
                edges = range(f, l)
321
            else:
322
                l -= e
323
                edges = chain(range(f, e), range(l))
324
            f = l
325
            # Find the first edge with the lowest reduced cost.
326
            i = min(edges, key=reduced_cost)
327
            c = reduced_cost(i)
328
            if c >= 0:
329
                # No entering edge found in the current block.
330
                m += 1
331
            else:
332
                # Entering edge found.
333
                if x[i] == 0:
334
                    p = S[i]
335
                    q = T[i]
336
                else:
337
                    p = T[i]
338
                    q = S[i]
339
                yield i, p, q
340
                m = 0
341
        # All edges have nonnegative reduced costs. The current flow is
342
        # optimal.
343

    
344
    def find_apex(p, q):
345
        """Find the lowest common ancestor of nodes p and q in the spanning
346
        tree.
347
        """
348
        size_p = size[p]
349
        size_q = size[q]
350
        while True:
351
            while size_p < size_q:
352
                p = parent[p]
353
                size_p = size[p]
354
            while size_p > size_q:
355
                q = parent[q]
356
                size_q = size[q]
357
            if size_p == size_q:
358
                if p != q:
359
                    p = parent[p]
360
                    size_p = size[p]
361
                    q = parent[q]
362
                    size_q = size[q]
363
                else:
364
                    return p
365

    
366
    def trace_path(p, w):
367
        """Returns the nodes and edges on the path from node p to its ancestor
368
        w.
369
        """
370
        Wn = [p]
371
        We = []
372
        while p != w:
373
            We.append(edge[p])
374
            p = parent[p]
375
            Wn.append(p)
376
        return Wn, We
377

    
378
    def find_cycle(i, p, q):
379
        """Returns the nodes and edges on the cycle containing edge i == (p, q)
380
        when the latter is added to the spanning tree.
381

382
        The cycle is oriented in the direction from p to q.
383
        """
384
        w = find_apex(p, q)
385
        Wn, We = trace_path(p, w)
386
        Wn.reverse()
387
        We.reverse()
388
        if We != [i]:
389
            We.append(i)
390
        WnR, WeR = trace_path(q, w)
391
        del WnR[-1]
392
        Wn += WnR
393
        We += WeR
394
        return Wn, We
395

    
396
    def residual_capacity(i, p):
397
        """Returns the residual capacity of an edge i in the direction away
398
        from its endpoint p.
399
        """
400
        return U[i] - x[i] if S[i] == p else x[i]
401

    
402
    def find_leaving_edge(Wn, We):
403
        """Returns the leaving edge in a cycle represented by Wn and We.
404
        """
405
        j, s = min(zip(reversed(We), reversed(Wn)),
406
                   key=lambda i_p: residual_capacity(*i_p))
407
        t = T[j] if S[j] == s else S[j]
408
        return j, s, t
409

    
410
    def augment_flow(Wn, We, f):
411
        """Augment f units of flow along a cycle represented by Wn and We.
412
        """
413
        for i, p in zip(We, Wn):
414
            if S[i] == p:
415
                x[i] += f
416
            else:
417
                x[i] -= f
418

    
419
    def trace_subtree(p):
420
        """Yield the nodes in the subtree rooted at a node p.
421
        """
422
        yield p
423
        l = last[p]
424
        while p != l:
425
            p = next[p]
426
            yield p
427

    
428
    def remove_edge(s, t):
429
        """Remove an edge (s, t) where parent[t] == s from the spanning tree.
430
        """
431
        size_t = size[t]
432
        prev_t = prev[t]
433
        last_t = last[t]
434
        next_last_t = next[last_t]
435
        # Remove (s, t).
436
        parent[t] = None
437
        edge[t] = None
438
        # Remove the subtree rooted at t from the depth-first thread.
439
        next[prev_t] = next_last_t
440
        prev[next_last_t] = prev_t
441
        next[last_t] = t
442
        prev[t] = last_t
443
        # Update the subtree sizes and last descendants of the (old) acenstors
444
        # of t.
445
        while s is not None:
446
            size[s] -= size_t
447
            if last[s] == last_t:
448
                last[s] = prev_t
449
            s = parent[s]
450

    
451
    def make_root(q):
452
        """Make a node q the root of its containing subtree.
453
        """
454
        ancestors = []
455
        while q is not None:
456
            ancestors.append(q)
457
            q = parent[q]
458
        ancestors.reverse()
459
        for p, q in zip(ancestors, islice(ancestors, 1, None)):
460
            size_p = size[p]
461
            last_p = last[p]
462
            prev_q = prev[q]
463
            last_q = last[q]
464
            next_last_q = next[last_q]
465
            # Make p a child of q.
466
            parent[p] = q
467
            parent[q] = None
468
            edge[p] = edge[q]
469
            edge[q] = None
470
            size[p] = size_p - size[q]
471
            size[q] = size_p
472
            # Remove the subtree rooted at q from the depth-first thread.
473
            next[prev_q] = next_last_q
474
            prev[next_last_q] = prev_q
475
            next[last_q] = q
476
            prev[q] = last_q
477
            if last_p == last_q:
478
                last[p] = prev_q
479
                last_p = prev_q
480
            # Add the remaining parts of the subtree rooted at p as a subtree
481
            # of q in the depth-first thread.
482
            prev[p] = last_q
483
            next[last_q] = p
484
            next[last_p] = q
485
            prev[q] = last_p
486
            last[q] = last_p
487

    
488
    def add_edge(i, p, q):
489
        """Add an edge (p, q) to the spanning tree where q is the root of a
490
        subtree.
491
        """
492
        last_p = last[p]
493
        next_last_p = next[last_p]
494
        size_q = size[q]
495
        last_q = last[q]
496
        # Make q a child of p.
497
        parent[q] = p
498
        edge[q] = i
499
        # Insert the subtree rooted at q into the depth-first thread.
500
        next[last_p] = q
501
        prev[q] = last_p
502
        prev[next_last_p] = last_q
503
        next[last_q] = next_last_p
504
        # Update the subtree sizes and last descendants of the (new) ancestors
505
        # of q.
506
        while p is not None:
507
            size[p] += size_q
508
            if last[p] == last_p:
509
                last[p] = last_q
510
            p = parent[p]
511

    
512
    def update_potentials(i, p, q):
513
        """Update the potentials of the nodes in the subtree rooted at a node
514
        q connected to its parent p by an edge i.
515
        """
516
        if q == T[i]:
517
            d = pi[p] - C[i] - pi[q]
518
        else:
519
            d = pi[p] + C[i] - pi[q]
520
        for q in trace_subtree(q):
521
            pi[q] += d
522

    
523
    # Pivot loop
524
    for i, p, q in find_entering_edges():
525
        Wn, We = find_cycle(i, p, q)
526
        j, s, t = find_leaving_edge(Wn, We)
527
        augment_flow(Wn, We, residual_capacity(j, s))
528
        if i != j:  # Do nothing more if the entering edge is the same as the
529
                    # the leaving edge.
530
            if parent[t] != s:
531
                # Ensure that s is the parent of t.
532
                s, t = t, s
533
            if We.index(i) > We.index(j):
534
                # Ensure that q is in the subtree rooted at t.
535
                p, q = q, p
536
            remove_edge(s, t)
537
            make_root(q)
538
            add_edge(i, p, q)
539
            update_potentials(i, p, q)
540

    
541
    ###########################################################################
542
    # Infeasibility and unboundedness detection
543
    ###########################################################################
544

    
545
    if any(x[i] != 0 for i in range(-n, 0)):
546
        raise nx.NetworkXUnfeasible('no flow satisfies all node demands')
547

    
548
    if (any(x[i] * 2 >= faux_inf for i in range(e)) or
549
        any(e[-1].get(capacity, inf) == inf and e[-1].get(weight, 0) < 0
550
            for e in nx.selfloop_edges(G, data=True))):
551
        raise nx.NetworkXUnbounded(
552
            'negative cycle with infinite capacity found')
553

    
554
    ###########################################################################
555
    # Flow cost calculation and flow dict construction
556
    ###########################################################################
557

    
558
    del x[e:]
559
    flow_cost = sum(c * x for c, x in zip(C, x))
560
    flow_dict = {n: {} for n in N}
561

    
562
    def add_entry(e):
563
        """Add a flow dict entry.
564
        """
565
        d = flow_dict[e[0]]
566
        for k in e[1:-2]:
567
            try:
568
                d = d[k]
569
            except KeyError:
570
                t = {}
571
                d[k] = t
572
                d = t
573
        d[e[-2]] = e[-1]
574

    
575
    S = (N[s] for s in S)  # Use original nodes.
576
    T = (N[t] for t in T)  # Use original nodes.
577
    if not multigraph:
578
        for e in zip(S, T, x):
579
            add_entry(e)
580
        edges = G.edges(data=True)
581
    else:
582
        for e in zip(S, T, K, x):
583
            add_entry(e)
584
        edges = G.edges(data=True, keys=True)
585
    for e in edges:
586
        if e[0] != e[1]:
587
            if e[-1].get(capacity, inf) == 0:
588
                add_entry(e[:-1] + (0,))
589
        else:
590
            c = e[-1].get(weight, 0)
591
            if c >= 0:
592
                add_entry(e[:-1] + (0,))
593
            else:
594
                u = e[-1][capacity]
595
                flow_cost += c * u
596
                add_entry(e[:-1] + (u,))
597

    
598
    return flow_cost, flow_dict