COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/cost_scaling.h @ 2577:2c6204d4b0f6

Last change on this file since 2577:2c6204d4b0f6 was 2577:2c6204d4b0f6, checked in by Peter Kovacs, 12 years ago

Add a cost scaling min cost flow algorithm.

Add a cost scaling algorithm, which is performing generalized
push-relabel operations. It is almost as efficient as the capacity
scaling algorithm, but slower than network simplex.

File size: 19.0 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  /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.
58  /// - \c CapacityMap::Value and \c SupplyMap::Value must be
59  ///   convertible to each other.
60  /// - All value types must be convertible to \c CostMap::Value, which
61  ///   must be signed type.
62  ///
63  /// \note Edge costs are multiplied with the number of nodes during
64  /// the algorithm so overflow problems may arise more easily than with
65  /// other minimum cost flow algorithms.
66  /// If it is available, <tt>long long int</tt> type is used instead of
67  /// <tt>long int</tt> in the inside computations.
68  ///
69  /// \author Peter Kovacs
70
71  template < typename Graph,
72             typename LowerMap = typename Graph::template EdgeMap<int>,
73             typename CapacityMap = typename Graph::template EdgeMap<int>,
74             typename CostMap = typename Graph::template EdgeMap<int>,
75             typename SupplyMap = typename Graph::template NodeMap<int> >
76  class CostScaling
77  {
78    GRAPH_TYPEDEFS(typename Graph);
79
80    typedef typename CapacityMap::Value Capacity;
81    typedef typename CostMap::Value Cost;
82    typedef typename SupplyMap::Value Supply;
83    typedef typename Graph::template EdgeMap<Capacity> CapacityEdgeMap;
84    typedef typename Graph::template NodeMap<Supply> SupplyNodeMap;
85
86    typedef ResGraphAdaptor< const Graph, Capacity,
87                             CapacityEdgeMap, CapacityEdgeMap > ResGraph;
88    typedef typename ResGraph::Edge ResEdge;
89
90#if defined __GNUC__ && !defined __STRICT_ANSI__
91    typedef long long int LCost;
92#else
93    typedef long int LCost;
94#endif
95    typedef typename Graph::template EdgeMap<LCost> LargeCostMap;
96
97  public:
98
99    /// The type of the flow map.
100    typedef CapacityEdgeMap FlowMap;
101    /// The type of the potential map.
102    typedef typename Graph::template NodeMap<LCost> PotentialMap;
103
104  private:
105
106    /// \brief Map adaptor class for handling residual edge costs.
107    ///
108    /// \ref ResidualCostMap is a map adaptor class for handling
109    /// residual edge costs.
110    class ResidualCostMap : public MapBase<ResEdge, LCost>
111    {
112    private:
113
114      const LargeCostMap &_cost_map;
115
116    public:
117
118      ///\e
119      ResidualCostMap(const LargeCostMap &cost_map) :
120        _cost_map(cost_map) {}
121
122      ///\e
123      LCost operator[](const ResEdge &e) const {
124        return ResGraph::forward(e) ?  _cost_map[e] : -_cost_map[e];
125      }
126
127    }; //class ResidualCostMap
128
129    /// \brief Map adaptor class for handling reduced edge costs.
130    ///
131    /// \ref ReducedCostMap is a map adaptor class for handling reduced
132    /// edge costs.
133    class ReducedCostMap : public MapBase<Edge, LCost>
134    {
135    private:
136
137      const Graph &_gr;
138      const LargeCostMap &_cost_map;
139      const PotentialMap &_pot_map;
140
141    public:
142
143      ///\e
144      ReducedCostMap( const Graph &gr,
145                      const LargeCostMap &cost_map,
146                      const PotentialMap &pot_map ) :
147        _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
148
149      ///\e
150      LCost operator[](const Edge &e) const {
151        return _cost_map[e] + _pot_map[_gr.source(e)]
152                            - _pot_map[_gr.target(e)];
153      }
154
155    }; //class ReducedCostMap
156
157  private:
158
159    // Scaling factor
160    static const int ALPHA = 4;
161
162    // Paramters for heuristics
163    static const int BF_HEURISTIC_EPSILON_BOUND    = 5000;
164    static const int BF_HEURISTIC_BOUND_FACTOR = 3;
165
166  private:
167
168    // The directed graph the algorithm runs on
169    const Graph &_graph;
170    // The original lower bound map
171    const LowerMap *_lower;
172    // The modified capacity map
173    CapacityEdgeMap _capacity;
174    // The original cost map
175    const CostMap &_orig_cost;
176    // The scaled cost map
177    LargeCostMap _cost;
178    // The modified supply map
179    SupplyNodeMap _supply;
180    bool _valid_supply;
181
182    // Edge map of the current flow
183    FlowMap _flow;
184    // Node map of the current potentials
185    PotentialMap _potential;
186
187    // The residual graph
188    ResGraph _res_graph;
189    // The residual cost map
190    ResidualCostMap _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 of the class (with lower bounds).
201    ///
202    /// General constructor of the class (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(graph, 0), _potential(graph, 0),
216      _res_graph(graph, _capacity, _flow), _res_cost(_cost),
217      _red_cost(graph, _cost, _potential), _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 of the class (without lower bounds).
235    ///
236    /// General constructor of the class (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(graph, 0), _potential(graph, 0),
248      _res_graph(graph, _capacity, _flow), _res_cost(_cost),
249      _red_cost(graph, _cost, _potential), _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 of the class (with lower bounds).
258    ///
259    /// Simple constructor of the class (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(graph, 0), _potential(graph, 0),
277      _res_graph(graph, _capacity, _flow), _res_cost(_cost),
278      _red_cost(graph, _cost, _potential), _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 of the class (without lower bounds).
296    ///
297    /// Simple constructor of the class (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(graph, 0), _potential(graph, 0),
313      _res_graph(graph, _capacity, _flow), _res_cost(_cost),
314      _red_cost(graph, _cost, _potential), _excess(graph, 0)
315    {
316      _supply[s] =  flow_value;
317      _supply[t] = -flow_value;
318      _valid_supply = true;
319    }
320
321    /// \brief Runs the algorithm.
322    ///
323    /// Runs the algorithm.
324    ///
325    /// \return \c true if a feasible flow can be found.
326    bool run() {
327      init() && start();
328    }
329
330    /// \brief Returns a const reference to the edge map storing the
331    /// found flow.
332    ///
333    /// Returns a const reference to the edge map storing the found flow.
334    ///
335    /// \pre \ref run() must be called before using this function.
336    const FlowMap& flowMap() const {
337      return _flow;
338    }
339
340    /// \brief Returns a const reference to the node map storing the
341    /// found potentials (the dual solution).
342    ///
343    /// Returns a const reference to the node map storing the found
344    /// potentials (the dual solution).
345    ///
346    /// \pre \ref run() must be called before using this function.
347    const PotentialMap& potentialMap() const {
348      return _potential;
349    }
350
351    /// \brief Returns the total cost of the found flow.
352    ///
353    /// Returns the total cost of the found flow. The complexity of the
354    /// function is \f$ O(e) \f$.
355    ///
356    /// \pre \ref run() must be called before using this function.
357    Cost totalCost() const {
358      Cost c = 0;
359      for (EdgeIt e(_graph); e != INVALID; ++e)
360        c += _flow[e] * _orig_cost[e];
361      return c;
362    }
363
364  private:
365
366    /// Initializes the algorithm.
367    bool init() {
368      if (!_valid_supply) return false;
369
370      // Initializing the scaled cost map and the epsilon parameter
371      Cost max_cost = 0;
372      int node_num = countNodes(_graph);
373      for (EdgeIt e(_graph); e != INVALID; ++e) {
374        _cost[e] = LCost(_orig_cost[e]) * node_num * ALPHA;
375        if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e];
376      }
377      _epsilon = max_cost * node_num;
378
379      // Finding a feasible flow using Circulation
380      Circulation< Graph, ConstMap<Edge, Capacity>, CapacityEdgeMap,
381                   SupplyMap >
382        circulation( _graph, constMap<Edge>((Capacity)0), _capacity,
383                     _supply );
384      return circulation.flowMap(_flow).run();
385    }
386
387
388    /// Executes the algorithm.
389    bool start() {
390      std::deque<Node> active_nodes;
391      typename Graph::template NodeMap<bool> hyper(_graph, false);
392
393      int node_num = countNodes(_graph);
394      for ( ; _epsilon >= 1; _epsilon = _epsilon < ALPHA && _epsilon > 1 ?
395                                        1 : _epsilon / ALPHA )
396      {
397        // Performing price refinement heuristic using Bellman-Ford
398        // algorithm
399        if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
400          typedef ShiftMap<ResidualCostMap> ShiftCostMap;
401          ShiftCostMap shift_cost(_res_cost, _epsilon);
402          BellmanFord<ResGraph, ShiftCostMap> bf(_res_graph, shift_cost);
403          bf.init(0);
404          bool done = false;
405          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
406          for (int i = 0; i < K && !done; ++i)
407            done = bf.processNextWeakRound();
408          if (done) {
409            for (NodeIt n(_graph); n != INVALID; ++n)
410              _potential[n] = bf.dist(n);
411            continue;
412          }
413        }
414
415        // Saturating edges not satisfying the optimality condition
416        Capacity delta;
417        for (EdgeIt e(_graph); e != INVALID; ++e) {
418          if (_capacity[e] - _flow[e] > 0 && _red_cost[e] < 0) {
419            delta = _capacity[e] - _flow[e];
420            _excess[_graph.source(e)] -= delta;
421            _excess[_graph.target(e)] += delta;
422            _flow[e] = _capacity[e];
423          }
424          if (_flow[e] > 0 && -_red_cost[e] < 0) {
425            _excess[_graph.target(e)] -= _flow[e];
426            _excess[_graph.source(e)] += _flow[e];
427            _flow[e] = 0;
428          }
429        }
430
431        // Finding active nodes (i.e. nodes with positive excess)
432        for (NodeIt n(_graph); n != INVALID; ++n)
433          if (_excess[n] > 0) active_nodes.push_back(n);
434
435        // Performing push and relabel operations
436        while (active_nodes.size() > 0) {
437          Node n = active_nodes[0], t;
438          bool relabel_enabled = true;
439
440          // Performing push operations if there are admissible edges
441          if (_excess[n] > 0) {
442            for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
443              if (_capacity[e] - _flow[e] > 0 && _red_cost[e] < 0) {
444                delta = _capacity[e] - _flow[e] <= _excess[n] ?
445                        _capacity[e] - _flow[e] : _excess[n];
446                t = _graph.target(e);
447
448                // Push-look-ahead heuristic
449                Capacity ahead = -_excess[t];
450                for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
451                  if (_capacity[oe] - _flow[oe] > 0 && _red_cost[oe] < 0)
452                    ahead += _capacity[oe] - _flow[oe];
453                }
454                for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
455                  if (_flow[ie] > 0 && -_red_cost[ie] < 0)
456                    ahead += _flow[ie];
457                }
458                if (ahead < 0) ahead = 0;
459
460                // Pushing flow along the edge
461                if (ahead < delta) {
462                  _flow[e] += ahead;
463                  _excess[n] -= ahead;
464                  _excess[t] += ahead;
465                  active_nodes.push_front(t);
466                  hyper[t] = true;
467                  relabel_enabled = false;
468                  break;
469                } else {
470                  _flow[e] += delta;
471                  _excess[n] -= delta;
472                  _excess[t] += delta;
473                  if (_excess[t] > 0 && _excess[t] <= delta)
474                    active_nodes.push_back(t);
475                }
476
477                if (_excess[n] == 0) break;
478              }
479            }
480          }
481
482          if (_excess[n] > 0) {
483            for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
484              if (_flow[e] > 0 && -_red_cost[e] < 0) {
485                delta = _flow[e] <= _excess[n] ? _flow[e] : _excess[n];
486                t = _graph.source(e);
487
488                // Push-look-ahead heuristic
489                Capacity ahead = -_excess[t];
490                for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
491                  if (_capacity[oe] - _flow[oe] > 0 && _red_cost[oe] < 0)
492                    ahead += _capacity[oe] - _flow[oe];
493                }
494                for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
495                  if (_flow[ie] > 0 && -_red_cost[ie] < 0)
496                    ahead += _flow[ie];
497                }
498                if (ahead < 0) ahead = 0;
499
500                // Pushing flow along the edge
501                if (ahead < delta) {
502                  _flow[e] -= ahead;
503                  _excess[n] -= ahead;
504                  _excess[t] += ahead;
505                  active_nodes.push_front(t);
506                  hyper[t] = true;
507                  relabel_enabled = false;
508                  break;
509                } else {
510                  _flow[e] -= delta;
511                  _excess[n] -= delta;
512                  _excess[t] += delta;
513                  if (_excess[t] > 0 && _excess[t] <= delta)
514                    active_nodes.push_back(t);
515                }
516
517                if (_excess[n] == 0) break;
518              }
519            }
520          }
521
522          if (relabel_enabled && (_excess[n] > 0 || hyper[n])) {
523            // Performing relabel operation if the node is still active
524            LCost min_red_cost = std::numeric_limits<LCost>::max();
525            for (OutEdgeIt oe(_graph, n); oe != INVALID; ++oe) {
526              if ( _capacity[oe] - _flow[oe] > 0 &&
527                   _red_cost[oe] < min_red_cost )
528                min_red_cost = _red_cost[oe];
529            }
530            for (InEdgeIt ie(_graph, n); ie != INVALID; ++ie) {
531              if (_flow[ie] > 0 && -_red_cost[ie] < min_red_cost)
532                min_red_cost = -_red_cost[ie];
533            }
534            _potential[n] -= min_red_cost + _epsilon;
535            hyper[n] = false;
536          }
537
538          // Removing active nodes with non-positive excess
539          while ( active_nodes.size() > 0 &&
540                  _excess[active_nodes[0]] <= 0 &&
541                  !hyper[active_nodes[0]] ) {
542            active_nodes.pop_front();
543          }
544        }
545      }
546
547      // Handling non-zero lower bounds
548      if (_lower) {
549        for (EdgeIt e(_graph); e != INVALID; ++e)
550          _flow[e] += (*_lower)[e];
551      }
552      return true;
553    }
554
555  }; //class CostScaling
556
557  ///@}
558
559} //namespace lemon
560
561#endif //LEMON_COST_SCALING_H
Note: See TracBrowser for help on using the repository browser.