COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/cost_scaling.h @ 839:f3bc4e9b5f3a

Last change on this file since 839:f3bc4e9b5f3a was 839:f3bc4e9b5f3a, checked in by Peter Kovacs <kpeter@…>, 10 years ago

New heuristics for MCF algorithms (#340)
and some implementation improvements.

  • A useful heuristic is added to NetworkSimplex? to make the initial pivots faster.
  • A powerful global update heuristic is added to CostScaling? and the implementation is reworked with various improvements.
  • Better relabeling in CostScaling? to improve numerical stability and make the code faster.
  • A small improvement is made in CapacityScaling? for better delta computation.
  • Add notes to the classes about the usage of vector<char> instead of vector<bool> for efficiency reasons.
File size: 41.3 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_algs
23/// \file
24/// \brief Cost scaling algorithm for finding a minimum cost flow.
25
26#include <vector>
27#include <deque>
28#include <limits>
29
30#include <lemon/core.h>
31#include <lemon/maps.h>
32#include <lemon/math.h>
33#include <lemon/static_graph.h>
34#include <lemon/circulation.h>
35#include <lemon/bellman_ford.h>
36
37namespace lemon {
38
39  /// \brief Default traits class of CostScaling algorithm.
40  ///
41  /// Default traits class of CostScaling algorithm.
42  /// \tparam GR Digraph type.
43  /// \tparam V The number type used for flow amounts, capacity bounds
44  /// and supply values. By default it is \c int.
45  /// \tparam C The number type used for costs and potentials.
46  /// By default it is the same as \c V.
47#ifdef DOXYGEN
48  template <typename GR, typename V = int, typename C = V>
49#else
50  template < typename GR, typename V = int, typename C = V,
51             bool integer = std::numeric_limits<C>::is_integer >
52#endif
53  struct CostScalingDefaultTraits
54  {
55    /// The type of the digraph
56    typedef GR Digraph;
57    /// The type of the flow amounts, capacity bounds and supply values
58    typedef V Value;
59    /// The type of the arc costs
60    typedef C Cost;
61
62    /// \brief The large cost type used for internal computations
63    ///
64    /// The large cost type used for internal computations.
65    /// It is \c long \c long if the \c Cost type is integer,
66    /// otherwise it is \c double.
67    /// \c Cost must be convertible to \c LargeCost.
68    typedef double LargeCost;
69  };
70
71  // Default traits class for integer cost types
72  template <typename GR, typename V, typename C>
73  struct CostScalingDefaultTraits<GR, V, C, true>
74  {
75    typedef GR Digraph;
76    typedef V Value;
77    typedef C Cost;
78#ifdef LEMON_HAVE_LONG_LONG
79    typedef long long LargeCost;
80#else
81    typedef long LargeCost;
82#endif
83  };
84
85
86  /// \addtogroup min_cost_flow_algs
87  /// @{
88
89  /// \brief Implementation of the Cost Scaling algorithm for
90  /// finding a \ref min_cost_flow "minimum cost flow".
91  ///
92  /// \ref CostScaling implements a cost scaling algorithm that performs
93  /// push/augment and relabel operations for finding a \ref min_cost_flow
94  /// "minimum cost flow" \ref amo93networkflows, \ref goldberg90approximation,
95  /// \ref goldberg97efficient, \ref bunnagel98efficient.
96  /// It is a highly efficient primal-dual solution method, which
97  /// can be viewed as the generalization of the \ref Preflow
98  /// "preflow push-relabel" algorithm for the maximum flow problem.
99  ///
100  /// Most of the parameters of the problem (except for the digraph)
101  /// can be given using separate functions, and the algorithm can be
102  /// executed using the \ref run() function. If some parameters are not
103  /// specified, then default values will be used.
104  ///
105  /// \tparam GR The digraph type the algorithm runs on.
106  /// \tparam V The number type used for flow amounts, capacity bounds
107  /// and supply values in the algorithm. By default it is \c int.
108  /// \tparam C The number type used for costs and potentials in the
109  /// algorithm. By default it is the same as \c V.
110  ///
111  /// \warning Both number types must be signed and all input data must
112  /// be integer.
113  /// \warning This algorithm does not support negative costs for such
114  /// arcs that have infinite upper bound.
115  ///
116  /// \note %CostScaling provides three different internal methods,
117  /// from which the most efficient one is used by default.
118  /// For more information, see \ref Method.
119#ifdef DOXYGEN
120  template <typename GR, typename V, typename C, typename TR>
121#else
122  template < typename GR, typename V = int, typename C = V,
123             typename TR = CostScalingDefaultTraits<GR, V, C> >
124#endif
125  class CostScaling
126  {
127  public:
128
129    /// The type of the digraph
130    typedef typename TR::Digraph Digraph;
131    /// The type of the flow amounts, capacity bounds and supply values
132    typedef typename TR::Value Value;
133    /// The type of the arc costs
134    typedef typename TR::Cost Cost;
135
136    /// \brief The large cost type
137    ///
138    /// The large cost type used for internal computations.
139    /// Using the \ref CostScalingDefaultTraits "default traits class",
140    /// it is \c long \c long if the \c Cost type is integer,
141    /// otherwise it is \c double.
142    typedef typename TR::LargeCost LargeCost;
143
144    /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
145    typedef TR Traits;
146
147  public:
148
149    /// \brief Problem type constants for the \c run() function.
150    ///
151    /// Enum type containing the problem type constants that can be
152    /// returned by the \ref run() function of the algorithm.
153    enum ProblemType {
154      /// The problem has no feasible solution (flow).
155      INFEASIBLE,
156      /// The problem has optimal solution (i.e. it is feasible and
157      /// bounded), and the algorithm has found optimal flow and node
158      /// potentials (primal and dual solutions).
159      OPTIMAL,
160      /// The digraph contains an arc of negative cost and infinite
161      /// upper bound. It means that the objective function is unbounded
162      /// on that arc, however, note that it could actually be bounded
163      /// over the feasible flows, but this algroithm cannot handle
164      /// these cases.
165      UNBOUNDED
166    };
167
168    /// \brief Constants for selecting the internal method.
169    ///
170    /// Enum type containing constants for selecting the internal method
171    /// for the \ref run() function.
172    ///
173    /// \ref CostScaling provides three internal methods that differ mainly
174    /// in their base operations, which are used in conjunction with the
175    /// relabel operation.
176    /// By default, the so called \ref PARTIAL_AUGMENT
177    /// "Partial Augment-Relabel" method is used, which proved to be
178    /// the most efficient and the most robust on various test inputs.
179    /// However, the other methods can be selected using the \ref run()
180    /// function with the proper parameter.
181    enum Method {
182      /// Local push operations are used, i.e. flow is moved only on one
183      /// admissible arc at once.
184      PUSH,
185      /// Augment operations are used, i.e. flow is moved on admissible
186      /// paths from a node with excess to a node with deficit.
187      AUGMENT,
188      /// Partial augment operations are used, i.e. flow is moved on
189      /// admissible paths started from a node with excess, but the
190      /// lengths of these paths are limited. This method can be viewed
191      /// as a combined version of the previous two operations.
192      PARTIAL_AUGMENT
193    };
194
195  private:
196
197    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
198
199    typedef std::vector<int> IntVector;
200    typedef std::vector<Value> ValueVector;
201    typedef std::vector<Cost> CostVector;
202    typedef std::vector<LargeCost> LargeCostVector;
203    typedef std::vector<char> BoolVector;
204    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
205
206  private:
207 
208    template <typename KT, typename VT>
209    class StaticVectorMap {
210    public:
211      typedef KT Key;
212      typedef VT Value;
213     
214      StaticVectorMap(std::vector<Value>& v) : _v(v) {}
215     
216      const Value& operator[](const Key& key) const {
217        return _v[StaticDigraph::id(key)];
218      }
219
220      Value& operator[](const Key& key) {
221        return _v[StaticDigraph::id(key)];
222      }
223     
224      void set(const Key& key, const Value& val) {
225        _v[StaticDigraph::id(key)] = val;
226      }
227
228    private:
229      std::vector<Value>& _v;
230    };
231
232    typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
233    typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
234
235  private:
236
237    // Data related to the underlying digraph
238    const GR &_graph;
239    int _node_num;
240    int _arc_num;
241    int _res_node_num;
242    int _res_arc_num;
243    int _root;
244
245    // Parameters of the problem
246    bool _have_lower;
247    Value _sum_supply;
248    int _sup_node_num;
249
250    // Data structures for storing the digraph
251    IntNodeMap _node_id;
252    IntArcMap _arc_idf;
253    IntArcMap _arc_idb;
254    IntVector _first_out;
255    BoolVector _forward;
256    IntVector _source;
257    IntVector _target;
258    IntVector _reverse;
259
260    // Node and arc data
261    ValueVector _lower;
262    ValueVector _upper;
263    CostVector _scost;
264    ValueVector _supply;
265
266    ValueVector _res_cap;
267    LargeCostVector _cost;
268    LargeCostVector _pi;
269    ValueVector _excess;
270    IntVector _next_out;
271    std::deque<int> _active_nodes;
272
273    // Data for scaling
274    LargeCost _epsilon;
275    int _alpha;
276
277    IntVector _buckets;
278    IntVector _bucket_next;
279    IntVector _bucket_prev;
280    IntVector _rank;
281    int _max_rank;
282 
283    // Data for a StaticDigraph structure
284    typedef std::pair<int, int> IntPair;
285    StaticDigraph _sgr;
286    std::vector<IntPair> _arc_vec;
287    std::vector<LargeCost> _cost_vec;
288    LargeCostArcMap _cost_map;
289    LargeCostNodeMap _pi_map;
290 
291  public:
292 
293    /// \brief Constant for infinite upper bounds (capacities).
294    ///
295    /// Constant for infinite upper bounds (capacities).
296    /// It is \c std::numeric_limits<Value>::infinity() if available,
297    /// \c std::numeric_limits<Value>::max() otherwise.
298    const Value INF;
299
300  public:
301
302    /// \name Named Template Parameters
303    /// @{
304
305    template <typename T>
306    struct SetLargeCostTraits : public Traits {
307      typedef T LargeCost;
308    };
309
310    /// \brief \ref named-templ-param "Named parameter" for setting
311    /// \c LargeCost type.
312    ///
313    /// \ref named-templ-param "Named parameter" for setting \c LargeCost
314    /// type, which is used for internal computations in the algorithm.
315    /// \c Cost must be convertible to \c LargeCost.
316    template <typename T>
317    struct SetLargeCost
318      : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {
319      typedef  CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;
320    };
321
322    /// @}
323
324  public:
325
326    /// \brief Constructor.
327    ///
328    /// The constructor of the class.
329    ///
330    /// \param graph The digraph the algorithm runs on.
331    CostScaling(const GR& graph) :
332      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
333      _cost_map(_cost_vec), _pi_map(_pi),
334      INF(std::numeric_limits<Value>::has_infinity ?
335          std::numeric_limits<Value>::infinity() :
336          std::numeric_limits<Value>::max())
337    {
338      // Check the number types
339      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
340        "The flow type of CostScaling must be signed");
341      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
342        "The cost type of CostScaling must be signed");
343
344      // Resize vectors
345      _node_num = countNodes(_graph);
346      _arc_num = countArcs(_graph);
347      _res_node_num = _node_num + 1;
348      _res_arc_num = 2 * (_arc_num + _node_num);
349      _root = _node_num;
350
351      _first_out.resize(_res_node_num + 1);
352      _forward.resize(_res_arc_num);
353      _source.resize(_res_arc_num);
354      _target.resize(_res_arc_num);
355      _reverse.resize(_res_arc_num);
356
357      _lower.resize(_res_arc_num);
358      _upper.resize(_res_arc_num);
359      _scost.resize(_res_arc_num);
360      _supply.resize(_res_node_num);
361     
362      _res_cap.resize(_res_arc_num);
363      _cost.resize(_res_arc_num);
364      _pi.resize(_res_node_num);
365      _excess.resize(_res_node_num);
366      _next_out.resize(_res_node_num);
367
368      _arc_vec.reserve(_res_arc_num);
369      _cost_vec.reserve(_res_arc_num);
370
371      // Copy the graph
372      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
373      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
374        _node_id[n] = i;
375      }
376      i = 0;
377      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
378        _first_out[i] = j;
379        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
380          _arc_idf[a] = j;
381          _forward[j] = true;
382          _source[j] = i;
383          _target[j] = _node_id[_graph.runningNode(a)];
384        }
385        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
386          _arc_idb[a] = j;
387          _forward[j] = false;
388          _source[j] = i;
389          _target[j] = _node_id[_graph.runningNode(a)];
390        }
391        _forward[j] = false;
392        _source[j] = i;
393        _target[j] = _root;
394        _reverse[j] = k;
395        _forward[k] = true;
396        _source[k] = _root;
397        _target[k] = i;
398        _reverse[k] = j;
399        ++j; ++k;
400      }
401      _first_out[i] = j;
402      _first_out[_res_node_num] = k;
403      for (ArcIt a(_graph); a != INVALID; ++a) {
404        int fi = _arc_idf[a];
405        int bi = _arc_idb[a];
406        _reverse[fi] = bi;
407        _reverse[bi] = fi;
408      }
409     
410      // Reset parameters
411      reset();
412    }
413
414    /// \name Parameters
415    /// The parameters of the algorithm can be specified using these
416    /// functions.
417
418    /// @{
419
420    /// \brief Set the lower bounds on the arcs.
421    ///
422    /// This function sets the lower bounds on the arcs.
423    /// If it is not used before calling \ref run(), the lower bounds
424    /// will be set to zero on all arcs.
425    ///
426    /// \param map An arc map storing the lower bounds.
427    /// Its \c Value type must be convertible to the \c Value type
428    /// of the algorithm.
429    ///
430    /// \return <tt>(*this)</tt>
431    template <typename LowerMap>
432    CostScaling& lowerMap(const LowerMap& map) {
433      _have_lower = true;
434      for (ArcIt a(_graph); a != INVALID; ++a) {
435        _lower[_arc_idf[a]] = map[a];
436        _lower[_arc_idb[a]] = map[a];
437      }
438      return *this;
439    }
440
441    /// \brief Set the upper bounds (capacities) on the arcs.
442    ///
443    /// This function sets the upper bounds (capacities) on the arcs.
444    /// If it is not used before calling \ref run(), the upper bounds
445    /// will be set to \ref INF on all arcs (i.e. the flow value will be
446    /// unbounded from above).
447    ///
448    /// \param map An arc map storing the upper bounds.
449    /// Its \c Value type must be convertible to the \c Value type
450    /// of the algorithm.
451    ///
452    /// \return <tt>(*this)</tt>
453    template<typename UpperMap>
454    CostScaling& upperMap(const UpperMap& map) {
455      for (ArcIt a(_graph); a != INVALID; ++a) {
456        _upper[_arc_idf[a]] = map[a];
457      }
458      return *this;
459    }
460
461    /// \brief Set the costs of the arcs.
462    ///
463    /// This function sets the costs of the arcs.
464    /// If it is not used before calling \ref run(), the costs
465    /// will be set to \c 1 on all arcs.
466    ///
467    /// \param map An arc map storing the costs.
468    /// Its \c Value type must be convertible to the \c Cost type
469    /// of the algorithm.
470    ///
471    /// \return <tt>(*this)</tt>
472    template<typename CostMap>
473    CostScaling& costMap(const CostMap& map) {
474      for (ArcIt a(_graph); a != INVALID; ++a) {
475        _scost[_arc_idf[a]] =  map[a];
476        _scost[_arc_idb[a]] = -map[a];
477      }
478      return *this;
479    }
480
481    /// \brief Set the supply values of the nodes.
482    ///
483    /// This function sets the supply values of the nodes.
484    /// If neither this function nor \ref stSupply() is used before
485    /// calling \ref run(), the supply of each node will be set to zero.
486    ///
487    /// \param map A node map storing the supply values.
488    /// Its \c Value type must be convertible to the \c Value type
489    /// of the algorithm.
490    ///
491    /// \return <tt>(*this)</tt>
492    template<typename SupplyMap>
493    CostScaling& supplyMap(const SupplyMap& map) {
494      for (NodeIt n(_graph); n != INVALID; ++n) {
495        _supply[_node_id[n]] = map[n];
496      }
497      return *this;
498    }
499
500    /// \brief Set single source and target nodes and a supply value.
501    ///
502    /// This function sets a single source node and a single target node
503    /// and the required flow value.
504    /// If neither this function nor \ref supplyMap() is used before
505    /// calling \ref run(), the supply of each node will be set to zero.
506    ///
507    /// Using this function has the same effect as using \ref supplyMap()
508    /// with such a map in which \c k is assigned to \c s, \c -k is
509    /// assigned to \c t and all other nodes have zero supply value.
510    ///
511    /// \param s The source node.
512    /// \param t The target node.
513    /// \param k The required amount of flow from node \c s to node \c t
514    /// (i.e. the supply of \c s and the demand of \c t).
515    ///
516    /// \return <tt>(*this)</tt>
517    CostScaling& stSupply(const Node& s, const Node& t, Value k) {
518      for (int i = 0; i != _res_node_num; ++i) {
519        _supply[i] = 0;
520      }
521      _supply[_node_id[s]] =  k;
522      _supply[_node_id[t]] = -k;
523      return *this;
524    }
525   
526    /// @}
527
528    /// \name Execution control
529    /// The algorithm can be executed using \ref run().
530
531    /// @{
532
533    /// \brief Run the algorithm.
534    ///
535    /// This function runs the algorithm.
536    /// The paramters can be specified using functions \ref lowerMap(),
537    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
538    /// For example,
539    /// \code
540    ///   CostScaling<ListDigraph> cs(graph);
541    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
542    ///     .supplyMap(sup).run();
543    /// \endcode
544    ///
545    /// This function can be called more than once. All the parameters
546    /// that have been given are kept for the next call, unless
547    /// \ref reset() is called, thus only the modified parameters
548    /// have to be set again. See \ref reset() for examples.
549    /// However, the underlying digraph must not be modified after this
550    /// class have been constructed, since it copies and extends the graph.
551    ///
552    /// \param method The internal method that will be used in the
553    /// algorithm. For more information, see \ref Method.
554    /// \param factor The cost scaling factor. It must be larger than one.
555    ///
556    /// \return \c INFEASIBLE if no feasible flow exists,
557    /// \n \c OPTIMAL if the problem has optimal solution
558    /// (i.e. it is feasible and bounded), and the algorithm has found
559    /// optimal flow and node potentials (primal and dual solutions),
560    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
561    /// and infinite upper bound. It means that the objective function
562    /// is unbounded on that arc, however, note that it could actually be
563    /// bounded over the feasible flows, but this algroithm cannot handle
564    /// these cases.
565    ///
566    /// \see ProblemType, Method
567    ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) {
568      _alpha = factor;
569      ProblemType pt = init();
570      if (pt != OPTIMAL) return pt;
571      start(method);
572      return OPTIMAL;
573    }
574
575    /// \brief Reset all the parameters that have been given before.
576    ///
577    /// This function resets all the paramaters that have been given
578    /// before using functions \ref lowerMap(), \ref upperMap(),
579    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
580    ///
581    /// It is useful for multiple run() calls. If this function is not
582    /// used, all the parameters given before are kept for the next
583    /// \ref run() call.
584    /// However, the underlying digraph must not be modified after this
585    /// class have been constructed, since it copies and extends the graph.
586    ///
587    /// For example,
588    /// \code
589    ///   CostScaling<ListDigraph> cs(graph);
590    ///
591    ///   // First run
592    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
593    ///     .supplyMap(sup).run();
594    ///
595    ///   // Run again with modified cost map (reset() is not called,
596    ///   // so only the cost map have to be set again)
597    ///   cost[e] += 100;
598    ///   cs.costMap(cost).run();
599    ///
600    ///   // Run again from scratch using reset()
601    ///   // (the lower bounds will be set to zero on all arcs)
602    ///   cs.reset();
603    ///   cs.upperMap(capacity).costMap(cost)
604    ///     .supplyMap(sup).run();
605    /// \endcode
606    ///
607    /// \return <tt>(*this)</tt>
608    CostScaling& reset() {
609      for (int i = 0; i != _res_node_num; ++i) {
610        _supply[i] = 0;
611      }
612      int limit = _first_out[_root];
613      for (int j = 0; j != limit; ++j) {
614        _lower[j] = 0;
615        _upper[j] = INF;
616        _scost[j] = _forward[j] ? 1 : -1;
617      }
618      for (int j = limit; j != _res_arc_num; ++j) {
619        _lower[j] = 0;
620        _upper[j] = INF;
621        _scost[j] = 0;
622        _scost[_reverse[j]] = 0;
623      }     
624      _have_lower = false;
625      return *this;
626    }
627
628    /// @}
629
630    /// \name Query Functions
631    /// The results of the algorithm can be obtained using these
632    /// functions.\n
633    /// The \ref run() function must be called before using them.
634
635    /// @{
636
637    /// \brief Return the total cost of the found flow.
638    ///
639    /// This function returns the total cost of the found flow.
640    /// Its complexity is O(e).
641    ///
642    /// \note The return type of the function can be specified as a
643    /// template parameter. For example,
644    /// \code
645    ///   cs.totalCost<double>();
646    /// \endcode
647    /// It is useful if the total cost cannot be stored in the \c Cost
648    /// type of the algorithm, which is the default return type of the
649    /// function.
650    ///
651    /// \pre \ref run() must be called before using this function.
652    template <typename Number>
653    Number totalCost() const {
654      Number c = 0;
655      for (ArcIt a(_graph); a != INVALID; ++a) {
656        int i = _arc_idb[a];
657        c += static_cast<Number>(_res_cap[i]) *
658             (-static_cast<Number>(_scost[i]));
659      }
660      return c;
661    }
662
663#ifndef DOXYGEN
664    Cost totalCost() const {
665      return totalCost<Cost>();
666    }
667#endif
668
669    /// \brief Return the flow on the given arc.
670    ///
671    /// This function returns the flow on the given arc.
672    ///
673    /// \pre \ref run() must be called before using this function.
674    Value flow(const Arc& a) const {
675      return _res_cap[_arc_idb[a]];
676    }
677
678    /// \brief Return the flow map (the primal solution).
679    ///
680    /// This function copies the flow value on each arc into the given
681    /// map. The \c Value type of the algorithm must be convertible to
682    /// the \c Value type of the map.
683    ///
684    /// \pre \ref run() must be called before using this function.
685    template <typename FlowMap>
686    void flowMap(FlowMap &map) const {
687      for (ArcIt a(_graph); a != INVALID; ++a) {
688        map.set(a, _res_cap[_arc_idb[a]]);
689      }
690    }
691
692    /// \brief Return the potential (dual value) of the given node.
693    ///
694    /// This function returns the potential (dual value) of the
695    /// given node.
696    ///
697    /// \pre \ref run() must be called before using this function.
698    Cost potential(const Node& n) const {
699      return static_cast<Cost>(_pi[_node_id[n]]);
700    }
701
702    /// \brief Return the potential map (the dual solution).
703    ///
704    /// This function copies the potential (dual value) of each node
705    /// into the given map.
706    /// The \c Cost type of the algorithm must be convertible to the
707    /// \c Value type of the map.
708    ///
709    /// \pre \ref run() must be called before using this function.
710    template <typename PotentialMap>
711    void potentialMap(PotentialMap &map) const {
712      for (NodeIt n(_graph); n != INVALID; ++n) {
713        map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
714      }
715    }
716
717    /// @}
718
719  private:
720
721    // Initialize the algorithm
722    ProblemType init() {
723      if (_res_node_num <= 1) return INFEASIBLE;
724
725      // Check the sum of supply values
726      _sum_supply = 0;
727      for (int i = 0; i != _root; ++i) {
728        _sum_supply += _supply[i];
729      }
730      if (_sum_supply > 0) return INFEASIBLE;
731     
732
733      // Initialize vectors
734      for (int i = 0; i != _res_node_num; ++i) {
735        _pi[i] = 0;
736        _excess[i] = _supply[i];
737      }
738     
739      // Remove infinite upper bounds and check negative arcs
740      const Value MAX = std::numeric_limits<Value>::max();
741      int last_out;
742      if (_have_lower) {
743        for (int i = 0; i != _root; ++i) {
744          last_out = _first_out[i+1];
745          for (int j = _first_out[i]; j != last_out; ++j) {
746            if (_forward[j]) {
747              Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
748              if (c >= MAX) return UNBOUNDED;
749              _excess[i] -= c;
750              _excess[_target[j]] += c;
751            }
752          }
753        }
754      } else {
755        for (int i = 0; i != _root; ++i) {
756          last_out = _first_out[i+1];
757          for (int j = _first_out[i]; j != last_out; ++j) {
758            if (_forward[j] && _scost[j] < 0) {
759              Value c = _upper[j];
760              if (c >= MAX) return UNBOUNDED;
761              _excess[i] -= c;
762              _excess[_target[j]] += c;
763            }
764          }
765        }
766      }
767      Value ex, max_cap = 0;
768      for (int i = 0; i != _res_node_num; ++i) {
769        ex = _excess[i];
770        _excess[i] = 0;
771        if (ex < 0) max_cap -= ex;
772      }
773      for (int j = 0; j != _res_arc_num; ++j) {
774        if (_upper[j] >= MAX) _upper[j] = max_cap;
775      }
776
777      // Initialize the large cost vector and the epsilon parameter
778      _epsilon = 0;
779      LargeCost lc;
780      for (int i = 0; i != _root; ++i) {
781        last_out = _first_out[i+1];
782        for (int j = _first_out[i]; j != last_out; ++j) {
783          lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha;
784          _cost[j] = lc;
785          if (lc > _epsilon) _epsilon = lc;
786        }
787      }
788      _epsilon /= _alpha;
789
790      // Initialize maps for Circulation and remove non-zero lower bounds
791      ConstMap<Arc, Value> low(0);
792      typedef typename Digraph::template ArcMap<Value> ValueArcMap;
793      typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
794      ValueArcMap cap(_graph), flow(_graph);
795      ValueNodeMap sup(_graph);
796      for (NodeIt n(_graph); n != INVALID; ++n) {
797        sup[n] = _supply[_node_id[n]];
798      }
799      if (_have_lower) {
800        for (ArcIt a(_graph); a != INVALID; ++a) {
801          int j = _arc_idf[a];
802          Value c = _lower[j];
803          cap[a] = _upper[j] - c;
804          sup[_graph.source(a)] -= c;
805          sup[_graph.target(a)] += c;
806        }
807      } else {
808        for (ArcIt a(_graph); a != INVALID; ++a) {
809          cap[a] = _upper[_arc_idf[a]];
810        }
811      }
812
813      _sup_node_num = 0;
814      for (NodeIt n(_graph); n != INVALID; ++n) {
815        if (sup[n] > 0) ++_sup_node_num;
816      }
817
818      // Find a feasible flow using Circulation
819      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
820        circ(_graph, low, cap, sup);
821      if (!circ.flowMap(flow).run()) return INFEASIBLE;
822
823      // Set residual capacities and handle GEQ supply type
824      if (_sum_supply < 0) {
825        for (ArcIt a(_graph); a != INVALID; ++a) {
826          Value fa = flow[a];
827          _res_cap[_arc_idf[a]] = cap[a] - fa;
828          _res_cap[_arc_idb[a]] = fa;
829          sup[_graph.source(a)] -= fa;
830          sup[_graph.target(a)] += fa;
831        }
832        for (NodeIt n(_graph); n != INVALID; ++n) {
833          _excess[_node_id[n]] = sup[n];
834        }
835        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
836          int u = _target[a];
837          int ra = _reverse[a];
838          _res_cap[a] = -_sum_supply + 1;
839          _res_cap[ra] = -_excess[u];
840          _cost[a] = 0;
841          _cost[ra] = 0;
842          _excess[u] = 0;
843        }
844      } else {
845        for (ArcIt a(_graph); a != INVALID; ++a) {
846          Value fa = flow[a];
847          _res_cap[_arc_idf[a]] = cap[a] - fa;
848          _res_cap[_arc_idb[a]] = fa;
849        }
850        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
851          int ra = _reverse[a];
852          _res_cap[a] = 0;
853          _res_cap[ra] = 0;
854          _cost[a] = 0;
855          _cost[ra] = 0;
856        }
857      }
858     
859      return OPTIMAL;
860    }
861
862    // Execute the algorithm and transform the results
863    void start(Method method) {
864      // Maximum path length for partial augment
865      const int MAX_PATH_LENGTH = 4;
866
867      // Initialize data structures for buckets     
868      _max_rank = _alpha * _res_node_num;
869      _buckets.resize(_max_rank);
870      _bucket_next.resize(_res_node_num + 1);
871      _bucket_prev.resize(_res_node_num + 1);
872      _rank.resize(_res_node_num + 1);
873 
874      // Execute the algorithm
875      switch (method) {
876        case PUSH:
877          startPush();
878          break;
879        case AUGMENT:
880          startAugment();
881          break;
882        case PARTIAL_AUGMENT:
883          startAugment(MAX_PATH_LENGTH);
884          break;
885      }
886
887      // Compute node potentials for the original costs
888      _arc_vec.clear();
889      _cost_vec.clear();
890      for (int j = 0; j != _res_arc_num; ++j) {
891        if (_res_cap[j] > 0) {
892          _arc_vec.push_back(IntPair(_source[j], _target[j]));
893          _cost_vec.push_back(_scost[j]);
894        }
895      }
896      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
897
898      typename BellmanFord<StaticDigraph, LargeCostArcMap>
899        ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);
900      bf.distMap(_pi_map);
901      bf.init(0);
902      bf.start();
903
904      // Handle non-zero lower bounds
905      if (_have_lower) {
906        int limit = _first_out[_root];
907        for (int j = 0; j != limit; ++j) {
908          if (!_forward[j]) _res_cap[j] += _lower[j];
909        }
910      }
911    }
912   
913    // Initialize a cost scaling phase
914    void initPhase() {
915      // Saturate arcs not satisfying the optimality condition
916      for (int u = 0; u != _res_node_num; ++u) {
917        int last_out = _first_out[u+1];
918        LargeCost pi_u = _pi[u];
919        for (int a = _first_out[u]; a != last_out; ++a) {
920          int v = _target[a];
921          if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) {
922            Value delta = _res_cap[a];
923            _excess[u] -= delta;
924            _excess[v] += delta;
925            _res_cap[a] = 0;
926            _res_cap[_reverse[a]] += delta;
927          }
928        }
929      }
930     
931      // Find active nodes (i.e. nodes with positive excess)
932      for (int u = 0; u != _res_node_num; ++u) {
933        if (_excess[u] > 0) _active_nodes.push_back(u);
934      }
935
936      // Initialize the next arcs
937      for (int u = 0; u != _res_node_num; ++u) {
938        _next_out[u] = _first_out[u];
939      }
940    }
941   
942    // Early termination heuristic
943    bool earlyTermination() {
944      const double EARLY_TERM_FACTOR = 3.0;
945
946      // Build a static residual graph
947      _arc_vec.clear();
948      _cost_vec.clear();
949      for (int j = 0; j != _res_arc_num; ++j) {
950        if (_res_cap[j] > 0) {
951          _arc_vec.push_back(IntPair(_source[j], _target[j]));
952          _cost_vec.push_back(_cost[j] + 1);
953        }
954      }
955      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
956
957      // Run Bellman-Ford algorithm to check if the current flow is optimal
958      BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
959      bf.init(0);
960      bool done = false;
961      int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num)));
962      for (int i = 0; i < K && !done; ++i) {
963        done = bf.processNextWeakRound();
964      }
965      return done;
966    }
967
968    // Global potential update heuristic
969    void globalUpdate() {
970      int bucket_end = _root + 1;
971   
972      // Initialize buckets
973      for (int r = 0; r != _max_rank; ++r) {
974        _buckets[r] = bucket_end;
975      }
976      Value total_excess = 0;
977      for (int i = 0; i != _res_node_num; ++i) {
978        if (_excess[i] < 0) {
979          _rank[i] = 0;
980          _bucket_next[i] = _buckets[0];
981          _bucket_prev[_buckets[0]] = i;
982          _buckets[0] = i;
983        } else {
984          total_excess += _excess[i];
985          _rank[i] = _max_rank;
986        }
987      }
988      if (total_excess == 0) return;
989
990      // Search the buckets
991      int r = 0;
992      for ( ; r != _max_rank; ++r) {
993        while (_buckets[r] != bucket_end) {
994          // Remove the first node from the current bucket
995          int u = _buckets[r];
996          _buckets[r] = _bucket_next[u];
997         
998          // Search the incomming arcs of u
999          LargeCost pi_u = _pi[u];
1000          int last_out = _first_out[u+1];
1001          for (int a = _first_out[u]; a != last_out; ++a) {
1002            int ra = _reverse[a];
1003            if (_res_cap[ra] > 0) {
1004              int v = _source[ra];
1005              int old_rank_v = _rank[v];
1006              if (r < old_rank_v) {
1007                // Compute the new rank of v
1008                LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
1009                int new_rank_v = old_rank_v;
1010                if (nrc < LargeCost(_max_rank))
1011                  new_rank_v = r + 1 + int(nrc);
1012                 
1013                // Change the rank of v
1014                if (new_rank_v < old_rank_v) {
1015                  _rank[v] = new_rank_v;
1016                  _next_out[v] = _first_out[v];
1017                 
1018                  // Remove v from its old bucket
1019                  if (old_rank_v < _max_rank) {
1020                    if (_buckets[old_rank_v] == v) {
1021                      _buckets[old_rank_v] = _bucket_next[v];
1022                    } else {
1023                      _bucket_next[_bucket_prev[v]] = _bucket_next[v];
1024                      _bucket_prev[_bucket_next[v]] = _bucket_prev[v];
1025                    }
1026                  }
1027                 
1028                  // Insert v to its new bucket
1029                  _bucket_next[v] = _buckets[new_rank_v];
1030                  _bucket_prev[_buckets[new_rank_v]] = v;
1031                  _buckets[new_rank_v] = v;
1032                }
1033              }
1034            }
1035          }
1036
1037          // Finish search if there are no more active nodes
1038          if (_excess[u] > 0) {
1039            total_excess -= _excess[u];
1040            if (total_excess <= 0) break;
1041          }
1042        }
1043        if (total_excess <= 0) break;
1044      }
1045     
1046      // Relabel nodes
1047      for (int u = 0; u != _res_node_num; ++u) {
1048        int k = std::min(_rank[u], r);
1049        if (k > 0) {
1050          _pi[u] -= _epsilon * k;
1051          _next_out[u] = _first_out[u];
1052        }
1053      }
1054    }
1055
1056    /// Execute the algorithm performing augment and relabel operations
1057    void startAugment(int max_length = std::numeric_limits<int>::max()) {
1058      // Paramters for heuristics
1059      const int EARLY_TERM_EPSILON_LIMIT = 1000;
1060      const double GLOBAL_UPDATE_FACTOR = 3.0;
1061
1062      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
1063        (_res_node_num + _sup_node_num * _sup_node_num));
1064      int next_update_limit = global_update_freq;
1065     
1066      int relabel_cnt = 0;
1067     
1068      // Perform cost scaling phases
1069      std::vector<int> path;
1070      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1071                                        1 : _epsilon / _alpha )
1072      {
1073        // Early termination heuristic
1074        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
1075          if (earlyTermination()) break;
1076        }
1077       
1078        // Initialize current phase
1079        initPhase();
1080       
1081        // Perform partial augment and relabel operations
1082        while (true) {
1083          // Select an active node (FIFO selection)
1084          while (_active_nodes.size() > 0 &&
1085                 _excess[_active_nodes.front()] <= 0) {
1086            _active_nodes.pop_front();
1087          }
1088          if (_active_nodes.size() == 0) break;
1089          int start = _active_nodes.front();
1090
1091          // Find an augmenting path from the start node
1092          path.clear();
1093          int tip = start;
1094          while (_excess[tip] >= 0 && int(path.size()) < max_length) {
1095            int u;
1096            LargeCost min_red_cost, rc, pi_tip = _pi[tip];
1097            int last_out = _first_out[tip+1];
1098            for (int a = _next_out[tip]; a != last_out; ++a) {
1099              u = _target[a];
1100              if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) {
1101                path.push_back(a);
1102                _next_out[tip] = a;
1103                tip = u;
1104                goto next_step;
1105              }
1106            }
1107
1108            // Relabel tip node
1109            min_red_cost = std::numeric_limits<LargeCost>::max();
1110            if (tip != start) {
1111              int ra = _reverse[path.back()];
1112              min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]];
1113            }
1114            for (int a = _first_out[tip]; a != last_out; ++a) {
1115              rc = _cost[a] + pi_tip - _pi[_target[a]];
1116              if (_res_cap[a] > 0 && rc < min_red_cost) {
1117                min_red_cost = rc;
1118              }
1119            }
1120            _pi[tip] -= min_red_cost + _epsilon;
1121            _next_out[tip] = _first_out[tip];
1122            ++relabel_cnt;
1123
1124            // Step back
1125            if (tip != start) {
1126              tip = _source[path.back()];
1127              path.pop_back();
1128            }
1129
1130          next_step: ;
1131          }
1132
1133          // Augment along the found path (as much flow as possible)
1134          Value delta;
1135          int pa, u, v = start;
1136          for (int i = 0; i != int(path.size()); ++i) {
1137            pa = path[i];
1138            u = v;
1139            v = _target[pa];
1140            delta = std::min(_res_cap[pa], _excess[u]);
1141            _res_cap[pa] -= delta;
1142            _res_cap[_reverse[pa]] += delta;
1143            _excess[u] -= delta;
1144            _excess[v] += delta;
1145            if (_excess[v] > 0 && _excess[v] <= delta)
1146              _active_nodes.push_back(v);
1147          }
1148
1149          // Global update heuristic
1150          if (relabel_cnt >= next_update_limit) {
1151            globalUpdate();
1152            next_update_limit += global_update_freq;
1153          }
1154        }
1155      }
1156    }
1157
1158    /// Execute the algorithm performing push and relabel operations
1159    void startPush() {
1160      // Paramters for heuristics
1161      const int EARLY_TERM_EPSILON_LIMIT = 1000;
1162      const double GLOBAL_UPDATE_FACTOR = 2.0;
1163
1164      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
1165        (_res_node_num + _sup_node_num * _sup_node_num));
1166      int next_update_limit = global_update_freq;
1167
1168      int relabel_cnt = 0;
1169     
1170      // Perform cost scaling phases
1171      BoolVector hyper(_res_node_num, false);
1172      LargeCostVector hyper_cost(_res_node_num);
1173      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1174                                        1 : _epsilon / _alpha )
1175      {
1176        // Early termination heuristic
1177        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
1178          if (earlyTermination()) break;
1179        }
1180       
1181        // Initialize current phase
1182        initPhase();
1183
1184        // Perform push and relabel operations
1185        while (_active_nodes.size() > 0) {
1186          LargeCost min_red_cost, rc, pi_n;
1187          Value delta;
1188          int n, t, a, last_out = _res_arc_num;
1189
1190        next_node:
1191          // Select an active node (FIFO selection)
1192          n = _active_nodes.front();
1193          last_out = _first_out[n+1];
1194          pi_n = _pi[n];
1195         
1196          // Perform push operations if there are admissible arcs
1197          if (_excess[n] > 0) {
1198            for (a = _next_out[n]; a != last_out; ++a) {
1199              if (_res_cap[a] > 0 &&
1200                  _cost[a] + pi_n - _pi[_target[a]] < 0) {
1201                delta = std::min(_res_cap[a], _excess[n]);
1202                t = _target[a];
1203
1204                // Push-look-ahead heuristic
1205                Value ahead = -_excess[t];
1206                int last_out_t = _first_out[t+1];
1207                LargeCost pi_t = _pi[t];
1208                for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
1209                  if (_res_cap[ta] > 0 &&
1210                      _cost[ta] + pi_t - _pi[_target[ta]] < 0)
1211                    ahead += _res_cap[ta];
1212                  if (ahead >= delta) break;
1213                }
1214                if (ahead < 0) ahead = 0;
1215
1216                // Push flow along the arc
1217                if (ahead < delta && !hyper[t]) {
1218                  _res_cap[a] -= ahead;
1219                  _res_cap[_reverse[a]] += ahead;
1220                  _excess[n] -= ahead;
1221                  _excess[t] += ahead;
1222                  _active_nodes.push_front(t);
1223                  hyper[t] = true;
1224                  hyper_cost[t] = _cost[a] + pi_n - pi_t;
1225                  _next_out[n] = a;
1226                  goto next_node;
1227                } else {
1228                  _res_cap[a] -= delta;
1229                  _res_cap[_reverse[a]] += delta;
1230                  _excess[n] -= delta;
1231                  _excess[t] += delta;
1232                  if (_excess[t] > 0 && _excess[t] <= delta)
1233                    _active_nodes.push_back(t);
1234                }
1235
1236                if (_excess[n] == 0) {
1237                  _next_out[n] = a;
1238                  goto remove_nodes;
1239                }
1240              }
1241            }
1242            _next_out[n] = a;
1243          }
1244
1245          // Relabel the node if it is still active (or hyper)
1246          if (_excess[n] > 0 || hyper[n]) {
1247             min_red_cost = hyper[n] ? -hyper_cost[n] :
1248               std::numeric_limits<LargeCost>::max();
1249            for (int a = _first_out[n]; a != last_out; ++a) {
1250              rc = _cost[a] + pi_n - _pi[_target[a]];
1251              if (_res_cap[a] > 0 && rc < min_red_cost) {
1252                min_red_cost = rc;
1253              }
1254            }
1255            _pi[n] -= min_red_cost + _epsilon;
1256            _next_out[n] = _first_out[n];
1257            hyper[n] = false;
1258            ++relabel_cnt;
1259          }
1260       
1261          // Remove nodes that are not active nor hyper
1262        remove_nodes:
1263          while ( _active_nodes.size() > 0 &&
1264                  _excess[_active_nodes.front()] <= 0 &&
1265                  !hyper[_active_nodes.front()] ) {
1266            _active_nodes.pop_front();
1267          }
1268         
1269          // Global update heuristic
1270          if (relabel_cnt >= next_update_limit) {
1271            globalUpdate();
1272            for (int u = 0; u != _res_node_num; ++u)
1273              hyper[u] = false;
1274            next_update_limit += global_update_freq;
1275          }
1276        }
1277      }
1278    }
1279
1280  }; //class CostScaling
1281
1282  ///@}
1283
1284} //namespace lemon
1285
1286#endif //LEMON_COST_SCALING_H
Note: See TracBrowser for help on using the repository browser.