COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/cost_scaling.h @ 2625:c51b320bc51c

Last change on this file since 2625:c51b320bc51c was 2625:c51b320bc51c, checked in by Peter Kovacs, 11 years ago

Major improvement in the cost scaling algorithm

  • Add a new variant that use the partial augment-relabel method.
  • Use this method instead of push-relabel by default.
  • Use the "Early Termination" heuristic instead of "Price Refinement".

Using the new method and heuristic the algorithm proved to be
2-2.5 times faster on all input files.

File size: 28.5 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2008
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_COST_SCALING_H
20#define LEMON_COST_SCALING_H
21
22/// \ingroup min_cost_flow
23/// \file
24/// \brief Cost scaling algorithm for finding a minimum cost flow.
25
26#include <deque>
27#include <lemon/graph_adaptor.h>
28#include <lemon/graph_utils.h>
29#include <lemon/maps.h>
30#include <lemon/math.h>
31
32#include <lemon/circulation.h>
33#include <lemon/bellman_ford.h>
34
35namespace lemon {
36
37  /// \addtogroup min_cost_flow
38  /// @{
39
40  /// \brief Implementation of the cost scaling algorithm for finding a
41  /// minimum cost flow.
42  ///
43  /// \ref CostScaling implements the cost scaling algorithm performing
44  /// augment/push and relabel operations for finding a minimum cost
45  /// flow.
46  ///
47  /// \tparam Graph The directed graph type the algorithm runs on.
48  /// \tparam LowerMap The type of the lower bound map.
49  /// \tparam CapacityMap The type of the capacity (upper bound) map.
50  /// \tparam CostMap The type of the cost (length) map.
51  /// \tparam SupplyMap The type of the supply map.
52  ///
53  /// \warning
54  /// - Edge capacities and costs should be \e non-negative \e integers.
55  /// - Supply values should be \e signed \e integers.
56  /// - The value types of the maps should be convertible to each other.
57  /// - \c CostMap::Value must be signed type.
58  ///
59  /// \note Edge costs are multiplied with the number of nodes during
60  /// the algorithm so overflow problems may arise more easily than with
61  /// other minimum cost flow algorithms.
62  /// If it is available, <tt>long long int</tt> type is used instead of
63  /// <tt>long int</tt> in the inside computations.
64  ///
65  /// \author Peter Kovacs
66  template < typename Graph,
67             typename LowerMap = typename Graph::template EdgeMap<int>,
68             typename CapacityMap = typename Graph::template EdgeMap<int>,
69             typename CostMap = typename Graph::template EdgeMap<int>,
70             typename SupplyMap = typename Graph::template NodeMap<int> >
71  class CostScaling
72  {
73    GRAPH_TYPEDEFS(typename Graph);
74
75    typedef typename CapacityMap::Value Capacity;
76    typedef typename CostMap::Value Cost;
77    typedef typename SupplyMap::Value Supply;
78    typedef typename Graph::template EdgeMap<Capacity> CapacityEdgeMap;
79    typedef typename Graph::template NodeMap<Supply> SupplyNodeMap;
80
81    typedef ResGraphAdaptor< const Graph, Capacity,
82                             CapacityEdgeMap, CapacityEdgeMap > ResGraph;
83    typedef typename ResGraph::Edge ResEdge;
84
85#if defined __GNUC__ && !defined __STRICT_ANSI__
86    typedef long long int LCost;
87#else
88    typedef long int LCost;
89#endif
90    typedef typename Graph::template EdgeMap<LCost> LargeCostMap;
91
92  public:
93
94    /// The type of the flow map.
95    typedef typename Graph::template EdgeMap<Capacity> FlowMap;
96    /// The type of the potential map.
97    typedef typename Graph::template NodeMap<LCost> PotentialMap;
98
99  private:
100
101    /// \brief Map adaptor class for handling residual edge costs.
102    ///
103    /// Map adaptor class for handling residual edge costs.
104    template <typename Map>
105    class ResidualCostMap : public MapBase<ResEdge, typename Map::Value>
106    {
107    private:
108
109      const Map &_cost_map;
110
111    public:
112
113      ///\e
114      ResidualCostMap(const Map &cost_map) :
115        _cost_map(cost_map) {}
116
117      ///\e
118      inline typename Map::Value operator[](const ResEdge &e) const {
119        return ResGraph::forward(e) ? _cost_map[e] : -_cost_map[e];
120      }
121
122    }; //class ResidualCostMap
123
124    /// \brief Map adaptor class for handling reduced edge costs.
125    ///
126    /// Map adaptor class for handling reduced edge costs.
127    class ReducedCostMap : public MapBase<Edge, LCost>
128    {
129    private:
130
131      const Graph &_gr;
132      const LargeCostMap &_cost_map;
133      const PotentialMap &_pot_map;
134
135    public:
136
137      ///\e
138      ReducedCostMap( const Graph &gr,
139                      const LargeCostMap &cost_map,
140                      const PotentialMap &pot_map ) :
141        _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
142
143      ///\e
144      inline LCost operator[](const Edge &e) const {
145        return _cost_map[e] + _pot_map[_gr.source(e)]
146                            - _pot_map[_gr.target(e)];
147      }
148
149    }; //class ReducedCostMap
150
151  private:
152
153    // The directed graph the algorithm runs on
154    const Graph &_graph;
155    // The original lower bound map
156    const LowerMap *_lower;
157    // The modified capacity map
158    CapacityEdgeMap _capacity;
159    // The original cost map
160    const CostMap &_orig_cost;
161    // The scaled cost map
162    LargeCostMap _cost;
163    // The modified supply map
164    SupplyNodeMap _supply;
165    bool _valid_supply;
166
167    // Edge map of the current flow
168    FlowMap *_flow;
169    bool _local_flow;
170    // Node map of the current potentials
171    PotentialMap *_potential;
172    bool _local_potential;
173
174    // The residual cost map
175    ResidualCostMap<LargeCostMap> _res_cost;
176    // The residual graph
177    ResGraph *_res_graph;
178    // The reduced cost map
179    ReducedCostMap *_red_cost;
180    // The excess map
181    SupplyNodeMap _excess;
182    // The epsilon parameter used for cost scaling
183    LCost _epsilon;
184    // The scaling factor
185    int _alpha;
186
187  public:
188
189    /// \brief General constructor (with lower bounds).
190    ///
191    /// General constructor (with lower bounds).
192    ///
193    /// \param graph The directed graph the algorithm runs on.
194    /// \param lower The lower bounds of the edges.
195    /// \param capacity The capacities (upper bounds) of the edges.
196    /// \param cost The cost (length) values of the edges.
197    /// \param supply The supply values of the nodes (signed).
198    CostScaling( const Graph &graph,
199                 const LowerMap &lower,
200                 const CapacityMap &capacity,
201                 const CostMap &cost,
202                 const SupplyMap &supply ) :
203      _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
204      _cost(graph), _supply(graph), _flow(NULL), _local_flow(false),
205      _potential(NULL), _local_potential(false), _res_cost(_cost),
206      _res_graph(NULL), _red_cost(NULL), _excess(graph, 0)
207    {
208      // Remove non-zero lower bounds
209      _capacity = subMap(capacity, lower);
210      Supply sum = 0;
211      for (NodeIt n(_graph); n != INVALID; ++n) {
212        Supply s = supply[n];
213        for (InEdgeIt e(_graph, n); e != INVALID; ++e)
214          s += lower[e];
215        for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
216          s -= lower[e];
217        _supply[n] = s;
218        sum += s;
219      }
220      _valid_supply = sum == 0;
221    }
222
223    /// \brief General constructor (without lower bounds).
224    ///
225    /// General constructor (without lower bounds).
226    ///
227    /// \param graph The directed graph the algorithm runs on.
228    /// \param capacity The capacities (upper bounds) of the edges.
229    /// \param cost The cost (length) values of the edges.
230    /// \param supply The supply values of the nodes (signed).
231    CostScaling( const Graph &graph,
232                 const CapacityMap &capacity,
233                 const CostMap &cost,
234                 const SupplyMap &supply ) :
235      _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
236      _cost(graph), _supply(supply), _flow(NULL), _local_flow(false),
237      _potential(NULL), _local_potential(false), _res_cost(_cost),
238      _res_graph(NULL), _red_cost(NULL), _excess(graph, 0)
239    {
240      // Check the sum of supply values
241      Supply sum = 0;
242      for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
243      _valid_supply = sum == 0;
244    }
245
246    /// \brief Simple constructor (with lower bounds).
247    ///
248    /// Simple constructor (with lower bounds).
249    ///
250    /// \param graph The directed graph the algorithm runs on.
251    /// \param lower The lower bounds of the edges.
252    /// \param capacity The capacities (upper bounds) of the edges.
253    /// \param cost The cost (length) values of the edges.
254    /// \param s The source node.
255    /// \param t The target node.
256    /// \param flow_value The required amount of flow from node \c s
257    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
258    CostScaling( const Graph &graph,
259                 const LowerMap &lower,
260                 const CapacityMap &capacity,
261                 const CostMap &cost,
262                 Node s, Node t,
263                 Supply flow_value ) :
264      _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
265      _cost(graph), _supply(graph), _flow(NULL), _local_flow(false),
266      _potential(NULL), _local_potential(false), _res_cost(_cost),
267      _res_graph(NULL), _red_cost(NULL), _excess(graph, 0)
268    {
269      // Remove nonzero lower bounds
270      _capacity = subMap(capacity, lower);
271      for (NodeIt n(_graph); n != INVALID; ++n) {
272        Supply sum = 0;
273        if (n == s) sum =  flow_value;
274        if (n == t) sum = -flow_value;
275        for (InEdgeIt e(_graph, n); e != INVALID; ++e)
276          sum += lower[e];
277        for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
278          sum -= lower[e];
279        _supply[n] = sum;
280      }
281      _valid_supply = true;
282    }
283
284    /// \brief Simple constructor (without lower bounds).
285    ///
286    /// Simple constructor (without lower bounds).
287    ///
288    /// \param graph The directed graph the algorithm runs on.
289    /// \param capacity The capacities (upper bounds) of the edges.
290    /// \param cost The cost (length) values of the edges.
291    /// \param s The source node.
292    /// \param t The target node.
293    /// \param flow_value The required amount of flow from node \c s
294    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
295    CostScaling( const Graph &graph,
296                 const CapacityMap &capacity,
297                 const CostMap &cost,
298                 Node s, Node t,
299                 Supply flow_value ) :
300      _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
301      _cost(graph), _supply(graph, 0), _flow(NULL), _local_flow(false),
302      _potential(NULL), _local_potential(false), _res_cost(_cost),
303      _res_graph(NULL), _red_cost(NULL), _excess(graph, 0)
304    {
305      _supply[s] =  flow_value;
306      _supply[t] = -flow_value;
307      _valid_supply = true;
308    }
309
310    /// Destructor.
311    ~CostScaling() {
312      if (_local_flow) delete _flow;
313      if (_local_potential) delete _potential;
314      delete _res_graph;
315      delete _red_cost;
316    }
317
318    /// \brief Set the flow map.
319    ///
320    /// Set the flow map.
321    ///
322    /// \return \c (*this)
323    CostScaling& flowMap(FlowMap &map) {
324      if (_local_flow) {
325        delete _flow;
326        _local_flow = false;
327      }
328      _flow = &map;
329      return *this;
330    }
331
332    /// \brief Set the potential map.
333    ///
334    /// Set the potential map.
335    ///
336    /// \return \c (*this)
337    CostScaling& potentialMap(PotentialMap &map) {
338      if (_local_potential) {
339        delete _potential;
340        _local_potential = false;
341      }
342      _potential = &map;
343      return *this;
344    }
345
346    /// \name Execution control
347
348    /// @{
349
350    /// \brief Run the algorithm.
351    ///
352    /// Run the algorithm.
353    ///
354    /// \param partial_augment By default the algorithm performs
355    /// partial augment and relabel operations in the cost scaling
356    /// phases. Set this parameter to \c false for using local push and
357    /// relabel operations instead.
358    ///
359    /// \return \c true if a feasible flow can be found.
360    bool run(bool partial_augment = true) {
361      if (partial_augment) {
362        return init() && startPartialAugment();
363      } else {
364        return init() && startPushRelabel();
365      }
366    }
367
368    /// @}
369
370    /// \name Query Functions
371    /// The result of the algorithm can be obtained using these
372    /// functions.\n
373    /// \ref lemon::CostScaling::run() "run()" must be called before
374    /// using them.
375
376    /// @{
377
378    /// \brief Return a const reference to the edge map storing the
379    /// found flow.
380    ///
381    /// Return a const reference to the edge map storing the found flow.
382    ///
383    /// \pre \ref run() must be called before using this function.
384    const FlowMap& flowMap() const {
385      return *_flow;
386    }
387
388    /// \brief Return a const reference to the node map storing the
389    /// found potentials (the dual solution).
390    ///
391    /// Return a const reference to the node map storing the found
392    /// potentials (the dual solution).
393    ///
394    /// \pre \ref run() must be called before using this function.
395    const PotentialMap& potentialMap() const {
396      return *_potential;
397    }
398
399    /// \brief Return the flow on the given edge.
400    ///
401    /// Return the flow on the given edge.
402    ///
403    /// \pre \ref run() must be called before using this function.
404    Capacity flow(const Edge& edge) const {
405      return (*_flow)[edge];
406    }
407
408    /// \brief Return the potential of the given node.
409    ///
410    /// Return the potential of the given node.
411    ///
412    /// \pre \ref run() must be called before using this function.
413    Cost potential(const Node& node) const {
414      return (*_potential)[node];
415    }
416
417    /// \brief Return the total cost of the found flow.
418    ///
419    /// Return the total cost of the found flow. The complexity of the
420    /// function is \f$ O(e) \f$.
421    ///
422    /// \pre \ref run() must be called before using this function.
423    Cost totalCost() const {
424      Cost c = 0;
425      for (EdgeIt e(_graph); e != INVALID; ++e)
426        c += (*_flow)[e] * _orig_cost[e];
427      return c;
428    }
429
430    /// @}
431
432  private:
433
434    /// Initialize the algorithm.
435    bool init() {
436      if (!_valid_supply) return false;
437      // The scaling factor
438      _alpha = 8;
439
440      // Initialize flow and potential maps
441      if (!_flow) {
442        _flow = new FlowMap(_graph);
443        _local_flow = true;
444      }
445      if (!_potential) {
446        _potential = new PotentialMap(_graph);
447        _local_potential = true;
448      }
449
450      _red_cost = new ReducedCostMap(_graph, _cost, *_potential);
451      _res_graph = new ResGraph(_graph, _capacity, *_flow);
452
453      // Initialize the scaled cost map and the epsilon parameter
454      Cost max_cost = 0;
455      int node_num = countNodes(_graph);
456      for (EdgeIt e(_graph); e != INVALID; ++e) {
457        _cost[e] = LCost(_orig_cost[e]) * node_num * _alpha;
458        if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e];
459      }
460      _epsilon = max_cost * node_num;
461
462      // Find a feasible flow using Circulation
463      Circulation< Graph, ConstMap<Edge, Capacity>, CapacityEdgeMap,
464                   SupplyMap >
465        circulation( _graph, constMap<Edge>(Capacity(0)), _capacity,
466                     _supply );
467      return circulation.flowMap(*_flow).run();
468    }
469
470    /// Execute the algorithm performing partial augmentation and
471    /// relabel operations.
472    bool startPartialAugment() {
473      // Paramters for heuristics
474      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
475      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
476      // Maximum augment path length
477      const int MAX_PATH_LENGTH = 4;
478
479      // Variables
480      typename Graph::template NodeMap<Edge> pred_edge(_graph);
481      typename Graph::template NodeMap<bool> forward(_graph);
482      typename Graph::template NodeMap<OutEdgeIt> next_out(_graph);
483      typename Graph::template NodeMap<InEdgeIt> next_in(_graph);
484      typename Graph::template NodeMap<bool> next_dir(_graph);
485      std::deque<Node> active_nodes;
486      std::vector<Node> path_nodes;
487
488      int node_num = countNodes(_graph);
489      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
490                                        1 : _epsilon / _alpha )
491      {
492        // "Early Termination" heuristic: use Bellman-Ford algorithm
493        // to check if the current flow is optimal
494        if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
495          typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
496          ShiftCostMap shift_cost(_res_cost, 1);
497          BellmanFord<ResGraph, ShiftCostMap> bf(*_res_graph, shift_cost);
498          bf.init(0);
499          bool done = false;
500          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
501          for (int i = 0; i < K && !done; ++i)
502            done = bf.processNextWeakRound();
503          if (done) break;
504        }
505
506        // Saturate edges not satisfying the optimality condition
507        Capacity delta;
508        for (EdgeIt e(_graph); e != INVALID; ++e) {
509          if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
510            delta = _capacity[e] - (*_flow)[e];
511            _excess[_graph.source(e)] -= delta;
512            _excess[_graph.target(e)] += delta;
513            (*_flow)[e] = _capacity[e];
514          }
515          if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
516            _excess[_graph.target(e)] -= (*_flow)[e];
517            _excess[_graph.source(e)] += (*_flow)[e];
518            (*_flow)[e] = 0;
519          }
520        }
521
522        // Find active nodes (i.e. nodes with positive excess)
523        for (NodeIt n(_graph); n != INVALID; ++n) {
524          if (_excess[n] > 0) active_nodes.push_back(n);
525        }
526
527        // Initialize the next edge maps
528        for (NodeIt n(_graph); n != INVALID; ++n) {
529          next_out[n] = OutEdgeIt(_graph, n);
530          next_in[n] = InEdgeIt(_graph, n);
531          next_dir[n] = true;
532        }
533
534        // Perform partial augment and relabel operations
535        while (active_nodes.size() > 0) {
536          // Select an active node (FIFO selection)
537          if (_excess[active_nodes[0]] <= 0) {
538            active_nodes.pop_front();
539            continue;
540          }
541          Node start = active_nodes[0];
542          path_nodes.clear();
543          path_nodes.push_back(start);
544
545          // Find an augmenting path from the start node
546          Node u, tip = start;
547          LCost min_red_cost;
548          while ( _excess[tip] >= 0 &&
549                  int(path_nodes.size()) <= MAX_PATH_LENGTH )
550          {
551            if (next_dir[tip]) {
552              for (OutEdgeIt e = next_out[tip]; e != INVALID; ++e) {
553                if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
554                  u = _graph.target(e);
555                  pred_edge[u] = e;
556                  forward[u] = true;
557                  next_out[tip] = e;
558                  tip = u;
559                  path_nodes.push_back(tip);
560                  goto next_step;
561                }
562              }
563              next_dir[tip] = false;
564            }
565            for (InEdgeIt e = next_in[tip]; e != INVALID; ++e) {
566              if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
567                u = _graph.source(e);
568                pred_edge[u] = e;
569                forward[u] = false;
570                next_in[tip] = e;
571                tip = u;
572                path_nodes.push_back(tip);
573                goto next_step;
574              }
575            }
576
577            // Relabel tip node
578            min_red_cost = std::numeric_limits<LCost>::max() / 2;
579            for (OutEdgeIt oe(_graph, tip); oe != INVALID; ++oe) {
580              if ( _capacity[oe] - (*_flow)[oe] > 0 &&
581                   (*_red_cost)[oe] < min_red_cost )
582                min_red_cost = (*_red_cost)[oe];
583            }
584            for (InEdgeIt ie(_graph, tip); ie != INVALID; ++ie) {
585              if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
586                min_red_cost = -(*_red_cost)[ie];
587            }
588            (*_potential)[tip] -= min_red_cost + _epsilon;
589
590            // Reset the next edge maps
591            next_out[tip] = OutEdgeIt(_graph, tip);
592            next_in[tip] = InEdgeIt(_graph, tip);
593            next_dir[tip] = true;
594
595            // Step back
596            if (tip != start) {
597              path_nodes.pop_back();
598              tip = path_nodes[path_nodes.size()-1];
599            }
600
601          next_step:
602            continue;
603          }
604
605          // Augment along the found path (as much flow as possible)
606          Capacity delta;
607          for (int i = 1; i < int(path_nodes.size()); ++i) {
608            u = path_nodes[i];
609            delta = forward[u] ?
610              _capacity[pred_edge[u]] - (*_flow)[pred_edge[u]] :
611              (*_flow)[pred_edge[u]];
612            delta = std::min(delta, _excess[path_nodes[i-1]]);
613            (*_flow)[pred_edge[u]] += forward[u] ? delta : -delta;
614            _excess[path_nodes[i-1]] -= delta;
615            _excess[u] += delta;
616            if (_excess[u] > 0 && _excess[u] <= delta) active_nodes.push_back(u);
617          }
618        }
619      }
620
621      // Compute node potentials for the original costs
622      ResidualCostMap<CostMap> res_cost(_orig_cost);
623      BellmanFord< ResGraph, ResidualCostMap<CostMap> >
624        bf(*_res_graph, res_cost);
625      bf.init(0); bf.start();
626      for (NodeIt n(_graph); n != INVALID; ++n)
627        (*_potential)[n] = bf.dist(n);
628
629      // Handle non-zero lower bounds
630      if (_lower) {
631        for (EdgeIt e(_graph); e != INVALID; ++e)
632          (*_flow)[e] += (*_lower)[e];
633      }
634      return true;
635    }
636
637    /// Execute the algorithm performing push and relabel operations.
638    bool startPushRelabel() {
639      // Paramters for heuristics
640      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
641      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
642
643      typename Graph::template NodeMap<bool> hyper(_graph, false);
644      typename Graph::template NodeMap<Edge> pred_edge(_graph);
645      typename Graph::template NodeMap<bool> forward(_graph);
646      typename Graph::template NodeMap<OutEdgeIt> next_out(_graph);
647      typename Graph::template NodeMap<InEdgeIt> next_in(_graph);
648      typename Graph::template NodeMap<bool> next_dir(_graph);
649      std::deque<Node> active_nodes;
650
651      int node_num = countNodes(_graph);
652      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
653                                        1 : _epsilon / _alpha )
654      {
655        // "Early Termination" heuristic: use Bellman-Ford algorithm
656        // to check if the current flow is optimal
657        if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
658          typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
659          ShiftCostMap shift_cost(_res_cost, 1);
660          BellmanFord<ResGraph, ShiftCostMap> bf(*_res_graph, shift_cost);
661          bf.init(0);
662          bool done = false;
663          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
664          for (int i = 0; i < K && !done; ++i)
665            done = bf.processNextWeakRound();
666          if (done) break;
667        }
668
669        // Saturate edges not satisfying the optimality condition
670        Capacity delta;
671        for (EdgeIt e(_graph); e != INVALID; ++e) {
672          if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
673            delta = _capacity[e] - (*_flow)[e];
674            _excess[_graph.source(e)] -= delta;
675            _excess[_graph.target(e)] += delta;
676            (*_flow)[e] = _capacity[e];
677          }
678          if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
679            _excess[_graph.target(e)] -= (*_flow)[e];
680            _excess[_graph.source(e)] += (*_flow)[e];
681            (*_flow)[e] = 0;
682          }
683        }
684
685        // Find active nodes (i.e. nodes with positive excess)
686        for (NodeIt n(_graph); n != INVALID; ++n) {
687          if (_excess[n] > 0) active_nodes.push_back(n);
688        }
689
690        // Initialize the next edge maps
691        for (NodeIt n(_graph); n != INVALID; ++n) {
692          next_out[n] = OutEdgeIt(_graph, n);
693          next_in[n] = InEdgeIt(_graph, n);
694          next_dir[n] = true;
695        }
696
697        // Perform push and relabel operations
698        while (active_nodes.size() > 0) {
699          // Select an active node (FIFO selection)
700          Node n = active_nodes[0], t;
701          bool relabel_enabled = true;
702
703          // Perform push operations if there are admissible edges
704          if (_excess[n] > 0 && next_dir[n]) {
705            OutEdgeIt e = next_out[n];
706            for ( ; e != INVALID; ++e) {
707              if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
708                delta = std::min(_capacity[e] - (*_flow)[e], _excess[n]);
709                t = _graph.target(e);
710
711                // Push-look-ahead heuristic
712                Capacity ahead = -_excess[t];
713                for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
714                  if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
715                    ahead += _capacity[oe] - (*_flow)[oe];
716                }
717                for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
718                  if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
719                    ahead += (*_flow)[ie];
720                }
721                if (ahead < 0) ahead = 0;
722
723                // Push flow along the edge
724                if (ahead < delta) {
725                  (*_flow)[e] += ahead;
726                  _excess[n] -= ahead;
727                  _excess[t] += ahead;
728                  active_nodes.push_front(t);
729                  hyper[t] = true;
730                  relabel_enabled = false;
731                  break;
732                } else {
733                  (*_flow)[e] += delta;
734                  _excess[n] -= delta;
735                  _excess[t] += delta;
736                  if (_excess[t] > 0 && _excess[t] <= delta)
737                    active_nodes.push_back(t);
738                }
739
740                if (_excess[n] == 0) break;
741              }
742            }
743            if (e != INVALID) {
744              next_out[n] = e;
745            } else {
746              next_dir[n] = false;
747            }
748          }
749
750          if (_excess[n] > 0 && !next_dir[n]) {
751            InEdgeIt e = next_in[n];
752            for ( ; e != INVALID; ++e) {
753              if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
754                delta = std::min((*_flow)[e], _excess[n]);
755                t = _graph.source(e);
756
757                // Push-look-ahead heuristic
758                Capacity ahead = -_excess[t];
759                for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
760                  if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
761                    ahead += _capacity[oe] - (*_flow)[oe];
762                }
763                for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
764                  if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
765                    ahead += (*_flow)[ie];
766                }
767                if (ahead < 0) ahead = 0;
768
769                // Push flow along the edge
770                if (ahead < delta) {
771                  (*_flow)[e] -= ahead;
772                  _excess[n] -= ahead;
773                  _excess[t] += ahead;
774                  active_nodes.push_front(t);
775                  hyper[t] = true;
776                  relabel_enabled = false;
777                  break;
778                } else {
779                  (*_flow)[e] -= delta;
780                  _excess[n] -= delta;
781                  _excess[t] += delta;
782                  if (_excess[t] > 0 && _excess[t] <= delta)
783                    active_nodes.push_back(t);
784                }
785
786                if (_excess[n] == 0) break;
787              }
788            }
789            next_in[n] = e;
790          }
791
792          // Relabel the node if it is still active (or hyper)
793          if (relabel_enabled && (_excess[n] > 0 || hyper[n])) {
794            LCost min_red_cost = std::numeric_limits<LCost>::max() / 2;
795            for (OutEdgeIt oe(_graph, n); oe != INVALID; ++oe) {
796              if ( _capacity[oe] - (*_flow)[oe] > 0 &&
797                   (*_red_cost)[oe] < min_red_cost )
798                min_red_cost = (*_red_cost)[oe];
799            }
800            for (InEdgeIt ie(_graph, n); ie != INVALID; ++ie) {
801              if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
802                min_red_cost = -(*_red_cost)[ie];
803            }
804            (*_potential)[n] -= min_red_cost + _epsilon;
805            hyper[n] = false;
806
807            // Reset the next edge maps
808            next_out[n] = OutEdgeIt(_graph, n);
809            next_in[n] = InEdgeIt(_graph, n);
810            next_dir[n] = true;
811          }
812
813          // Remove nodes that are not active nor hyper
814          while ( active_nodes.size() > 0 &&
815                  _excess[active_nodes[0]] <= 0 &&
816                  !hyper[active_nodes[0]] ) {
817            active_nodes.pop_front();
818          }
819        }
820      }
821
822      // Compute node potentials for the original costs
823      ResidualCostMap<CostMap> res_cost(_orig_cost);
824      BellmanFord< ResGraph, ResidualCostMap<CostMap> >
825        bf(*_res_graph, res_cost);
826      bf.init(0); bf.start();
827      for (NodeIt n(_graph); n != INVALID; ++n)
828        (*_potential)[n] = bf.dist(n);
829
830      // Handle non-zero lower bounds
831      if (_lower) {
832        for (EdgeIt e(_graph); e != INVALID; ++e)
833          (*_flow)[e] += (*_lower)[e];
834      }
835      return true;
836    }
837
838  }; //class CostScaling
839
840  ///@}
841
842} //namespace lemon
843
844#endif //LEMON_COST_SCALING_H
Note: See TracBrowser for help on using the repository browser.