COIN-OR::LEMON - Graph Library

source: lemon/lemon/cost_scaling.h @ 898:75c97c3786d6

Last change on this file since 898:75c97c3786d6 was 898:75c97c3786d6, checked in by Peter Kovacs <kpeter@…>, 14 years ago

Handle graph changes in the MCF algorithms (#327)

The reset() functions are renamed to resetParams() and the new reset()
functions handle the graph chnages, as well.

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