COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/cost_scaling.h @ 2623:90defb96ee61

Last change on this file since 2623:90defb96ee61 was 2623:90defb96ee61, checked in by Peter Kovacs, 11 years ago

Add missing pointer initializing in min cost flow classes

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