COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/cost_scaling.h @ 821:072ec8120958

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

Small bug fixes (#180)

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