COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/cost_scaling.h @ 809:22bb98ca0101

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

Entirely rework CostScaling? (#180)

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