COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/capacity_scaling.h

Last change on this file was 2629:84354c78b068, checked in by Peter Kovacs, 15 years ago

Improved constructors for min cost flow classes
Removing the non-zero lower bounds is faster

File size: 22.2 KB
RevLine 
[2440]1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
[2553]5 * Copyright (C) 2003-2008
[2440]6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
12 *
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_CAPACITY_SCALING_H
20#define LEMON_CAPACITY_SCALING_H
21
22/// \ingroup min_cost_flow
23///
24/// \file
[2574]25/// \brief Capacity scaling algorithm for finding a minimum cost flow.
26
27#include <vector>
[2535]28#include <lemon/bin_heap.h>
[2457]29
[2440]30namespace lemon {
31
32  /// \addtogroup min_cost_flow
33  /// @{
34
[2574]35  /// \brief Implementation of the capacity scaling algorithm for
36  /// finding a minimum cost flow.
[2440]37  ///
[2535]38  /// \ref CapacityScaling implements the capacity scaling version
39  /// of the successive shortest path algorithm for finding a minimum
40  /// cost flow.
[2440]41  ///
[2574]42  /// \tparam Graph The directed graph type the algorithm runs on.
43  /// \tparam LowerMap The type of the lower bound map.
44  /// \tparam CapacityMap The type of the capacity (upper bound) map.
45  /// \tparam CostMap The type of the cost (length) map.
46  /// \tparam SupplyMap The type of the supply map.
[2440]47  ///
48  /// \warning
[2574]49  /// - Edge capacities and costs should be \e non-negative \e integers.
50  /// - Supply values should be \e signed \e integers.
[2581]51  /// - The value types of the maps should be convertible to each other.
52  /// - \c CostMap::Value must be signed type.
[2440]53  ///
54  /// \author Peter Kovacs
[2533]55  template < typename Graph,
[2535]56             typename LowerMap = typename Graph::template EdgeMap<int>,
[2574]57             typename CapacityMap = typename Graph::template EdgeMap<int>,
[2535]58             typename CostMap = typename Graph::template EdgeMap<int>,
[2574]59             typename SupplyMap = typename Graph::template NodeMap<int> >
[2440]60  class CapacityScaling
61  {
[2556]62    GRAPH_TYPEDEFS(typename Graph);
[2440]63
64    typedef typename CapacityMap::Value Capacity;
65    typedef typename CostMap::Value Cost;
66    typedef typename SupplyMap::Value Supply;
[2556]67    typedef typename Graph::template EdgeMap<Capacity> CapacityEdgeMap;
68    typedef typename Graph::template NodeMap<Supply> SupplyNodeMap;
[2535]69    typedef typename Graph::template NodeMap<Edge> PredMap;
[2440]70
71  public:
72
[2556]73    /// The type of the flow map.
74    typedef typename Graph::template EdgeMap<Capacity> FlowMap;
75    /// The type of the potential map.
[2440]76    typedef typename Graph::template NodeMap<Cost> PotentialMap;
77
[2574]78  private:
[2440]79
[2535]80    /// \brief Special implementation of the \ref Dijkstra algorithm
[2574]81    /// for finding shortest paths in the residual network.
82    ///
83    /// \ref ResidualDijkstra is a special implementation of the
84    /// \ref Dijkstra algorithm for finding shortest paths in the
85    /// residual network of the graph with respect to the reduced edge
86    /// costs and modifying the node potentials according to the
87    /// distance of the nodes.
[2535]88    class ResidualDijkstra
[2440]89    {
[2535]90      typedef typename Graph::template NodeMap<int> HeapCrossRef;
91      typedef BinHeap<Cost, HeapCrossRef> Heap;
92
[2574]93    private:
[2535]94
[2574]95      // The directed graph the algorithm runs on
96      const Graph &_graph;
[2535]97
[2574]98      // The main maps
99      const FlowMap &_flow;
100      const CapacityEdgeMap &_res_cap;
101      const CostMap &_cost;
102      const SupplyNodeMap &_excess;
103      PotentialMap &_potential;
[2535]104
[2574]105      // The distance map
[2588]106      PotentialMap _dist;
[2574]107      // The pred edge map
108      PredMap &_pred;
109      // The processed (i.e. permanently labeled) nodes
110      std::vector<Node> _proc_nodes;
[2440]111
112    public:
113
[2581]114      /// Constructor.
[2574]115      ResidualDijkstra( const Graph &graph,
116                        const FlowMap &flow,
117                        const CapacityEdgeMap &res_cap,
118                        const CostMap &cost,
119                        const SupplyMap &excess,
120                        PotentialMap &potential,
121                        PredMap &pred ) :
122        _graph(graph), _flow(flow), _res_cap(res_cap), _cost(cost),
123        _excess(excess), _potential(potential), _dist(graph),
124        _pred(pred)
[2535]125      {}
[2440]126
[2620]127      /// Run the algorithm from the given source node.
[2588]128      Node run(Node s, Capacity delta = 1) {
[2574]129        HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
[2535]130        Heap heap(heap_cross_ref);
131        heap.push(s, 0);
[2574]132        _pred[s] = INVALID;
133        _proc_nodes.clear();
[2535]134
135        // Processing nodes
[2574]136        while (!heap.empty() && _excess[heap.top()] > -delta) {
[2535]137          Node u = heap.top(), v;
[2574]138          Cost d = heap.prio() + _potential[u], nd;
139          _dist[u] = heap.prio();
[2535]140          heap.pop();
[2574]141          _proc_nodes.push_back(u);
[2535]142
143          // Traversing outgoing edges
[2574]144          for (OutEdgeIt e(_graph, u); e != INVALID; ++e) {
145            if (_res_cap[e] >= delta) {
146              v = _graph.target(e);
[2535]147              switch(heap.state(v)) {
148              case Heap::PRE_HEAP:
[2574]149                heap.push(v, d + _cost[e] - _potential[v]);
150                _pred[v] = e;
[2535]151                break;
152              case Heap::IN_HEAP:
[2574]153                nd = d + _cost[e] - _potential[v];
[2535]154                if (nd < heap[v]) {
155                  heap.decrease(v, nd);
[2574]156                  _pred[v] = e;
[2535]157                }
158                break;
159              case Heap::POST_HEAP:
160                break;
161              }
162            }
163          }
164
165          // Traversing incoming edges
[2574]166          for (InEdgeIt e(_graph, u); e != INVALID; ++e) {
167            if (_flow[e] >= delta) {
168              v = _graph.source(e);
[2535]169              switch(heap.state(v)) {
170              case Heap::PRE_HEAP:
[2574]171                heap.push(v, d - _cost[e] - _potential[v]);
172                _pred[v] = e;
[2535]173                break;
174              case Heap::IN_HEAP:
[2574]175                nd = d - _cost[e] - _potential[v];
[2535]176                if (nd < heap[v]) {
177                  heap.decrease(v, nd);
[2574]178                  _pred[v] = e;
[2535]179                }
180                break;
181              case Heap::POST_HEAP:
182                break;
183              }
184            }
185          }
186        }
187        if (heap.empty()) return INVALID;
188
189        // Updating potentials of processed nodes
190        Node t = heap.top();
[2574]191        Cost t_dist = heap.prio();
192        for (int i = 0; i < int(_proc_nodes.size()); ++i)
193          _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
[2535]194
195        return t;
[2440]196      }
197
[2535]198    }; //class ResidualDijkstra
[2440]199
[2574]200  private:
[2440]201
[2574]202    // The directed graph the algorithm runs on
203    const Graph &_graph;
204    // The original lower bound map
205    const LowerMap *_lower;
206    // The modified capacity map
207    CapacityEdgeMap _capacity;
208    // The original cost map
209    const CostMap &_cost;
210    // The modified supply map
211    SupplyNodeMap _supply;
212    bool _valid_supply;
[2440]213
[2574]214    // Edge map of the current flow
[2581]215    FlowMap *_flow;
216    bool _local_flow;
[2574]217    // Node map of the current potentials
[2581]218    PotentialMap *_potential;
219    bool _local_potential;
[2440]220
[2574]221    // The residual capacity map
222    CapacityEdgeMap _res_cap;
223    // The excess map
224    SupplyNodeMap _excess;
225    // The excess nodes (i.e. nodes with positive excess)
226    std::vector<Node> _excess_nodes;
227    // The deficit nodes (i.e. nodes with negative excess)
228    std::vector<Node> _deficit_nodes;
[2440]229
[2574]230    // The delta parameter used for capacity scaling
231    Capacity _delta;
232    // The maximum number of phases
233    int _phase_num;
[2440]234
[2574]235    // The pred edge map
236    PredMap _pred;
237    // Implementation of the Dijkstra algorithm for finding augmenting
238    // shortest paths in the residual network
[2581]239    ResidualDijkstra *_dijkstra;
[2440]240
[2581]241  public:
[2440]242
[2581]243    /// \brief General constructor (with lower bounds).
[2440]244    ///
[2581]245    /// General constructor (with lower bounds).
[2440]246    ///
[2574]247    /// \param graph The directed graph the algorithm runs on.
248    /// \param lower The lower bounds of the edges.
249    /// \param capacity The capacities (upper bounds) of the edges.
250    /// \param cost The cost (length) values of the edges.
251    /// \param supply The supply values of the nodes (signed).
252    CapacityScaling( const Graph &graph,
253                     const LowerMap &lower,
254                     const CapacityMap &capacity,
255                     const CostMap &cost,
256                     const SupplyMap &supply ) :
[2629]257      _graph(graph), _lower(&lower), _capacity(capacity), _cost(cost),
258      _supply(supply), _flow(NULL), _local_flow(false),
[2623]259      _potential(NULL), _local_potential(false),
[2629]260      _res_cap(capacity), _excess(supply), _pred(graph), _dijkstra(NULL)
[2440]261    {
[2629]262      // Check the sum of supply values
[2440]263      Supply sum = 0;
[2629]264      for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
265      _valid_supply = sum == 0;
266
267      // Remove non-zero lower bounds
268      typename LowerMap::Value lcap;
269      for (EdgeIt e(_graph); e != INVALID; ++e) {
270        if ((lcap = lower[e]) != 0) {
271          _capacity[e] -= lcap;
272          _res_cap[e] -= lcap;
273          _supply[_graph.source(e)] -= lcap;
274          _supply[_graph.target(e)] += lcap;
275          _excess[_graph.source(e)] -= lcap;
276          _excess[_graph.target(e)] += lcap;
277        }
[2440]278      }
279    }
280
[2581]281    /// \brief General constructor (without lower bounds).
[2440]282    ///
[2581]283    /// General constructor (without lower bounds).
[2440]284    ///
[2574]285    /// \param graph The directed graph the algorithm runs on.
286    /// \param capacity The capacities (upper bounds) of the edges.
287    /// \param cost The cost (length) values of the edges.
288    /// \param supply The supply values of the nodes (signed).
289    CapacityScaling( const Graph &graph,
290                     const CapacityMap &capacity,
291                     const CostMap &cost,
292                     const SupplyMap &supply ) :
293      _graph(graph), _lower(NULL), _capacity(capacity), _cost(cost),
[2623]294      _supply(supply), _flow(NULL), _local_flow(false),
295      _potential(NULL), _local_potential(false),
[2629]296      _res_cap(capacity), _excess(supply), _pred(graph), _dijkstra(NULL)
[2440]297    {
[2629]298      // Check the sum of supply values
[2440]299      Supply sum = 0;
[2574]300      for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
301      _valid_supply = sum == 0;
[2440]302    }
303
[2581]304    /// \brief Simple constructor (with lower bounds).
[2440]305    ///
[2581]306    /// Simple constructor (with lower bounds).
[2440]307    ///
[2574]308    /// \param graph The directed graph the algorithm runs on.
309    /// \param lower The lower bounds of the edges.
310    /// \param capacity The capacities (upper bounds) of the edges.
311    /// \param cost The cost (length) values of the edges.
312    /// \param s The source node.
313    /// \param t The target node.
314    /// \param flow_value The required amount of flow from node \c s
315    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
316    CapacityScaling( const Graph &graph,
317                     const LowerMap &lower,
318                     const CapacityMap &capacity,
319                     const CostMap &cost,
320                     Node s, Node t,
321                     Supply flow_value ) :
[2629]322      _graph(graph), _lower(&lower), _capacity(capacity), _cost(cost),
323      _supply(graph, 0), _flow(NULL), _local_flow(false),
[2623]324      _potential(NULL), _local_potential(false),
[2629]325      _res_cap(capacity), _excess(graph, 0), _pred(graph), _dijkstra(NULL)
[2440]326    {
[2629]327      // Remove non-zero lower bounds
328      _supply[s] = _excess[s] =  flow_value;
329      _supply[t] = _excess[t] = -flow_value;
330      typename LowerMap::Value lcap;
331      for (EdgeIt e(_graph); e != INVALID; ++e) {
332        if ((lcap = lower[e]) != 0) {
333          _capacity[e] -= lcap;
334          _res_cap[e] -= lcap;
335          _supply[_graph.source(e)] -= lcap;
336          _supply[_graph.target(e)] += lcap;
337          _excess[_graph.source(e)] -= lcap;
338          _excess[_graph.target(e)] += lcap;
339        }
[2440]340      }
[2574]341      _valid_supply = true;
[2440]342    }
343
[2581]344    /// \brief Simple constructor (without lower bounds).
[2440]345    ///
[2581]346    /// Simple constructor (without lower bounds).
[2440]347    ///
[2574]348    /// \param graph The directed graph the algorithm runs on.
349    /// \param capacity The capacities (upper bounds) of the edges.
350    /// \param cost The cost (length) values of the edges.
351    /// \param s The source node.
352    /// \param t The target node.
353    /// \param flow_value The required amount of flow from node \c s
354    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
355    CapacityScaling( const Graph &graph,
356                     const CapacityMap &capacity,
357                     const CostMap &cost,
358                     Node s, Node t,
359                     Supply flow_value ) :
360      _graph(graph), _lower(NULL), _capacity(capacity), _cost(cost),
[2623]361      _supply(graph, 0), _flow(NULL), _local_flow(false),
362      _potential(NULL), _local_potential(false),
[2629]363      _res_cap(capacity), _excess(graph, 0), _pred(graph), _dijkstra(NULL)
[2440]364    {
[2629]365      _supply[s] = _excess[s] =  flow_value;
366      _supply[t] = _excess[t] = -flow_value;
[2574]367      _valid_supply = true;
[2440]368    }
369
[2581]370    /// Destructor.
371    ~CapacityScaling() {
372      if (_local_flow) delete _flow;
373      if (_local_potential) delete _potential;
374      delete _dijkstra;
375    }
376
[2620]377    /// \brief Set the flow map.
[2581]378    ///
[2620]379    /// Set the flow map.
[2581]380    ///
381    /// \return \c (*this)
382    CapacityScaling& flowMap(FlowMap &map) {
383      if (_local_flow) {
384        delete _flow;
385        _local_flow = false;
386      }
387      _flow = &map;
388      return *this;
389    }
390
[2620]391    /// \brief Set the potential map.
[2581]392    ///
[2620]393    /// Set the potential map.
[2581]394    ///
395    /// \return \c (*this)
396    CapacityScaling& potentialMap(PotentialMap &map) {
397      if (_local_potential) {
398        delete _potential;
399        _local_potential = false;
400      }
401      _potential = &map;
402      return *this;
403    }
404
405    /// \name Execution control
406
407    /// @{
408
[2620]409    /// \brief Run the algorithm.
[2556]410    ///
[2620]411    /// This function runs the algorithm.
[2556]412    ///
[2574]413    /// \param scaling Enable or disable capacity scaling.
[2556]414    /// If the maximum edge capacity and/or the amount of total supply
[2574]415    /// is rather small, the algorithm could be slightly faster without
[2556]416    /// scaling.
417    ///
418    /// \return \c true if a feasible flow can be found.
[2574]419    bool run(bool scaling = true) {
420      return init(scaling) && start();
[2556]421    }
422
[2581]423    /// @}
424
425    /// \name Query Functions
[2620]426    /// The results of the algorithm can be obtained using these
427    /// functions.\n
428    /// \ref lemon::CapacityScaling::run() "run()" must be called before
429    /// using them.
[2581]430
431    /// @{
432
[2620]433    /// \brief Return a const reference to the edge map storing the
[2574]434    /// found flow.
[2440]435    ///
[2620]436    /// Return a const reference to the edge map storing the found flow.
[2440]437    ///
438    /// \pre \ref run() must be called before using this function.
439    const FlowMap& flowMap() const {
[2581]440      return *_flow;
[2440]441    }
442
[2620]443    /// \brief Return a const reference to the node map storing the
[2574]444    /// found potentials (the dual solution).
[2440]445    ///
[2620]446    /// Return a const reference to the node map storing the found
[2574]447    /// potentials (the dual solution).
[2440]448    ///
449    /// \pre \ref run() must be called before using this function.
450    const PotentialMap& potentialMap() const {
[2581]451      return *_potential;
452    }
453
[2620]454    /// \brief Return the flow on the given edge.
[2581]455    ///
[2620]456    /// Return the flow on the given edge.
[2581]457    ///
458    /// \pre \ref run() must be called before using this function.
459    Capacity flow(const Edge& edge) const {
460      return (*_flow)[edge];
461    }
462
[2620]463    /// \brief Return the potential of the given node.
[2581]464    ///
[2620]465    /// Return the potential of the given node.
[2581]466    ///
467    /// \pre \ref run() must be called before using this function.
468    Cost potential(const Node& node) const {
469      return (*_potential)[node];
[2440]470    }
471
[2620]472    /// \brief Return the total cost of the found flow.
[2440]473    ///
[2620]474    /// Return the total cost of the found flow. The complexity of the
[2440]475    /// function is \f$ O(e) \f$.
476    ///
477    /// \pre \ref run() must be called before using this function.
478    Cost totalCost() const {
479      Cost c = 0;
[2574]480      for (EdgeIt e(_graph); e != INVALID; ++e)
[2581]481        c += (*_flow)[e] * _cost[e];
[2440]482      return c;
483    }
484
[2581]485    /// @}
486
[2574]487  private:
[2440]488
[2620]489    /// Initialize the algorithm.
[2574]490    bool init(bool scaling) {
491      if (!_valid_supply) return false;
[2581]492
493      // Initializing maps
494      if (!_flow) {
495        _flow = new FlowMap(_graph);
496        _local_flow = true;
497      }
498      if (!_potential) {
499        _potential = new PotentialMap(_graph);
500        _local_potential = true;
501      }
502      for (EdgeIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
503      for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
[2440]504
[2581]505      _dijkstra = new ResidualDijkstra( _graph, *_flow, _res_cap, _cost,
506                                        _excess, *_potential, _pred );
507
508      // Initializing delta value
[2574]509      if (scaling) {
[2535]510        // With scaling
511        Supply max_sup = 0, max_dem = 0;
[2574]512        for (NodeIt n(_graph); n != INVALID; ++n) {
513          if ( _supply[n] > max_sup) max_sup =  _supply[n];
514          if (-_supply[n] > max_dem) max_dem = -_supply[n];
[2535]515        }
[2588]516        Capacity max_cap = 0;
517        for (EdgeIt e(_graph); e != INVALID; ++e) {
518          if (_capacity[e] > max_cap) max_cap = _capacity[e];
519        }
520        max_sup = std::min(std::min(max_sup, max_dem), max_cap);
[2574]521        _phase_num = 0;
522        for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2)
523          ++_phase_num;
[2535]524      } else {
525        // Without scaling
[2574]526        _delta = 1;
[2440]527      }
[2581]528
[2440]529      return true;
530    }
531
[2535]532    bool start() {
[2574]533      if (_delta > 1)
[2535]534        return startWithScaling();
535      else
536        return startWithoutScaling();
537    }
538
[2620]539    /// Execute the capacity scaling algorithm.
[2535]540    bool startWithScaling() {
541      // Processing capacity scaling phases
542      Node s, t;
543      int phase_cnt = 0;
544      int factor = 4;
545      while (true) {
546        // Saturating all edges not satisfying the optimality condition
[2574]547        for (EdgeIt e(_graph); e != INVALID; ++e) {
548          Node u = _graph.source(e), v = _graph.target(e);
[2581]549          Cost c = _cost[e] + (*_potential)[u] - (*_potential)[v];
[2574]550          if (c < 0 && _res_cap[e] >= _delta) {
551            _excess[u] -= _res_cap[e];
552            _excess[v] += _res_cap[e];
[2581]553            (*_flow)[e] = _capacity[e];
[2574]554            _res_cap[e] = 0;
[2535]555          }
[2581]556          else if (c > 0 && (*_flow)[e] >= _delta) {
557            _excess[u] += (*_flow)[e];
558            _excess[v] -= (*_flow)[e];
559            (*_flow)[e] = 0;
[2574]560            _res_cap[e] = _capacity[e];
[2535]561          }
562        }
563
564        // Finding excess nodes and deficit nodes
[2574]565        _excess_nodes.clear();
566        _deficit_nodes.clear();
567        for (NodeIt n(_graph); n != INVALID; ++n) {
568          if (_excess[n] >=  _delta) _excess_nodes.push_back(n);
569          if (_excess[n] <= -_delta) _deficit_nodes.push_back(n);
[2535]570        }
[2620]571        int next_node = 0, next_def_node = 0;
[2535]572
573        // Finding augmenting shortest paths
[2574]574        while (next_node < int(_excess_nodes.size())) {
[2535]575          // Checking deficit nodes
[2574]576          if (_delta > 1) {
[2535]577            bool delta_deficit = false;
[2620]578            for ( ; next_def_node < int(_deficit_nodes.size());
579                    ++next_def_node ) {
580              if (_excess[_deficit_nodes[next_def_node]] <= -_delta) {
[2535]581                delta_deficit = true;
582                break;
583              }
584            }
585            if (!delta_deficit) break;
586          }
587
588          // Running Dijkstra
[2574]589          s = _excess_nodes[next_node];
[2581]590          if ((t = _dijkstra->run(s, _delta)) == INVALID) {
[2574]591            if (_delta > 1) {
[2535]592              ++next_node;
593              continue;
594            }
595            return false;
596          }
597
598          // Augmenting along a shortest path from s to t.
[2588]599          Capacity d = std::min(_excess[s], -_excess[t]);
[2535]600          Node u = t;
601          Edge e;
[2574]602          if (d > _delta) {
603            while ((e = _pred[u]) != INVALID) {
[2535]604              Capacity rc;
[2574]605              if (u == _graph.target(e)) {
606                rc = _res_cap[e];
607                u = _graph.source(e);
[2535]608              } else {
[2581]609                rc = (*_flow)[e];
[2574]610                u = _graph.target(e);
[2535]611              }
612              if (rc < d) d = rc;
613            }
614          }
615          u = t;
[2574]616          while ((e = _pred[u]) != INVALID) {
617            if (u == _graph.target(e)) {
[2581]618              (*_flow)[e] += d;
[2574]619              _res_cap[e] -= d;
620              u = _graph.source(e);
[2535]621            } else {
[2581]622              (*_flow)[e] -= d;
[2574]623              _res_cap[e] += d;
624              u = _graph.target(e);
[2535]625            }
626          }
[2574]627          _excess[s] -= d;
628          _excess[t] += d;
[2535]629
[2574]630          if (_excess[s] < _delta) ++next_node;
[2535]631        }
632
[2574]633        if (_delta == 1) break;
634        if (++phase_cnt > _phase_num / 4) factor = 2;
635        _delta = _delta <= factor ? 1 : _delta / factor;
[2535]636      }
637
[2556]638      // Handling non-zero lower bounds
[2574]639      if (_lower) {
640        for (EdgeIt e(_graph); e != INVALID; ++e)
[2581]641          (*_flow)[e] += (*_lower)[e];
[2535]642      }
643      return true;
644    }
645
[2620]646    /// Execute the successive shortest path algorithm.
[2535]647    bool startWithoutScaling() {
[2440]648      // Finding excess nodes
[2574]649      for (NodeIt n(_graph); n != INVALID; ++n)
650        if (_excess[n] > 0) _excess_nodes.push_back(n);
651      if (_excess_nodes.size() == 0) return true;
[2556]652      int next_node = 0;
[2440]653
[2457]654      // Finding shortest paths
[2535]655      Node s, t;
[2574]656      while ( _excess[_excess_nodes[next_node]] > 0 ||
657              ++next_node < int(_excess_nodes.size()) )
[2440]658      {
[2535]659        // Running Dijkstra
[2574]660        s = _excess_nodes[next_node];
[2589]661        if ((t = _dijkstra->run(s)) == INVALID) return false;
[2440]662
[2535]663        // Augmenting along a shortest path from s to t
[2588]664        Capacity d = std::min(_excess[s], -_excess[t]);
[2535]665        Node u = t;
666        Edge e;
[2588]667        if (d > 1) {
668          while ((e = _pred[u]) != INVALID) {
669            Capacity rc;
670            if (u == _graph.target(e)) {
671              rc = _res_cap[e];
672              u = _graph.source(e);
673            } else {
674              rc = (*_flow)[e];
675              u = _graph.target(e);
676            }
677            if (rc < d) d = rc;
[2535]678          }
679        }
680        u = t;
[2574]681        while ((e = _pred[u]) != INVALID) {
682          if (u == _graph.target(e)) {
[2581]683            (*_flow)[e] += d;
[2574]684            _res_cap[e] -= d;
685            u = _graph.source(e);
[2535]686          } else {
[2581]687            (*_flow)[e] -= d;
[2574]688            _res_cap[e] += d;
689            u = _graph.target(e);
[2535]690          }
691        }
[2574]692        _excess[s] -= d;
693        _excess[t] += d;
[2440]694      }
695
[2556]696      // Handling non-zero lower bounds
[2574]697      if (_lower) {
698        for (EdgeIt e(_graph); e != INVALID; ++e)
[2581]699          (*_flow)[e] += (*_lower)[e];
[2440]700      }
701      return true;
702    }
703
704  }; //class CapacityScaling
705
706  ///@}
707
708} //namespace lemon
709
710#endif //LEMON_CAPACITY_SCALING_H
Note: See TracBrowser for help on using the repository browser.