COIN-OR::LEMON - Graph Library

source: lemon/lemon/cost_scaling.h @ 1047:ddd3c0d3d9bf

Last change on this file since 1047:ddd3c0d3d9bf was 1047:ddd3c0d3d9bf, checked in by Peter Kovacs <kpeter@…>, 13 years ago

Implement the scaling Price Refinement heuristic in CostScaling? (#417)
instead of Early Termination.

These two heuristics are similar, but the newer one is faster
and not only makes it possible to skip some epsilon phases, but
it can improve the performance of the other phases, as well.

File size: 49.9 KB
Line 
1/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library.
4 *
5 * Copyright (C) 2003-2010
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  /// In general, \ref NetworkSimplex and \ref CostScaling are the fastest
101  /// implementations available in LEMON for this problem.
102  ///
103  /// Most of the parameters of the problem (except for the digraph)
104  /// can be given using separate functions, and the algorithm can be
105  /// executed using the \ref run() function. If some parameters are not
106  /// specified, then default values will be used.
107  ///
108  /// \tparam GR The digraph type the algorithm runs on.
109  /// \tparam V The number type used for flow amounts, capacity bounds
110  /// and supply values in the algorithm. By default, it is \c int.
111  /// \tparam C The number type used for costs and potentials in the
112  /// algorithm. By default, it is the same as \c V.
113  /// \tparam TR The traits class that defines various types used by the
114  /// algorithm. By default, it is \ref CostScalingDefaultTraits
115  /// "CostScalingDefaultTraits<GR, V, C>".
116  /// In most cases, this parameter should not be set directly,
117  /// consider to use the named template parameters instead.
118  ///
119  /// \warning Both \c V and \c C must be signed number types.
120  /// \warning All input data (capacities, supply values, and costs) must
121  /// be integer.
122  /// \warning This algorithm does not support negative costs for
123  /// arcs having infinite upper bound.
124  ///
125  /// \note %CostScaling provides three different internal methods,
126  /// from which the most efficient one is used by default.
127  /// For more information, see \ref Method.
128#ifdef DOXYGEN
129  template <typename GR, typename V, typename C, typename TR>
130#else
131  template < typename GR, typename V = int, typename C = V,
132             typename TR = CostScalingDefaultTraits<GR, V, C> >
133#endif
134  class CostScaling
135  {
136  public:
137
138    /// The type of the digraph
139    typedef typename TR::Digraph Digraph;
140    /// The type of the flow amounts, capacity bounds and supply values
141    typedef typename TR::Value Value;
142    /// The type of the arc costs
143    typedef typename TR::Cost Cost;
144
145    /// \brief The large cost type
146    ///
147    /// The large cost type used for internal computations.
148    /// By default, it is \c long \c long if the \c Cost type is integer,
149    /// otherwise it is \c double.
150    typedef typename TR::LargeCost LargeCost;
151
152    /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
153    typedef TR Traits;
154
155  public:
156
157    /// \brief Problem type constants for the \c run() function.
158    ///
159    /// Enum type containing the problem type constants that can be
160    /// returned by the \ref run() function of the algorithm.
161    enum ProblemType {
162      /// The problem has no feasible solution (flow).
163      INFEASIBLE,
164      /// The problem has optimal solution (i.e. it is feasible and
165      /// bounded), and the algorithm has found optimal flow and node
166      /// potentials (primal and dual solutions).
167      OPTIMAL,
168      /// The digraph contains an arc of negative cost and infinite
169      /// upper bound. It means that the objective function is unbounded
170      /// on that arc, however, note that it could actually be bounded
171      /// over the feasible flows, but this algroithm cannot handle
172      /// these cases.
173      UNBOUNDED
174    };
175
176    /// \brief Constants for selecting the internal method.
177    ///
178    /// Enum type containing constants for selecting the internal method
179    /// for the \ref run() function.
180    ///
181    /// \ref CostScaling provides three internal methods that differ mainly
182    /// in their base operations, which are used in conjunction with the
183    /// relabel operation.
184    /// By default, the so called \ref PARTIAL_AUGMENT
185    /// "Partial Augment-Relabel" method is used, which turned out to be
186    /// the most efficient and the most robust on various test inputs.
187    /// However, the other methods can be selected using the \ref run()
188    /// function with the proper parameter.
189    enum Method {
190      /// Local push operations are used, i.e. flow is moved only on one
191      /// admissible arc at once.
192      PUSH,
193      /// Augment operations are used, i.e. flow is moved on admissible
194      /// paths from a node with excess to a node with deficit.
195      AUGMENT,
196      /// Partial augment operations are used, i.e. flow is moved on
197      /// admissible paths started from a node with excess, but the
198      /// lengths of these paths are limited. This method can be viewed
199      /// as a combined version of the previous two operations.
200      PARTIAL_AUGMENT
201    };
202
203  private:
204
205    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
206
207    typedef std::vector<int> IntVector;
208    typedef std::vector<Value> ValueVector;
209    typedef std::vector<Cost> CostVector;
210    typedef std::vector<LargeCost> LargeCostVector;
211    typedef std::vector<char> BoolVector;
212    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
213
214  private:
215
216    template <typename KT, typename VT>
217    class StaticVectorMap {
218    public:
219      typedef KT Key;
220      typedef VT Value;
221
222      StaticVectorMap(std::vector<Value>& v) : _v(v) {}
223
224      const Value& operator[](const Key& key) const {
225        return _v[StaticDigraph::id(key)];
226      }
227
228      Value& operator[](const Key& key) {
229        return _v[StaticDigraph::id(key)];
230      }
231
232      void set(const Key& key, const Value& val) {
233        _v[StaticDigraph::id(key)] = val;
234      }
235
236    private:
237      std::vector<Value>& _v;
238    };
239
240    typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
241    typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
242
243  private:
244
245    // Data related to the underlying digraph
246    const GR &_graph;
247    int _node_num;
248    int _arc_num;
249    int _res_node_num;
250    int _res_arc_num;
251    int _root;
252
253    // Parameters of the problem
254    bool _have_lower;
255    Value _sum_supply;
256    int _sup_node_num;
257
258    // Data structures for storing the digraph
259    IntNodeMap _node_id;
260    IntArcMap _arc_idf;
261    IntArcMap _arc_idb;
262    IntVector _first_out;
263    BoolVector _forward;
264    IntVector _source;
265    IntVector _target;
266    IntVector _reverse;
267
268    // Node and arc data
269    ValueVector _lower;
270    ValueVector _upper;
271    CostVector _scost;
272    ValueVector _supply;
273
274    ValueVector _res_cap;
275    LargeCostVector _cost;
276    LargeCostVector _pi;
277    ValueVector _excess;
278    IntVector _next_out;
279    std::deque<int> _active_nodes;
280
281    // Data for scaling
282    LargeCost _epsilon;
283    int _alpha;
284
285    IntVector _buckets;
286    IntVector _bucket_next;
287    IntVector _bucket_prev;
288    IntVector _rank;
289    int _max_rank;
290
291    // Data for a StaticDigraph structure
292    typedef std::pair<int, int> IntPair;
293    StaticDigraph _sgr;
294    std::vector<IntPair> _arc_vec;
295    std::vector<LargeCost> _cost_vec;
296    LargeCostArcMap _cost_map;
297    LargeCostNodeMap _pi_map;
298
299  public:
300
301    /// \brief Constant for infinite upper bounds (capacities).
302    ///
303    /// Constant for infinite upper bounds (capacities).
304    /// It is \c std::numeric_limits<Value>::infinity() if available,
305    /// \c std::numeric_limits<Value>::max() otherwise.
306    const Value INF;
307
308  public:
309
310    /// \name Named Template Parameters
311    /// @{
312
313    template <typename T>
314    struct SetLargeCostTraits : public Traits {
315      typedef T LargeCost;
316    };
317
318    /// \brief \ref named-templ-param "Named parameter" for setting
319    /// \c LargeCost type.
320    ///
321    /// \ref named-templ-param "Named parameter" for setting \c LargeCost
322    /// type, which is used for internal computations in the algorithm.
323    /// \c Cost must be convertible to \c LargeCost.
324    template <typename T>
325    struct SetLargeCost
326      : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {
327      typedef  CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;
328    };
329
330    /// @}
331
332  protected:
333
334    CostScaling() {}
335
336  public:
337
338    /// \brief Constructor.
339    ///
340    /// The constructor of the class.
341    ///
342    /// \param graph The digraph the algorithm runs on.
343    CostScaling(const GR& graph) :
344      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
345      _cost_map(_cost_vec), _pi_map(_pi),
346      INF(std::numeric_limits<Value>::has_infinity ?
347          std::numeric_limits<Value>::infinity() :
348          std::numeric_limits<Value>::max())
349    {
350      // Check the number types
351      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
352        "The flow type of CostScaling must be signed");
353      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
354        "The cost type of CostScaling must be signed");
355
356      // Reset data structures
357      reset();
358    }
359
360    /// \name Parameters
361    /// The parameters of the algorithm can be specified using these
362    /// functions.
363
364    /// @{
365
366    /// \brief Set the lower bounds on the arcs.
367    ///
368    /// This function sets the lower bounds on the arcs.
369    /// If it is not used before calling \ref run(), the lower bounds
370    /// will be set to zero on all arcs.
371    ///
372    /// \param map An arc map storing the lower bounds.
373    /// Its \c Value type must be convertible to the \c Value type
374    /// of the algorithm.
375    ///
376    /// \return <tt>(*this)</tt>
377    template <typename LowerMap>
378    CostScaling& lowerMap(const LowerMap& map) {
379      _have_lower = true;
380      for (ArcIt a(_graph); a != INVALID; ++a) {
381        _lower[_arc_idf[a]] = map[a];
382        _lower[_arc_idb[a]] = map[a];
383      }
384      return *this;
385    }
386
387    /// \brief Set the upper bounds (capacities) on the arcs.
388    ///
389    /// This function sets the upper bounds (capacities) on the arcs.
390    /// If it is not used before calling \ref run(), the upper bounds
391    /// will be set to \ref INF on all arcs (i.e. the flow value will be
392    /// unbounded from above).
393    ///
394    /// \param map An arc map storing the upper bounds.
395    /// Its \c Value type must be convertible to the \c Value type
396    /// of the algorithm.
397    ///
398    /// \return <tt>(*this)</tt>
399    template<typename UpperMap>
400    CostScaling& upperMap(const UpperMap& map) {
401      for (ArcIt a(_graph); a != INVALID; ++a) {
402        _upper[_arc_idf[a]] = map[a];
403      }
404      return *this;
405    }
406
407    /// \brief Set the costs of the arcs.
408    ///
409    /// This function sets the costs of the arcs.
410    /// If it is not used before calling \ref run(), the costs
411    /// will be set to \c 1 on all arcs.
412    ///
413    /// \param map An arc map storing the costs.
414    /// Its \c Value type must be convertible to the \c Cost type
415    /// of the algorithm.
416    ///
417    /// \return <tt>(*this)</tt>
418    template<typename CostMap>
419    CostScaling& costMap(const CostMap& map) {
420      for (ArcIt a(_graph); a != INVALID; ++a) {
421        _scost[_arc_idf[a]] =  map[a];
422        _scost[_arc_idb[a]] = -map[a];
423      }
424      return *this;
425    }
426
427    /// \brief Set the supply values of the nodes.
428    ///
429    /// This function sets the supply values of the nodes.
430    /// If neither this function nor \ref stSupply() is used before
431    /// calling \ref run(), the supply of each node will be set to zero.
432    ///
433    /// \param map A node map storing the supply values.
434    /// Its \c Value type must be convertible to the \c Value type
435    /// of the algorithm.
436    ///
437    /// \return <tt>(*this)</tt>
438    template<typename SupplyMap>
439    CostScaling& supplyMap(const SupplyMap& map) {
440      for (NodeIt n(_graph); n != INVALID; ++n) {
441        _supply[_node_id[n]] = map[n];
442      }
443      return *this;
444    }
445
446    /// \brief Set single source and target nodes and a supply value.
447    ///
448    /// This function sets a single source node and a single target node
449    /// and the required flow value.
450    /// If neither this function nor \ref supplyMap() is used before
451    /// calling \ref run(), the supply of each node will be set to zero.
452    ///
453    /// Using this function has the same effect as using \ref supplyMap()
454    /// with a map in which \c k is assigned to \c s, \c -k is
455    /// assigned to \c t and all other nodes have zero supply value.
456    ///
457    /// \param s The source node.
458    /// \param t The target node.
459    /// \param k The required amount of flow from node \c s to node \c t
460    /// (i.e. the supply of \c s and the demand of \c t).
461    ///
462    /// \return <tt>(*this)</tt>
463    CostScaling& stSupply(const Node& s, const Node& t, Value k) {
464      for (int i = 0; i != _res_node_num; ++i) {
465        _supply[i] = 0;
466      }
467      _supply[_node_id[s]] =  k;
468      _supply[_node_id[t]] = -k;
469      return *this;
470    }
471
472    /// @}
473
474    /// \name Execution control
475    /// The algorithm can be executed using \ref run().
476
477    /// @{
478
479    /// \brief Run the algorithm.
480    ///
481    /// This function runs the algorithm.
482    /// The paramters can be specified using functions \ref lowerMap(),
483    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
484    /// For example,
485    /// \code
486    ///   CostScaling<ListDigraph> cs(graph);
487    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
488    ///     .supplyMap(sup).run();
489    /// \endcode
490    ///
491    /// This function can be called more than once. All the given parameters
492    /// are kept for the next call, unless \ref resetParams() or \ref reset()
493    /// is used, thus only the modified parameters have to be set again.
494    /// If the underlying digraph was also modified after the construction
495    /// of the class (or the last \ref reset() call), then the \ref reset()
496    /// function must be called.
497    ///
498    /// \param method The internal method that will be used in the
499    /// algorithm. For more information, see \ref Method.
500    /// \param factor The cost scaling factor. It must be larger than one.
501    ///
502    /// \return \c INFEASIBLE if no feasible flow exists,
503    /// \n \c OPTIMAL if the problem has optimal solution
504    /// (i.e. it is feasible and bounded), and the algorithm has found
505    /// optimal flow and node potentials (primal and dual solutions),
506    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
507    /// and infinite upper bound. It means that the objective function
508    /// is unbounded on that arc, however, note that it could actually be
509    /// bounded over the feasible flows, but this algroithm cannot handle
510    /// these cases.
511    ///
512    /// \see ProblemType, Method
513    /// \see resetParams(), reset()
514    ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) {
515      _alpha = factor;
516      ProblemType pt = init();
517      if (pt != OPTIMAL) return pt;
518      start(method);
519      return OPTIMAL;
520    }
521
522    /// \brief Reset all the parameters that have been given before.
523    ///
524    /// This function resets all the paramaters that have been given
525    /// before using functions \ref lowerMap(), \ref upperMap(),
526    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
527    ///
528    /// It is useful for multiple \ref run() calls. Basically, all the given
529    /// parameters are kept for the next \ref run() call, unless
530    /// \ref resetParams() or \ref reset() is used.
531    /// If the underlying digraph was also modified after the construction
532    /// of the class or the last \ref reset() call, then the \ref reset()
533    /// function must be used, otherwise \ref resetParams() is sufficient.
534    ///
535    /// For example,
536    /// \code
537    ///   CostScaling<ListDigraph> cs(graph);
538    ///
539    ///   // First run
540    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
541    ///     .supplyMap(sup).run();
542    ///
543    ///   // Run again with modified cost map (resetParams() is not called,
544    ///   // so only the cost map have to be set again)
545    ///   cost[e] += 100;
546    ///   cs.costMap(cost).run();
547    ///
548    ///   // Run again from scratch using resetParams()
549    ///   // (the lower bounds will be set to zero on all arcs)
550    ///   cs.resetParams();
551    ///   cs.upperMap(capacity).costMap(cost)
552    ///     .supplyMap(sup).run();
553    /// \endcode
554    ///
555    /// \return <tt>(*this)</tt>
556    ///
557    /// \see reset(), run()
558    CostScaling& resetParams() {
559      for (int i = 0; i != _res_node_num; ++i) {
560        _supply[i] = 0;
561      }
562      int limit = _first_out[_root];
563      for (int j = 0; j != limit; ++j) {
564        _lower[j] = 0;
565        _upper[j] = INF;
566        _scost[j] = _forward[j] ? 1 : -1;
567      }
568      for (int j = limit; j != _res_arc_num; ++j) {
569        _lower[j] = 0;
570        _upper[j] = INF;
571        _scost[j] = 0;
572        _scost[_reverse[j]] = 0;
573      }
574      _have_lower = false;
575      return *this;
576    }
577
578    /// \brief Reset the internal data structures and all the parameters
579    /// that have been given before.
580    ///
581    /// This function resets the internal data structures and all the
582    /// paramaters that have been given before using functions \ref lowerMap(),
583    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
584    ///
585    /// It is useful for multiple \ref run() calls. By default, all the given
586    /// parameters are kept for the next \ref run() call, unless
587    /// \ref resetParams() or \ref reset() is used.
588    /// If the underlying digraph was also modified after the construction
589    /// of the class or the last \ref reset() call, then the \ref reset()
590    /// function must be used, otherwise \ref resetParams() is sufficient.
591    ///
592    /// See \ref resetParams() for examples.
593    ///
594    /// \return <tt>(*this)</tt>
595    ///
596    /// \see resetParams(), run()
597    CostScaling& reset() {
598      // Resize vectors
599      _node_num = countNodes(_graph);
600      _arc_num = countArcs(_graph);
601      _res_node_num = _node_num + 1;
602      _res_arc_num = 2 * (_arc_num + _node_num);
603      _root = _node_num;
604
605      _first_out.resize(_res_node_num + 1);
606      _forward.resize(_res_arc_num);
607      _source.resize(_res_arc_num);
608      _target.resize(_res_arc_num);
609      _reverse.resize(_res_arc_num);
610
611      _lower.resize(_res_arc_num);
612      _upper.resize(_res_arc_num);
613      _scost.resize(_res_arc_num);
614      _supply.resize(_res_node_num);
615
616      _res_cap.resize(_res_arc_num);
617      _cost.resize(_res_arc_num);
618      _pi.resize(_res_node_num);
619      _excess.resize(_res_node_num);
620      _next_out.resize(_res_node_num);
621
622      _arc_vec.reserve(_res_arc_num);
623      _cost_vec.reserve(_res_arc_num);
624
625      // Copy the graph
626      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
627      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
628        _node_id[n] = i;
629      }
630      i = 0;
631      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
632        _first_out[i] = j;
633        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
634          _arc_idf[a] = j;
635          _forward[j] = true;
636          _source[j] = i;
637          _target[j] = _node_id[_graph.runningNode(a)];
638        }
639        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
640          _arc_idb[a] = j;
641          _forward[j] = false;
642          _source[j] = i;
643          _target[j] = _node_id[_graph.runningNode(a)];
644        }
645        _forward[j] = false;
646        _source[j] = i;
647        _target[j] = _root;
648        _reverse[j] = k;
649        _forward[k] = true;
650        _source[k] = _root;
651        _target[k] = i;
652        _reverse[k] = j;
653        ++j; ++k;
654      }
655      _first_out[i] = j;
656      _first_out[_res_node_num] = k;
657      for (ArcIt a(_graph); a != INVALID; ++a) {
658        int fi = _arc_idf[a];
659        int bi = _arc_idb[a];
660        _reverse[fi] = bi;
661        _reverse[bi] = fi;
662      }
663
664      // Reset parameters
665      resetParams();
666      return *this;
667    }
668
669    /// @}
670
671    /// \name Query Functions
672    /// The results of the algorithm can be obtained using these
673    /// functions.\n
674    /// The \ref run() function must be called before using them.
675
676    /// @{
677
678    /// \brief Return the total cost of the found flow.
679    ///
680    /// This function returns the total cost of the found flow.
681    /// Its complexity is O(e).
682    ///
683    /// \note The return type of the function can be specified as a
684    /// template parameter. For example,
685    /// \code
686    ///   cs.totalCost<double>();
687    /// \endcode
688    /// It is useful if the total cost cannot be stored in the \c Cost
689    /// type of the algorithm, which is the default return type of the
690    /// function.
691    ///
692    /// \pre \ref run() must be called before using this function.
693    template <typename Number>
694    Number totalCost() const {
695      Number c = 0;
696      for (ArcIt a(_graph); a != INVALID; ++a) {
697        int i = _arc_idb[a];
698        c += static_cast<Number>(_res_cap[i]) *
699             (-static_cast<Number>(_scost[i]));
700      }
701      return c;
702    }
703
704#ifndef DOXYGEN
705    Cost totalCost() const {
706      return totalCost<Cost>();
707    }
708#endif
709
710    /// \brief Return the flow on the given arc.
711    ///
712    /// This function returns the flow on the given arc.
713    ///
714    /// \pre \ref run() must be called before using this function.
715    Value flow(const Arc& a) const {
716      return _res_cap[_arc_idb[a]];
717    }
718
719    /// \brief Return the flow map (the primal solution).
720    ///
721    /// This function copies the flow value on each arc into the given
722    /// map. The \c Value type of the algorithm must be convertible to
723    /// the \c Value type of the map.
724    ///
725    /// \pre \ref run() must be called before using this function.
726    template <typename FlowMap>
727    void flowMap(FlowMap &map) const {
728      for (ArcIt a(_graph); a != INVALID; ++a) {
729        map.set(a, _res_cap[_arc_idb[a]]);
730      }
731    }
732
733    /// \brief Return the potential (dual value) of the given node.
734    ///
735    /// This function returns the potential (dual value) of the
736    /// given node.
737    ///
738    /// \pre \ref run() must be called before using this function.
739    Cost potential(const Node& n) const {
740      return static_cast<Cost>(_pi[_node_id[n]]);
741    }
742
743    /// \brief Return the potential map (the dual solution).
744    ///
745    /// This function copies the potential (dual value) of each node
746    /// into the given map.
747    /// The \c Cost type of the algorithm must be convertible to the
748    /// \c Value type of the map.
749    ///
750    /// \pre \ref run() must be called before using this function.
751    template <typename PotentialMap>
752    void potentialMap(PotentialMap &map) const {
753      for (NodeIt n(_graph); n != INVALID; ++n) {
754        map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
755      }
756    }
757
758    /// @}
759
760  private:
761
762    // Initialize the algorithm
763    ProblemType init() {
764      if (_res_node_num <= 1) return INFEASIBLE;
765
766      // Check the sum of supply values
767      _sum_supply = 0;
768      for (int i = 0; i != _root; ++i) {
769        _sum_supply += _supply[i];
770      }
771      if (_sum_supply > 0) return INFEASIBLE;
772
773
774      // Initialize vectors
775      for (int i = 0; i != _res_node_num; ++i) {
776        _pi[i] = 0;
777        _excess[i] = _supply[i];
778      }
779
780      // Remove infinite upper bounds and check negative arcs
781      const Value MAX = std::numeric_limits<Value>::max();
782      int last_out;
783      if (_have_lower) {
784        for (int i = 0; i != _root; ++i) {
785          last_out = _first_out[i+1];
786          for (int j = _first_out[i]; j != last_out; ++j) {
787            if (_forward[j]) {
788              Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
789              if (c >= MAX) return UNBOUNDED;
790              _excess[i] -= c;
791              _excess[_target[j]] += c;
792            }
793          }
794        }
795      } else {
796        for (int i = 0; i != _root; ++i) {
797          last_out = _first_out[i+1];
798          for (int j = _first_out[i]; j != last_out; ++j) {
799            if (_forward[j] && _scost[j] < 0) {
800              Value c = _upper[j];
801              if (c >= MAX) return UNBOUNDED;
802              _excess[i] -= c;
803              _excess[_target[j]] += c;
804            }
805          }
806        }
807      }
808      Value ex, max_cap = 0;
809      for (int i = 0; i != _res_node_num; ++i) {
810        ex = _excess[i];
811        _excess[i] = 0;
812        if (ex < 0) max_cap -= ex;
813      }
814      for (int j = 0; j != _res_arc_num; ++j) {
815        if (_upper[j] >= MAX) _upper[j] = max_cap;
816      }
817
818      // Initialize the large cost vector and the epsilon parameter
819      _epsilon = 0;
820      LargeCost lc;
821      for (int i = 0; i != _root; ++i) {
822        last_out = _first_out[i+1];
823        for (int j = _first_out[i]; j != last_out; ++j) {
824          lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha;
825          _cost[j] = lc;
826          if (lc > _epsilon) _epsilon = lc;
827        }
828      }
829      _epsilon /= _alpha;
830
831      // Initialize maps for Circulation and remove non-zero lower bounds
832      ConstMap<Arc, Value> low(0);
833      typedef typename Digraph::template ArcMap<Value> ValueArcMap;
834      typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
835      ValueArcMap cap(_graph), flow(_graph);
836      ValueNodeMap sup(_graph);
837      for (NodeIt n(_graph); n != INVALID; ++n) {
838        sup[n] = _supply[_node_id[n]];
839      }
840      if (_have_lower) {
841        for (ArcIt a(_graph); a != INVALID; ++a) {
842          int j = _arc_idf[a];
843          Value c = _lower[j];
844          cap[a] = _upper[j] - c;
845          sup[_graph.source(a)] -= c;
846          sup[_graph.target(a)] += c;
847        }
848      } else {
849        for (ArcIt a(_graph); a != INVALID; ++a) {
850          cap[a] = _upper[_arc_idf[a]];
851        }
852      }
853
854      _sup_node_num = 0;
855      for (NodeIt n(_graph); n != INVALID; ++n) {
856        if (sup[n] > 0) ++_sup_node_num;
857      }
858
859      // Find a feasible flow using Circulation
860      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
861        circ(_graph, low, cap, sup);
862      if (!circ.flowMap(flow).run()) return INFEASIBLE;
863
864      // Set residual capacities and handle GEQ supply type
865      if (_sum_supply < 0) {
866        for (ArcIt a(_graph); a != INVALID; ++a) {
867          Value fa = flow[a];
868          _res_cap[_arc_idf[a]] = cap[a] - fa;
869          _res_cap[_arc_idb[a]] = fa;
870          sup[_graph.source(a)] -= fa;
871          sup[_graph.target(a)] += fa;
872        }
873        for (NodeIt n(_graph); n != INVALID; ++n) {
874          _excess[_node_id[n]] = sup[n];
875        }
876        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
877          int u = _target[a];
878          int ra = _reverse[a];
879          _res_cap[a] = -_sum_supply + 1;
880          _res_cap[ra] = -_excess[u];
881          _cost[a] = 0;
882          _cost[ra] = 0;
883          _excess[u] = 0;
884        }
885      } else {
886        for (ArcIt a(_graph); a != INVALID; ++a) {
887          Value fa = flow[a];
888          _res_cap[_arc_idf[a]] = cap[a] - fa;
889          _res_cap[_arc_idb[a]] = fa;
890        }
891        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
892          int ra = _reverse[a];
893          _res_cap[a] = 0;
894          _res_cap[ra] = 0;
895          _cost[a] = 0;
896          _cost[ra] = 0;
897        }
898      }
899
900      // Initialize data structures for buckets
901      _max_rank = _alpha * _res_node_num;
902      _buckets.resize(_max_rank);
903      _bucket_next.resize(_res_node_num + 1);
904      _bucket_prev.resize(_res_node_num + 1);
905      _rank.resize(_res_node_num + 1);
906
907      return OPTIMAL;
908    }
909
910    // Execute the algorithm and transform the results
911    void start(Method method) {
912      const int MAX_PARTIAL_PATH_LENGTH = 4;
913
914      switch (method) {
915        case PUSH:
916          startPush();
917          break;
918        case AUGMENT:
919          startAugment(_res_node_num - 1);
920          break;
921        case PARTIAL_AUGMENT:
922          startAugment(MAX_PARTIAL_PATH_LENGTH);
923          break;
924      }
925
926      // Compute node potentials for the original costs
927      _arc_vec.clear();
928      _cost_vec.clear();
929      for (int j = 0; j != _res_arc_num; ++j) {
930        if (_res_cap[j] > 0) {
931          _arc_vec.push_back(IntPair(_source[j], _target[j]));
932          _cost_vec.push_back(_scost[j]);
933        }
934      }
935      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
936
937      typename BellmanFord<StaticDigraph, LargeCostArcMap>
938        ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);
939      bf.distMap(_pi_map);
940      bf.init(0);
941      bf.start();
942
943      // Handle non-zero lower bounds
944      if (_have_lower) {
945        int limit = _first_out[_root];
946        for (int j = 0; j != limit; ++j) {
947          if (!_forward[j]) _res_cap[j] += _lower[j];
948        }
949      }
950    }
951
952    // Initialize a cost scaling phase
953    void initPhase() {
954      // Saturate arcs not satisfying the optimality condition
955      for (int u = 0; u != _res_node_num; ++u) {
956        int last_out = _first_out[u+1];
957        LargeCost pi_u = _pi[u];
958        for (int a = _first_out[u]; a != last_out; ++a) {
959          Value delta = _res_cap[a];
960          if (delta > 0) {
961            int v = _target[a];
962            if (_cost[a] + pi_u - _pi[v] < 0) {
963              _excess[u] -= delta;
964              _excess[v] += delta;
965              _res_cap[a] = 0;
966              _res_cap[_reverse[a]] += delta;
967            }
968          }
969        }
970      }
971
972      // Find active nodes (i.e. nodes with positive excess)
973      for (int u = 0; u != _res_node_num; ++u) {
974        if (_excess[u] > 0) _active_nodes.push_back(u);
975      }
976
977      // Initialize the next arcs
978      for (int u = 0; u != _res_node_num; ++u) {
979        _next_out[u] = _first_out[u];
980      }
981    }
982
983    // Price (potential) refinement heuristic
984    bool priceRefinement() {
985
986      // Stack for stroing the topological order
987      IntVector stack(_res_node_num);
988      int stack_top;
989
990      // Perform phases
991      while (topologicalSort(stack, stack_top)) {
992
993        // Compute node ranks in the acyclic admissible network and
994        // store the nodes in buckets
995        for (int i = 0; i != _res_node_num; ++i) {
996          _rank[i] = 0;
997        }
998        const int bucket_end = _root + 1;
999        for (int r = 0; r != _max_rank; ++r) {
1000          _buckets[r] = bucket_end;
1001        }
1002        int top_rank = 0;
1003        for ( ; stack_top >= 0; --stack_top) {
1004          int u = stack[stack_top], v;
1005          int rank_u = _rank[u];
1006
1007          LargeCost rc, pi_u = _pi[u];
1008          int last_out = _first_out[u+1];
1009          for (int a = _first_out[u]; a != last_out; ++a) {
1010            if (_res_cap[a] > 0) {
1011              v = _target[a];
1012              rc = _cost[a] + pi_u - _pi[v];
1013              if (rc < 0) {
1014                LargeCost nrc = static_cast<LargeCost>((-rc - 0.5) / _epsilon);
1015                if (nrc < LargeCost(_max_rank)) {
1016                  int new_rank_v = rank_u + static_cast<int>(nrc);
1017                  if (new_rank_v > _rank[v]) {
1018                    _rank[v] = new_rank_v;
1019                  }
1020                }
1021              }
1022            }
1023          }
1024
1025          if (rank_u > 0) {
1026            top_rank = std::max(top_rank, rank_u);
1027            int bfirst = _buckets[rank_u];
1028            _bucket_next[u] = bfirst;
1029            _bucket_prev[bfirst] = u;
1030            _buckets[rank_u] = u;
1031          }
1032        }
1033
1034        // Check if the current flow is epsilon-optimal
1035        if (top_rank == 0) {
1036          return true;
1037        }
1038
1039        // Process buckets in top-down order
1040        for (int rank = top_rank; rank > 0; --rank) {
1041          while (_buckets[rank] != bucket_end) {
1042            // Remove the first node from the current bucket
1043            int u = _buckets[rank];
1044            _buckets[rank] = _bucket_next[u];
1045
1046            // Search the outgoing arcs of u
1047            LargeCost rc, pi_u = _pi[u];
1048            int last_out = _first_out[u+1];
1049            int v, old_rank_v, new_rank_v;
1050            for (int a = _first_out[u]; a != last_out; ++a) {
1051              if (_res_cap[a] > 0) {
1052                v = _target[a];
1053                old_rank_v = _rank[v];
1054
1055                if (old_rank_v < rank) {
1056
1057                  // Compute the new rank of node v
1058                  rc = _cost[a] + pi_u - _pi[v];
1059                  if (rc < 0) {
1060                    new_rank_v = rank;
1061                  } else {
1062                    LargeCost nrc = rc / _epsilon;
1063                    new_rank_v = 0;
1064                    if (nrc < LargeCost(_max_rank)) {
1065                      new_rank_v = rank - 1 - static_cast<int>(nrc);
1066                    }
1067                  }
1068
1069                  // Change the rank of node v
1070                  if (new_rank_v > old_rank_v) {
1071                    _rank[v] = new_rank_v;
1072
1073                    // Remove v from its old bucket
1074                    if (old_rank_v > 0) {
1075                      if (_buckets[old_rank_v] == v) {
1076                        _buckets[old_rank_v] = _bucket_next[v];
1077                      } else {
1078                        int pv = _bucket_prev[v], nv = _bucket_next[v];
1079                        _bucket_next[pv] = nv;
1080                        _bucket_prev[nv] = pv;
1081                      }
1082                    }
1083
1084                    // Insert v into its new bucket
1085                    int nv = _buckets[new_rank_v];
1086                    _bucket_next[v] = nv;
1087                    _bucket_prev[nv] = v;
1088                    _buckets[new_rank_v] = v;
1089                  }
1090                }
1091              }
1092            }
1093
1094            // Refine potential of node u
1095            _pi[u] -= rank * _epsilon;
1096          }
1097        }
1098
1099      }
1100
1101      return false;
1102    }
1103
1104    // Find and cancel cycles in the admissible network and
1105    // determine topological order using DFS
1106    bool topologicalSort(IntVector &stack, int &stack_top) {
1107      const int MAX_CYCLE_CANCEL = 1;
1108
1109      BoolVector reached(_res_node_num, false);
1110      BoolVector processed(_res_node_num, false);
1111      IntVector pred(_res_node_num);
1112      for (int i = 0; i != _res_node_num; ++i) {
1113        _next_out[i] = _first_out[i];
1114      }
1115      stack_top = -1;
1116
1117      int cycle_cnt = 0;
1118      for (int start = 0; start != _res_node_num; ++start) {
1119        if (reached[start]) continue;
1120
1121        // Start DFS search from this start node
1122        pred[start] = -1;
1123        int tip = start, v;
1124        while (true) {
1125          // Check the outgoing arcs of the current tip node
1126          reached[tip] = true;
1127          LargeCost pi_tip = _pi[tip];
1128          int a, last_out = _first_out[tip+1];
1129          for (a = _next_out[tip]; a != last_out; ++a) {
1130            if (_res_cap[a] > 0) {
1131              v = _target[a];
1132              if (_cost[a] + pi_tip - _pi[v] < 0) {
1133                if (!reached[v]) {
1134                  // A new node is reached
1135                  reached[v] = true;
1136                  pred[v] = tip;
1137                  _next_out[tip] = a;
1138                  tip = v;
1139                  a = _next_out[tip];
1140                  last_out = _first_out[tip+1];
1141                  break;
1142                }
1143                else if (!processed[v]) {
1144                  // A cycle is found
1145                  ++cycle_cnt;
1146                  _next_out[tip] = a;
1147
1148                  // Find the minimum residual capacity along the cycle
1149                  Value d, delta = _res_cap[a];
1150                  int u, delta_node = tip;
1151                  for (u = tip; u != v; ) {
1152                    u = pred[u];
1153                    d = _res_cap[_next_out[u]];
1154                    if (d <= delta) {
1155                      delta = d;
1156                      delta_node = u;
1157                    }
1158                  }
1159
1160                  // Augment along the cycle
1161                  _res_cap[a] -= delta;
1162                  _res_cap[_reverse[a]] += delta;
1163                  for (u = tip; u != v; ) {
1164                    u = pred[u];
1165                    int ca = _next_out[u];
1166                    _res_cap[ca] -= delta;
1167                    _res_cap[_reverse[ca]] += delta;
1168                  }
1169
1170                  // Check the maximum number of cycle canceling
1171                  if (cycle_cnt >= MAX_CYCLE_CANCEL) {
1172                    return false;
1173                  }
1174
1175                  // Roll back search to delta_node
1176                  if (delta_node != tip) {
1177                    for (u = tip; u != delta_node; u = pred[u]) {
1178                      reached[u] = false;
1179                    }
1180                    tip = delta_node;
1181                    a = _next_out[tip] + 1;
1182                    last_out = _first_out[tip+1];
1183                    break;
1184                  }
1185                }
1186              }
1187            }
1188          }
1189
1190          // Step back to the previous node
1191          if (a == last_out) {
1192            processed[tip] = true;
1193            stack[++stack_top] = tip;
1194            tip = pred[tip];
1195            if (tip < 0) {
1196              // Finish DFS from the current start node
1197              break;
1198            }
1199            ++_next_out[tip];
1200          }
1201        }
1202
1203      }
1204
1205      return (cycle_cnt == 0);
1206    }
1207
1208    // Global potential update heuristic
1209    void globalUpdate() {
1210      const int bucket_end = _root + 1;
1211
1212      // Initialize buckets
1213      for (int r = 0; r != _max_rank; ++r) {
1214        _buckets[r] = bucket_end;
1215      }
1216      Value total_excess = 0;
1217      int b0 = bucket_end;
1218      for (int i = 0; i != _res_node_num; ++i) {
1219        if (_excess[i] < 0) {
1220          _rank[i] = 0;
1221          _bucket_next[i] = b0;
1222          _bucket_prev[b0] = i;
1223          b0 = i;
1224        } else {
1225          total_excess += _excess[i];
1226          _rank[i] = _max_rank;
1227        }
1228      }
1229      if (total_excess == 0) return;
1230      _buckets[0] = b0;
1231
1232      // Search the buckets
1233      int r = 0;
1234      for ( ; r != _max_rank; ++r) {
1235        while (_buckets[r] != bucket_end) {
1236          // Remove the first node from the current bucket
1237          int u = _buckets[r];
1238          _buckets[r] = _bucket_next[u];
1239
1240          // Search the incomming arcs of u
1241          LargeCost pi_u = _pi[u];
1242          int last_out = _first_out[u+1];
1243          for (int a = _first_out[u]; a != last_out; ++a) {
1244            int ra = _reverse[a];
1245            if (_res_cap[ra] > 0) {
1246              int v = _source[ra];
1247              int old_rank_v = _rank[v];
1248              if (r < old_rank_v) {
1249                // Compute the new rank of v
1250                LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
1251                int new_rank_v = old_rank_v;
1252                if (nrc < LargeCost(_max_rank)) {
1253                  new_rank_v = r + 1 + static_cast<int>(nrc);
1254                }
1255
1256                // Change the rank of v
1257                if (new_rank_v < old_rank_v) {
1258                  _rank[v] = new_rank_v;
1259                  _next_out[v] = _first_out[v];
1260
1261                  // Remove v from its old bucket
1262                  if (old_rank_v < _max_rank) {
1263                    if (_buckets[old_rank_v] == v) {
1264                      _buckets[old_rank_v] = _bucket_next[v];
1265                    } else {
1266                      int pv = _bucket_prev[v], nv = _bucket_next[v];
1267                      _bucket_next[pv] = nv;
1268                      _bucket_prev[nv] = pv;
1269                    }
1270                  }
1271
1272                  // Insert v into its new bucket
1273                  int nv = _buckets[new_rank_v];
1274                  _bucket_next[v] = nv;
1275                  _bucket_prev[nv] = v;
1276                  _buckets[new_rank_v] = v;
1277                }
1278              }
1279            }
1280          }
1281
1282          // Finish search if there are no more active nodes
1283          if (_excess[u] > 0) {
1284            total_excess -= _excess[u];
1285            if (total_excess <= 0) break;
1286          }
1287        }
1288        if (total_excess <= 0) break;
1289      }
1290
1291      // Relabel nodes
1292      for (int u = 0; u != _res_node_num; ++u) {
1293        int k = std::min(_rank[u], r);
1294        if (k > 0) {
1295          _pi[u] -= _epsilon * k;
1296          _next_out[u] = _first_out[u];
1297        }
1298      }
1299    }
1300
1301    /// Execute the algorithm performing augment and relabel operations
1302    void startAugment(int max_length) {
1303      // Paramters for heuristics
1304      const int PRICE_REFINEMENT_LIMIT = 2;
1305      const double GLOBAL_UPDATE_FACTOR = 1.0;
1306      const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
1307        (_res_node_num + _sup_node_num * _sup_node_num));
1308      int next_global_update_limit = global_update_skip;
1309
1310      // Perform cost scaling phases
1311      IntVector path;
1312      BoolVector path_arc(_res_arc_num, false);
1313      int relabel_cnt = 0;
1314      int eps_phase_cnt = 0;
1315      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1316                                        1 : _epsilon / _alpha )
1317      {
1318        ++eps_phase_cnt;
1319
1320        // Price refinement heuristic
1321        if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) {
1322          if (priceRefinement()) continue;
1323        }
1324
1325        // Initialize current phase
1326        initPhase();
1327
1328        // Perform partial augment and relabel operations
1329        while (true) {
1330          // Select an active node (FIFO selection)
1331          while (_active_nodes.size() > 0 &&
1332                 _excess[_active_nodes.front()] <= 0) {
1333            _active_nodes.pop_front();
1334          }
1335          if (_active_nodes.size() == 0) break;
1336          int start = _active_nodes.front();
1337
1338          // Find an augmenting path from the start node
1339          int tip = start;
1340          while (int(path.size()) < max_length && _excess[tip] >= 0) {
1341            int u;
1342            LargeCost rc, min_red_cost = std::numeric_limits<LargeCost>::max();
1343            LargeCost pi_tip = _pi[tip];
1344            int last_out = _first_out[tip+1];
1345            for (int a = _next_out[tip]; a != last_out; ++a) {
1346              if (_res_cap[a] > 0) {
1347                u = _target[a];
1348                rc = _cost[a] + pi_tip - _pi[u];
1349                if (rc < 0) {
1350                  path.push_back(a);
1351                  _next_out[tip] = a;
1352                  if (path_arc[a]) {
1353                    goto augment;   // a cycle is found, stop path search
1354                  }
1355                  tip = u;
1356                  path_arc[a] = true;
1357                  goto next_step;
1358                }
1359                else if (rc < min_red_cost) {
1360                  min_red_cost = rc;
1361                }
1362              }
1363            }
1364
1365            // Relabel tip node
1366            if (tip != start) {
1367              int ra = _reverse[path.back()];
1368              min_red_cost =
1369                std::min(min_red_cost, _cost[ra] + pi_tip - _pi[_target[ra]]);
1370            }
1371            last_out = _next_out[tip];
1372            for (int a = _first_out[tip]; a != last_out; ++a) {
1373              if (_res_cap[a] > 0) {
1374                rc = _cost[a] + pi_tip - _pi[_target[a]];
1375                if (rc < min_red_cost) {
1376                  min_red_cost = rc;
1377                }
1378              }
1379            }
1380            _pi[tip] -= min_red_cost + _epsilon;
1381            _next_out[tip] = _first_out[tip];
1382            ++relabel_cnt;
1383
1384            // Step back
1385            if (tip != start) {
1386              int pa = path.back();
1387              path_arc[pa] = false;
1388              tip = _source[pa];
1389              path.pop_back();
1390            }
1391
1392          next_step: ;
1393          }
1394
1395          // Augment along the found path (as much flow as possible)
1396        augment:
1397          Value delta;
1398          int pa, u, v = start;
1399          for (int i = 0; i != int(path.size()); ++i) {
1400            pa = path[i];
1401            u = v;
1402            v = _target[pa];
1403            path_arc[pa] = false;
1404            delta = std::min(_res_cap[pa], _excess[u]);
1405            _res_cap[pa] -= delta;
1406            _res_cap[_reverse[pa]] += delta;
1407            _excess[u] -= delta;
1408            _excess[v] += delta;
1409            if (_excess[v] > 0 && _excess[v] <= delta) {
1410              _active_nodes.push_back(v);
1411            }
1412          }
1413          path.clear();
1414
1415          // Global update heuristic
1416          if (relabel_cnt >= next_global_update_limit) {
1417            globalUpdate();
1418            next_global_update_limit += global_update_skip;
1419          }
1420        }
1421
1422      }
1423
1424    }
1425
1426    /// Execute the algorithm performing push and relabel operations
1427    void startPush() {
1428      // Paramters for heuristics
1429      const int PRICE_REFINEMENT_LIMIT = 2;
1430      const double GLOBAL_UPDATE_FACTOR = 2.0;
1431
1432      const int global_update_skip = static_cast<int>(GLOBAL_UPDATE_FACTOR *
1433        (_res_node_num + _sup_node_num * _sup_node_num));
1434      int next_global_update_limit = global_update_skip;
1435
1436      // Perform cost scaling phases
1437      BoolVector hyper(_res_node_num, false);
1438      LargeCostVector hyper_cost(_res_node_num);
1439      int relabel_cnt = 0;
1440      int eps_phase_cnt = 0;
1441      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1442                                        1 : _epsilon / _alpha )
1443      {
1444        ++eps_phase_cnt;
1445
1446        // Price refinement heuristic
1447        if (eps_phase_cnt >= PRICE_REFINEMENT_LIMIT) {
1448          if (priceRefinement()) continue;
1449        }
1450
1451        // Initialize current phase
1452        initPhase();
1453
1454        // Perform push and relabel operations
1455        while (_active_nodes.size() > 0) {
1456          LargeCost min_red_cost, rc, pi_n;
1457          Value delta;
1458          int n, t, a, last_out = _res_arc_num;
1459
1460        next_node:
1461          // Select an active node (FIFO selection)
1462          n = _active_nodes.front();
1463          last_out = _first_out[n+1];
1464          pi_n = _pi[n];
1465
1466          // Perform push operations if there are admissible arcs
1467          if (_excess[n] > 0) {
1468            for (a = _next_out[n]; a != last_out; ++a) {
1469              if (_res_cap[a] > 0 &&
1470                  _cost[a] + pi_n - _pi[_target[a]] < 0) {
1471                delta = std::min(_res_cap[a], _excess[n]);
1472                t = _target[a];
1473
1474                // Push-look-ahead heuristic
1475                Value ahead = -_excess[t];
1476                int last_out_t = _first_out[t+1];
1477                LargeCost pi_t = _pi[t];
1478                for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
1479                  if (_res_cap[ta] > 0 &&
1480                      _cost[ta] + pi_t - _pi[_target[ta]] < 0)
1481                    ahead += _res_cap[ta];
1482                  if (ahead >= delta) break;
1483                }
1484                if (ahead < 0) ahead = 0;
1485
1486                // Push flow along the arc
1487                if (ahead < delta && !hyper[t]) {
1488                  _res_cap[a] -= ahead;
1489                  _res_cap[_reverse[a]] += ahead;
1490                  _excess[n] -= ahead;
1491                  _excess[t] += ahead;
1492                  _active_nodes.push_front(t);
1493                  hyper[t] = true;
1494                  hyper_cost[t] = _cost[a] + pi_n - pi_t;
1495                  _next_out[n] = a;
1496                  goto next_node;
1497                } else {
1498                  _res_cap[a] -= delta;
1499                  _res_cap[_reverse[a]] += delta;
1500                  _excess[n] -= delta;
1501                  _excess[t] += delta;
1502                  if (_excess[t] > 0 && _excess[t] <= delta)
1503                    _active_nodes.push_back(t);
1504                }
1505
1506                if (_excess[n] == 0) {
1507                  _next_out[n] = a;
1508                  goto remove_nodes;
1509                }
1510              }
1511            }
1512            _next_out[n] = a;
1513          }
1514
1515          // Relabel the node if it is still active (or hyper)
1516          if (_excess[n] > 0 || hyper[n]) {
1517             min_red_cost = hyper[n] ? -hyper_cost[n] :
1518               std::numeric_limits<LargeCost>::max();
1519            for (int a = _first_out[n]; a != last_out; ++a) {
1520              if (_res_cap[a] > 0) {
1521                rc = _cost[a] + pi_n - _pi[_target[a]];
1522                if (rc < min_red_cost) {
1523                  min_red_cost = rc;
1524                }
1525              }
1526            }
1527            _pi[n] -= min_red_cost + _epsilon;
1528            _next_out[n] = _first_out[n];
1529            hyper[n] = false;
1530            ++relabel_cnt;
1531          }
1532
1533          // Remove nodes that are not active nor hyper
1534        remove_nodes:
1535          while ( _active_nodes.size() > 0 &&
1536                  _excess[_active_nodes.front()] <= 0 &&
1537                  !hyper[_active_nodes.front()] ) {
1538            _active_nodes.pop_front();
1539          }
1540
1541          // Global update heuristic
1542          if (relabel_cnt >= next_global_update_limit) {
1543            globalUpdate();
1544            for (int u = 0; u != _res_node_num; ++u)
1545              hyper[u] = false;
1546            next_global_update_limit += global_update_skip;
1547          }
1548        }
1549      }
1550    }
1551
1552  }; //class CostScaling
1553
1554  ///@}
1555
1556} //namespace lemon
1557
1558#endif //LEMON_COST_SCALING_H
Note: See TracBrowser for help on using the repository browser.