COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/cost_scaling.h @ 2588:4d3bc1d04c1d

Last change on this file since 2588:4d3bc1d04c1d was 2588:4d3bc1d04c1d, checked in by Peter Kovacs, 11 years ago

Small improvements in min cost flow files.

File size: 21.2 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///
24/// \file
25/// \brief Cost scaling algorithm for finding a minimum cost flow.
26
27#include <deque>
28#include <lemon/graph_adaptor.h>
29#include <lemon/graph_utils.h>
30#include <lemon/maps.h>
31#include <lemon/math.h>
32
33#include <lemon/circulation.h>
34#include <lemon/bellman_ford.h>
35
36namespace lemon {
37
38  /// \addtogroup min_cost_flow
39  /// @{
40
41  /// \brief Implementation of the cost scaling algorithm for finding a
42  /// minimum cost flow.
43  ///
44  /// \ref CostScaling implements the cost scaling algorithm performing
45  /// generalized push-relabel operations for finding a minimum cost
46  /// flow.
47  ///
48  /// \tparam Graph The directed graph type the algorithm runs on.
49  /// \tparam LowerMap The type of the lower bound map.
50  /// \tparam CapacityMap The type of the capacity (upper bound) map.
51  /// \tparam CostMap The type of the cost (length) map.
52  /// \tparam SupplyMap The type of the supply map.
53  ///
54  /// \warning
55  /// - Edge capacities and costs should be \e non-negative \e integers.
56  /// - Supply values should be \e signed \e integers.
57  /// - The value types of the maps should be convertible to each other.
58  /// - \c CostMap::Value must be signed type.
59  ///
60  /// \note Edge costs are multiplied with the number of nodes during
61  /// the algorithm so overflow problems may arise more easily than with
62  /// other minimum cost flow algorithms.
63  /// If it is available, <tt>long long int</tt> type is used instead of
64  /// <tt>long int</tt> in the inside computations.
65  ///
66  /// \author Peter Kovacs
67
68  template < typename Graph,
69             typename LowerMap = typename Graph::template EdgeMap<int>,
70             typename CapacityMap = typename Graph::template EdgeMap<int>,
71             typename CostMap = typename Graph::template EdgeMap<int>,
72             typename SupplyMap = typename Graph::template NodeMap<int> >
73  class CostScaling
74  {
75    GRAPH_TYPEDEFS(typename Graph);
76
77    typedef typename CapacityMap::Value Capacity;
78    typedef typename CostMap::Value Cost;
79    typedef typename SupplyMap::Value Supply;
80    typedef typename Graph::template EdgeMap<Capacity> CapacityEdgeMap;
81    typedef typename Graph::template NodeMap<Supply> SupplyNodeMap;
82
83    typedef ResGraphAdaptor< const Graph, Capacity,
84                             CapacityEdgeMap, CapacityEdgeMap > ResGraph;
85    typedef typename ResGraph::Edge ResEdge;
86
87#if defined __GNUC__ && !defined __STRICT_ANSI__
88    typedef long long int LCost;
89#else
90    typedef long int LCost;
91#endif
92    typedef typename Graph::template EdgeMap<LCost> LargeCostMap;
93
94  public:
95
96    /// The type of the flow map.
97    typedef typename Graph::template EdgeMap<Capacity> FlowMap;
98    /// The type of the potential map.
99    typedef typename Graph::template NodeMap<LCost> PotentialMap;
100
101  private:
102
103    /// \brief Map adaptor class for handling residual edge costs.
104    ///
105    /// \ref ResidualCostMap is a map adaptor class for handling
106    /// residual edge costs.
107    template <typename Map>
108    class ResidualCostMap : public MapBase<ResEdge, typename Map::Value>
109    {
110    private:
111
112      const Map &_cost_map;
113
114    public:
115
116      ///\e
117      ResidualCostMap(const Map &cost_map) :
118        _cost_map(cost_map) {}
119
120      ///\e
121      typename Map::Value operator[](const ResEdge &e) const {
122        return ResGraph::forward(e) ?  _cost_map[e] : -_cost_map[e];
123      }
124
125    }; //class ResidualCostMap
126
127    /// \brief Map adaptor class for handling reduced edge costs.
128    ///
129    /// \ref ReducedCostMap is a map adaptor class for handling reduced
130    /// edge costs.
131    class ReducedCostMap : public MapBase<Edge, LCost>
132    {
133    private:
134
135      const Graph &_gr;
136      const LargeCostMap &_cost_map;
137      const PotentialMap &_pot_map;
138
139    public:
140
141      ///\e
142      ReducedCostMap( const Graph &gr,
143                      const LargeCostMap &cost_map,
144                      const PotentialMap &pot_map ) :
145        _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
146
147      ///\e
148      LCost operator[](const Edge &e) const {
149        return _cost_map[e] + _pot_map[_gr.source(e)]
150                            - _pot_map[_gr.target(e)];
151      }
152
153    }; //class ReducedCostMap
154
155  private:
156
157    // Scaling factor
158    static const int ALPHA = 4;
159
160    // Paramters for heuristics
161    static const int BF_HEURISTIC_EPSILON_BOUND = 5000;
162    static const int BF_HEURISTIC_BOUND_FACTOR  = 3;
163
164  private:
165
166    // The directed graph the algorithm runs on
167    const Graph &_graph;
168    // The original lower bound map
169    const LowerMap *_lower;
170    // The modified capacity map
171    CapacityEdgeMap _capacity;
172    // The original cost map
173    const CostMap &_orig_cost;
174    // The scaled cost map
175    LargeCostMap _cost;
176    // The modified supply map
177    SupplyNodeMap _supply;
178    bool _valid_supply;
179
180    // Edge map of the current flow
181    FlowMap *_flow;
182    bool _local_flow;
183    // Node map of the current potentials
184    PotentialMap *_potential;
185    bool _local_potential;
186
187    // The residual graph
188    ResGraph *_res_graph;
189    // The residual cost map
190    ResidualCostMap<LargeCostMap> _res_cost;
191    // The reduced cost map
192    ReducedCostMap *_red_cost;
193    // The excess map
194    SupplyNodeMap _excess;
195    // The epsilon parameter used for cost scaling
196    LCost _epsilon;
197
198  public:
199
200    /// \brief General constructor (with lower bounds).
201    ///
202    /// General constructor (with lower bounds).
203    ///
204    /// \param graph The directed graph the algorithm runs on.
205    /// \param lower The lower bounds of the edges.
206    /// \param capacity The capacities (upper bounds) of the edges.
207    /// \param cost The cost (length) values of the edges.
208    /// \param supply The supply values of the nodes (signed).
209    CostScaling( const Graph &graph,
210                 const LowerMap &lower,
211                 const CapacityMap &capacity,
212                 const CostMap &cost,
213                 const SupplyMap &supply ) :
214      _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
215      _cost(graph), _supply(graph), _flow(0), _local_flow(false),
216      _potential(0), _local_potential(false), _res_cost(_cost),
217      _excess(graph, 0)
218    {
219      // Removing non-zero lower bounds
220      _capacity = subMap(capacity, lower);
221      Supply sum = 0;
222      for (NodeIt n(_graph); n != INVALID; ++n) {
223        Supply s = supply[n];
224        for (InEdgeIt e(_graph, n); e != INVALID; ++e)
225          s += lower[e];
226        for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
227          s -= lower[e];
228        _supply[n] = s;
229        sum += s;
230      }
231      _valid_supply = sum == 0;
232    }
233
234    /// \brief General constructor (without lower bounds).
235    ///
236    /// General constructor (without lower bounds).
237    ///
238    /// \param graph The directed graph the algorithm runs on.
239    /// \param capacity The capacities (upper bounds) of the edges.
240    /// \param cost The cost (length) values of the edges.
241    /// \param supply The supply values of the nodes (signed).
242    CostScaling( const Graph &graph,
243                 const CapacityMap &capacity,
244                 const CostMap &cost,
245                 const SupplyMap &supply ) :
246      _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
247      _cost(graph), _supply(supply), _flow(0), _local_flow(false),
248      _potential(0), _local_potential(false), _res_cost(_cost),
249      _excess(graph, 0)
250    {
251      // Checking the sum of supply values
252      Supply sum = 0;
253      for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
254      _valid_supply = sum == 0;
255    }
256
257    /// \brief Simple constructor (with lower bounds).
258    ///
259    /// Simple constructor (with lower bounds).
260    ///
261    /// \param graph The directed graph the algorithm runs on.
262    /// \param lower The lower bounds of the edges.
263    /// \param capacity The capacities (upper bounds) of the edges.
264    /// \param cost The cost (length) values of the edges.
265    /// \param s The source node.
266    /// \param t The target node.
267    /// \param flow_value The required amount of flow from node \c s
268    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
269    CostScaling( const Graph &graph,
270                 const LowerMap &lower,
271                 const CapacityMap &capacity,
272                 const CostMap &cost,
273                 Node s, Node t,
274                 Supply flow_value ) :
275      _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
276      _cost(graph), _supply(graph), _flow(0), _local_flow(false),
277      _potential(0), _local_potential(false), _res_cost(_cost),
278      _excess(graph, 0)
279    {
280      // Removing nonzero lower bounds
281      _capacity = subMap(capacity, lower);
282      for (NodeIt n(_graph); n != INVALID; ++n) {
283        Supply sum = 0;
284        if (n == s) sum =  flow_value;
285        if (n == t) sum = -flow_value;
286        for (InEdgeIt e(_graph, n); e != INVALID; ++e)
287          sum += lower[e];
288        for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
289          sum -= lower[e];
290        _supply[n] = sum;
291      }
292      _valid_supply = true;
293    }
294
295    /// \brief Simple constructor (without lower bounds).
296    ///
297    /// Simple constructor (without lower bounds).
298    ///
299    /// \param graph The directed graph the algorithm runs on.
300    /// \param capacity The capacities (upper bounds) of the edges.
301    /// \param cost The cost (length) values of the edges.
302    /// \param s The source node.
303    /// \param t The target node.
304    /// \param flow_value The required amount of flow from node \c s
305    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
306    CostScaling( const Graph &graph,
307                 const CapacityMap &capacity,
308                 const CostMap &cost,
309                 Node s, Node t,
310                 Supply flow_value ) :
311      _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
312      _cost(graph), _supply(graph, 0), _flow(0), _local_flow(false),
313      _potential(0), _local_potential(false), _res_cost(_cost),
314      _excess(graph, 0)
315    {
316      _supply[s] =  flow_value;
317      _supply[t] = -flow_value;
318      _valid_supply = true;
319    }
320
321    /// Destructor.
322    ~CostScaling() {
323      if (_local_flow) delete _flow;
324      if (_local_potential) delete _potential;
325      delete _res_graph;
326      delete _red_cost;
327    }
328
329    /// \brief Sets the flow map.
330    ///
331    /// Sets the flow map.
332    ///
333    /// \return \c (*this)
334    CostScaling& flowMap(FlowMap &map) {
335      if (_local_flow) {
336        delete _flow;
337        _local_flow = false;
338      }
339      _flow = &map;
340      return *this;
341    }
342
343    /// \brief Sets the potential map.
344    ///
345    /// Sets the potential map.
346    ///
347    /// \return \c (*this)
348    CostScaling& potentialMap(PotentialMap &map) {
349      if (_local_potential) {
350        delete _potential;
351        _local_potential = false;
352      }
353      _potential = &map;
354      return *this;
355    }
356
357    /// \name Execution control
358    /// The only way to execute the algorithm is to call the run()
359    /// function.
360
361    /// @{
362
363    /// \brief Runs the algorithm.
364    ///
365    /// Runs the algorithm.
366    ///
367    /// \return \c true if a feasible flow can be found.
368    bool run() {
369      return init() && start();
370    }
371
372    /// @}
373
374    /// \name Query Functions
375    /// The result of the algorithm can be obtained using these
376    /// functions.
377    /// \n run() must be called before using them.
378
379    /// @{
380
381    /// \brief Returns a const reference to the edge map storing the
382    /// found flow.
383    ///
384    /// Returns a const reference to the edge map storing the found flow.
385    ///
386    /// \pre \ref run() must be called before using this function.
387    const FlowMap& flowMap() const {
388      return *_flow;
389    }
390
391    /// \brief Returns a const reference to the node map storing the
392    /// found potentials (the dual solution).
393    ///
394    /// Returns a const reference to the node map storing the found
395    /// potentials (the dual solution).
396    ///
397    /// \pre \ref run() must be called before using this function.
398    const PotentialMap& potentialMap() const {
399      return *_potential;
400    }
401
402    /// \brief Returns the flow on the given edge.
403    ///
404    /// Returns the flow on the given edge.
405    ///
406    /// \pre \ref run() must be called before using this function.
407    Capacity flow(const Edge& edge) const {
408      return (*_flow)[edge];
409    }
410
411    /// \brief Returns the potential of the given node.
412    ///
413    /// Returns the potential of the given node.
414    ///
415    /// \pre \ref run() must be called before using this function.
416    Cost potential(const Node& node) const {
417      return (*_potential)[node];
418    }
419
420    /// \brief Returns the total cost of the found flow.
421    ///
422    /// Returns the total cost of the found flow. The complexity of the
423    /// function is \f$ O(e) \f$.
424    ///
425    /// \pre \ref run() must be called before using this function.
426    Cost totalCost() const {
427      Cost c = 0;
428      for (EdgeIt e(_graph); e != INVALID; ++e)
429        c += (*_flow)[e] * _orig_cost[e];
430      return c;
431    }
432
433    /// @}
434
435  private:
436
437    /// Initializes the algorithm.
438    bool init() {
439      if (!_valid_supply) return false;
440
441      // Initializing flow and potential maps
442      if (!_flow) {
443        _flow = new FlowMap(_graph);
444        _local_flow = true;
445      }
446      if (!_potential) {
447        _potential = new PotentialMap(_graph);
448        _local_potential = true;
449      }
450
451      _red_cost = new ReducedCostMap(_graph, _cost, *_potential);
452      _res_graph = new ResGraph(_graph, _capacity, *_flow);
453
454      // Initializing the scaled cost map and the epsilon parameter
455      Cost max_cost = 0;
456      int node_num = countNodes(_graph);
457      for (EdgeIt e(_graph); e != INVALID; ++e) {
458        _cost[e] = LCost(_orig_cost[e]) * node_num * ALPHA;
459        if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e];
460      }
461      _epsilon = max_cost * node_num;
462
463      // Finding a feasible flow using Circulation
464      Circulation< Graph, ConstMap<Edge, Capacity>, CapacityEdgeMap,
465                   SupplyMap >
466        circulation( _graph, constMap<Edge>(Capacity(0)), _capacity,
467                     _supply );
468      return circulation.flowMap(*_flow).run();
469    }
470
471
472    /// Executes the algorithm.
473    bool start() {
474      std::deque<Node> active_nodes;
475      typename Graph::template NodeMap<bool> hyper(_graph, false);
476
477      int node_num = countNodes(_graph);
478      for ( ; _epsilon >= 1; _epsilon = _epsilon < ALPHA && _epsilon > 1 ?
479                                        1 : _epsilon / ALPHA )
480      {
481        // Performing price refinement heuristic using Bellman-Ford
482        // algorithm
483        if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
484          typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
485          ShiftCostMap shift_cost(_res_cost, _epsilon);
486          BellmanFord<ResGraph, ShiftCostMap> bf(*_res_graph, shift_cost);
487          bf.init(0);
488          bool done = false;
489          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
490          for (int i = 0; i < K && !done; ++i)
491            done = bf.processNextWeakRound();
492          if (done) {
493            for (NodeIt n(_graph); n != INVALID; ++n)
494              (*_potential)[n] = bf.dist(n);
495            continue;
496          }
497        }
498
499        // Saturating edges not satisfying the optimality condition
500        Capacity delta;
501        for (EdgeIt e(_graph); e != INVALID; ++e) {
502          if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
503            delta = _capacity[e] - (*_flow)[e];
504            _excess[_graph.source(e)] -= delta;
505            _excess[_graph.target(e)] += delta;
506            (*_flow)[e] = _capacity[e];
507          }
508          if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
509            _excess[_graph.target(e)] -= (*_flow)[e];
510            _excess[_graph.source(e)] += (*_flow)[e];
511            (*_flow)[e] = 0;
512          }
513        }
514
515        // Finding active nodes (i.e. nodes with positive excess)
516        for (NodeIt n(_graph); n != INVALID; ++n)
517          if (_excess[n] > 0) active_nodes.push_back(n);
518
519        // Performing push and relabel operations
520        while (active_nodes.size() > 0) {
521          Node n = active_nodes[0], t;
522          bool relabel_enabled = true;
523
524          // Performing push operations if there are admissible edges
525          if (_excess[n] > 0) {
526            for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
527              if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
528                delta = _capacity[e] - (*_flow)[e] <= _excess[n] ?
529                        _capacity[e] - (*_flow)[e] : _excess[n];
530                t = _graph.target(e);
531
532                // Push-look-ahead heuristic
533                Capacity ahead = -_excess[t];
534                for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
535                  if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
536                    ahead += _capacity[oe] - (*_flow)[oe];
537                }
538                for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
539                  if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
540                    ahead += (*_flow)[ie];
541                }
542                if (ahead < 0) ahead = 0;
543
544                // Pushing flow along the edge
545                if (ahead < delta) {
546                  (*_flow)[e] += ahead;
547                  _excess[n] -= ahead;
548                  _excess[t] += ahead;
549                  active_nodes.push_front(t);
550                  hyper[t] = true;
551                  relabel_enabled = false;
552                  break;
553                } else {
554                  (*_flow)[e] += delta;
555                  _excess[n] -= delta;
556                  _excess[t] += delta;
557                  if (_excess[t] > 0 && _excess[t] <= delta)
558                    active_nodes.push_back(t);
559                }
560
561                if (_excess[n] == 0) break;
562              }
563            }
564          }
565
566          if (_excess[n] > 0) {
567            for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
568              if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
569                delta = (*_flow)[e] <= _excess[n] ? (*_flow)[e] : _excess[n];
570                t = _graph.source(e);
571
572                // Push-look-ahead heuristic
573                Capacity ahead = -_excess[t];
574                for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
575                  if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
576                    ahead += _capacity[oe] - (*_flow)[oe];
577                }
578                for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
579                  if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
580                    ahead += (*_flow)[ie];
581                }
582                if (ahead < 0) ahead = 0;
583
584                // Pushing flow along the edge
585                if (ahead < delta) {
586                  (*_flow)[e] -= ahead;
587                  _excess[n] -= ahead;
588                  _excess[t] += ahead;
589                  active_nodes.push_front(t);
590                  hyper[t] = true;
591                  relabel_enabled = false;
592                  break;
593                } else {
594                  (*_flow)[e] -= delta;
595                  _excess[n] -= delta;
596                  _excess[t] += delta;
597                  if (_excess[t] > 0 && _excess[t] <= delta)
598                    active_nodes.push_back(t);
599                }
600
601                if (_excess[n] == 0) break;
602              }
603            }
604          }
605
606          if (relabel_enabled && (_excess[n] > 0 || hyper[n])) {
607            // Performing relabel operation if the node is still active
608            LCost min_red_cost = std::numeric_limits<LCost>::max();
609            for (OutEdgeIt oe(_graph, n); oe != INVALID; ++oe) {
610              if ( _capacity[oe] - (*_flow)[oe] > 0 &&
611                   (*_red_cost)[oe] < min_red_cost )
612                min_red_cost = (*_red_cost)[oe];
613            }
614            for (InEdgeIt ie(_graph, n); ie != INVALID; ++ie) {
615              if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
616                min_red_cost = -(*_red_cost)[ie];
617            }
618            (*_potential)[n] -= min_red_cost + _epsilon;
619            hyper[n] = false;
620          }
621
622          // Removing active nodes with non-positive excess
623          while ( active_nodes.size() > 0 &&
624                  _excess[active_nodes[0]] <= 0 &&
625                  !hyper[active_nodes[0]] ) {
626            active_nodes.pop_front();
627          }
628        }
629      }
630
631      // Computing node potentials for the original costs
632      ResidualCostMap<CostMap> res_cost(_orig_cost);
633      BellmanFord< ResGraph, ResidualCostMap<CostMap> >
634        bf(*_res_graph, res_cost);
635      bf.init(0); bf.start();
636      for (NodeIt n(_graph); n != INVALID; ++n)
637        (*_potential)[n] = bf.dist(n);
638
639      // Handling non-zero lower bounds
640      if (_lower) {
641        for (EdgeIt e(_graph); e != INVALID; ++e)
642          (*_flow)[e] += (*_lower)[e];
643      }
644      return true;
645    }
646
647  }; //class CostScaling
648
649  ///@}
650
651} //namespace lemon
652
653#endif //LEMON_COST_SCALING_H
Note: See TracBrowser for help on using the repository browser.