COIN-OR::LEMON - Graph Library

source: lemon/lemon/cycle_canceling.h @ 1216:dbaf21739390

Last change on this file since 1216:dbaf21739390 was 1179:f6f6896a4724, checked in by Peter Kovacs <kpeter@…>, 11 years ago

Ensure strongly polynomial running time for CycleCanceling? (#436)
The number of iterations performed by Howard's algorithm is limited.
If the limit is reached, a strongly polynomial implementation,
HartmannOrlinMmc? is executed to find a minimum mean cycle.
This iteration limit is typically not reached, thus the combined
method is practically equivalent to Howard's algorithm, while it
also ensures the strongly polynomial time bound.

File size: 40.2 KB
Line 
1/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library.
4 *
5 * Copyright (C) 2003-2010
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
12 *
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_CYCLE_CANCELING_H
20#define LEMON_CYCLE_CANCELING_H
21
22/// \ingroup min_cost_flow_algs
23/// \file
24/// \brief Cycle-canceling algorithms for finding a minimum cost flow.
25
26#include <vector>
27#include <limits>
28
29#include <lemon/core.h>
30#include <lemon/maps.h>
31#include <lemon/path.h>
32#include <lemon/math.h>
33#include <lemon/static_graph.h>
34#include <lemon/adaptors.h>
35#include <lemon/circulation.h>
36#include <lemon/bellman_ford.h>
37#include <lemon/howard_mmc.h>
38#include <lemon/hartmann_orlin_mmc.h>
39
40namespace lemon {
41
42  /// \addtogroup min_cost_flow_algs
43  /// @{
44
45  /// \brief Implementation of cycle-canceling algorithms for
46  /// finding a \ref min_cost_flow "minimum cost flow".
47  ///
48  /// \ref CycleCanceling implements three different cycle-canceling
49  /// algorithms for finding a \ref min_cost_flow "minimum cost flow"
50  /// \ref amo93networkflows, \ref klein67primal,
51  /// \ref goldberg89cyclecanceling.
52  /// The most efficent one is the \ref CANCEL_AND_TIGHTEN
53  /// "Cancel-and-Tighten" algorithm, thus it is the default method.
54  /// It runs in strongly polynomial time, but in practice, it is typically
55  /// orders of magnitude slower than the scaling algorithms and
56  /// \ref NetworkSimplex.
57  /// (For more information, see \ref min_cost_flow_algs "the module page".)
58  ///
59  /// Most of the parameters of the problem (except for the digraph)
60  /// can be given using separate functions, and the algorithm can be
61  /// executed using the \ref run() function. If some parameters are not
62  /// specified, then default values will be used.
63  ///
64  /// \tparam GR The digraph type the algorithm runs on.
65  /// \tparam V The number type used for flow amounts, capacity bounds
66  /// and supply values in the algorithm. By default, it is \c int.
67  /// \tparam C The number type used for costs and potentials in the
68  /// algorithm. By default, it is the same as \c V.
69  ///
70  /// \warning Both \c V and \c C must be signed number types.
71  /// \warning All input data (capacities, supply values, and costs) must
72  /// be integer.
73  /// \warning This algorithm does not support negative costs for
74  /// arcs having infinite upper bound.
75  ///
76  /// \note For more information about the three available methods,
77  /// see \ref Method.
78#ifdef DOXYGEN
79  template <typename GR, typename V, typename C>
80#else
81  template <typename GR, typename V = int, typename C = V>
82#endif
83  class CycleCanceling
84  {
85  public:
86
87    /// The type of the digraph
88    typedef GR Digraph;
89    /// The type of the flow amounts, capacity bounds and supply values
90    typedef V Value;
91    /// The type of the arc costs
92    typedef C Cost;
93
94  public:
95
96    /// \brief Problem type constants for the \c run() function.
97    ///
98    /// Enum type containing the problem type constants that can be
99    /// returned by the \ref run() function of the algorithm.
100    enum ProblemType {
101      /// The problem has no feasible solution (flow).
102      INFEASIBLE,
103      /// The problem has optimal solution (i.e. it is feasible and
104      /// bounded), and the algorithm has found optimal flow and node
105      /// potentials (primal and dual solutions).
106      OPTIMAL,
107      /// The digraph contains an arc of negative cost and infinite
108      /// upper bound. It means that the objective function is unbounded
109      /// on that arc, however, note that it could actually be bounded
110      /// over the feasible flows, but this algroithm cannot handle
111      /// these cases.
112      UNBOUNDED
113    };
114
115    /// \brief Constants for selecting the used method.
116    ///
117    /// Enum type containing constants for selecting the used method
118    /// for the \ref run() function.
119    ///
120    /// \ref CycleCanceling provides three different cycle-canceling
121    /// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel-and-Tighten"
122    /// is used, which is by far the most efficient and the most robust.
123    /// However, the other methods can be selected using the \ref run()
124    /// function with the proper parameter.
125    enum Method {
126      /// A simple cycle-canceling method, which uses the
127      /// \ref BellmanFord "Bellman-Ford" algorithm for detecting negative
128      /// cycles in the residual network.
129      /// The number of Bellman-Ford iterations is bounded by a successively
130      /// increased limit.
131      SIMPLE_CYCLE_CANCELING,
132      /// The "Minimum Mean Cycle-Canceling" algorithm, which is a
133      /// well-known strongly polynomial method
134      /// \ref goldberg89cyclecanceling. It improves along a
135      /// \ref min_mean_cycle "minimum mean cycle" in each iteration.
136      /// Its running time complexity is O(n<sup>2</sup>e<sup>3</sup>log(n)).
137      MINIMUM_MEAN_CYCLE_CANCELING,
138      /// The "Cancel-and-Tighten" algorithm, which can be viewed as an
139      /// improved version of the previous method
140      /// \ref goldberg89cyclecanceling.
141      /// It is faster both in theory and in practice, its running time
142      /// complexity is O(n<sup>2</sup>e<sup>2</sup>log(n)).
143      CANCEL_AND_TIGHTEN
144    };
145
146  private:
147
148    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
149
150    typedef std::vector<int> IntVector;
151    typedef std::vector<double> DoubleVector;
152    typedef std::vector<Value> ValueVector;
153    typedef std::vector<Cost> CostVector;
154    typedef std::vector<char> BoolVector;
155    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
156
157  private:
158
159    template <typename KT, typename VT>
160    class StaticVectorMap {
161    public:
162      typedef KT Key;
163      typedef VT Value;
164
165      StaticVectorMap(std::vector<Value>& v) : _v(v) {}
166
167      const Value& operator[](const Key& key) const {
168        return _v[StaticDigraph::id(key)];
169      }
170
171      Value& operator[](const Key& key) {
172        return _v[StaticDigraph::id(key)];
173      }
174
175      void set(const Key& key, const Value& val) {
176        _v[StaticDigraph::id(key)] = val;
177      }
178
179    private:
180      std::vector<Value>& _v;
181    };
182
183    typedef StaticVectorMap<StaticDigraph::Node, Cost> CostNodeMap;
184    typedef StaticVectorMap<StaticDigraph::Arc, Cost> CostArcMap;
185
186  private:
187
188
189    // Data related to the underlying digraph
190    const GR &_graph;
191    int _node_num;
192    int _arc_num;
193    int _res_node_num;
194    int _res_arc_num;
195    int _root;
196
197    // Parameters of the problem
198    bool _have_lower;
199    Value _sum_supply;
200
201    // Data structures for storing the digraph
202    IntNodeMap _node_id;
203    IntArcMap _arc_idf;
204    IntArcMap _arc_idb;
205    IntVector _first_out;
206    BoolVector _forward;
207    IntVector _source;
208    IntVector _target;
209    IntVector _reverse;
210
211    // Node and arc data
212    ValueVector _lower;
213    ValueVector _upper;
214    CostVector _cost;
215    ValueVector _supply;
216
217    ValueVector _res_cap;
218    CostVector _pi;
219
220    // Data for a StaticDigraph structure
221    typedef std::pair<int, int> IntPair;
222    StaticDigraph _sgr;
223    std::vector<IntPair> _arc_vec;
224    std::vector<Cost> _cost_vec;
225    IntVector _id_vec;
226    CostArcMap _cost_map;
227    CostNodeMap _pi_map;
228
229  public:
230
231    /// \brief Constant for infinite upper bounds (capacities).
232    ///
233    /// Constant for infinite upper bounds (capacities).
234    /// It is \c std::numeric_limits<Value>::infinity() if available,
235    /// \c std::numeric_limits<Value>::max() otherwise.
236    const Value INF;
237
238  public:
239
240    /// \brief Constructor.
241    ///
242    /// The constructor of the class.
243    ///
244    /// \param graph The digraph the algorithm runs on.
245    CycleCanceling(const GR& graph) :
246      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
247      _cost_map(_cost_vec), _pi_map(_pi),
248      INF(std::numeric_limits<Value>::has_infinity ?
249          std::numeric_limits<Value>::infinity() :
250          std::numeric_limits<Value>::max())
251    {
252      // Check the number types
253      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
254        "The flow type of CycleCanceling must be signed");
255      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
256        "The cost type of CycleCanceling must be signed");
257
258      // Reset data structures
259      reset();
260    }
261
262    /// \name Parameters
263    /// The parameters of the algorithm can be specified using these
264    /// functions.
265
266    /// @{
267
268    /// \brief Set the lower bounds on the arcs.
269    ///
270    /// This function sets the lower bounds on the arcs.
271    /// If it is not used before calling \ref run(), the lower bounds
272    /// will be set to zero on all arcs.
273    ///
274    /// \param map An arc map storing the lower bounds.
275    /// Its \c Value type must be convertible to the \c Value type
276    /// of the algorithm.
277    ///
278    /// \return <tt>(*this)</tt>
279    template <typename LowerMap>
280    CycleCanceling& lowerMap(const LowerMap& map) {
281      _have_lower = true;
282      for (ArcIt a(_graph); a != INVALID; ++a) {
283        _lower[_arc_idf[a]] = map[a];
284        _lower[_arc_idb[a]] = map[a];
285      }
286      return *this;
287    }
288
289    /// \brief Set the upper bounds (capacities) on the arcs.
290    ///
291    /// This function sets the upper bounds (capacities) on the arcs.
292    /// If it is not used before calling \ref run(), the upper bounds
293    /// will be set to \ref INF on all arcs (i.e. the flow value will be
294    /// unbounded from above).
295    ///
296    /// \param map An arc map storing the upper bounds.
297    /// Its \c Value type must be convertible to the \c Value type
298    /// of the algorithm.
299    ///
300    /// \return <tt>(*this)</tt>
301    template<typename UpperMap>
302    CycleCanceling& upperMap(const UpperMap& map) {
303      for (ArcIt a(_graph); a != INVALID; ++a) {
304        _upper[_arc_idf[a]] = map[a];
305      }
306      return *this;
307    }
308
309    /// \brief Set the costs of the arcs.
310    ///
311    /// This function sets the costs of the arcs.
312    /// If it is not used before calling \ref run(), the costs
313    /// will be set to \c 1 on all arcs.
314    ///
315    /// \param map An arc map storing the costs.
316    /// Its \c Value type must be convertible to the \c Cost type
317    /// of the algorithm.
318    ///
319    /// \return <tt>(*this)</tt>
320    template<typename CostMap>
321    CycleCanceling& costMap(const CostMap& map) {
322      for (ArcIt a(_graph); a != INVALID; ++a) {
323        _cost[_arc_idf[a]] =  map[a];
324        _cost[_arc_idb[a]] = -map[a];
325      }
326      return *this;
327    }
328
329    /// \brief Set the supply values of the nodes.
330    ///
331    /// This function sets the supply values of the nodes.
332    /// If neither this function nor \ref stSupply() is used before
333    /// calling \ref run(), the supply of each node will be set to zero.
334    ///
335    /// \param map A node map storing the supply values.
336    /// Its \c Value type must be convertible to the \c Value type
337    /// of the algorithm.
338    ///
339    /// \return <tt>(*this)</tt>
340    template<typename SupplyMap>
341    CycleCanceling& supplyMap(const SupplyMap& map) {
342      for (NodeIt n(_graph); n != INVALID; ++n) {
343        _supply[_node_id[n]] = map[n];
344      }
345      return *this;
346    }
347
348    /// \brief Set single source and target nodes and a supply value.
349    ///
350    /// This function sets a single source node and a single target node
351    /// and the required flow value.
352    /// If neither this function nor \ref supplyMap() is used before
353    /// calling \ref run(), the supply of each node will be set to zero.
354    ///
355    /// Using this function has the same effect as using \ref supplyMap()
356    /// with a map in which \c k is assigned to \c s, \c -k is
357    /// assigned to \c t and all other nodes have zero supply value.
358    ///
359    /// \param s The source node.
360    /// \param t The target node.
361    /// \param k The required amount of flow from node \c s to node \c t
362    /// (i.e. the supply of \c s and the demand of \c t).
363    ///
364    /// \return <tt>(*this)</tt>
365    CycleCanceling& stSupply(const Node& s, const Node& t, Value k) {
366      for (int i = 0; i != _res_node_num; ++i) {
367        _supply[i] = 0;
368      }
369      _supply[_node_id[s]] =  k;
370      _supply[_node_id[t]] = -k;
371      return *this;
372    }
373
374    /// @}
375
376    /// \name Execution control
377    /// The algorithm can be executed using \ref run().
378
379    /// @{
380
381    /// \brief Run the algorithm.
382    ///
383    /// This function runs the algorithm.
384    /// The paramters can be specified using functions \ref lowerMap(),
385    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
386    /// For example,
387    /// \code
388    ///   CycleCanceling<ListDigraph> cc(graph);
389    ///   cc.lowerMap(lower).upperMap(upper).costMap(cost)
390    ///     .supplyMap(sup).run();
391    /// \endcode
392    ///
393    /// This function can be called more than once. All the given parameters
394    /// are kept for the next call, unless \ref resetParams() or \ref reset()
395    /// is used, thus only the modified parameters have to be set again.
396    /// If the underlying digraph was also modified after the construction
397    /// of the class (or the last \ref reset() call), then the \ref reset()
398    /// function must be called.
399    ///
400    /// \param method The cycle-canceling method that will be used.
401    /// For more information, see \ref Method.
402    ///
403    /// \return \c INFEASIBLE if no feasible flow exists,
404    /// \n \c OPTIMAL if the problem has optimal solution
405    /// (i.e. it is feasible and bounded), and the algorithm has found
406    /// optimal flow and node potentials (primal and dual solutions),
407    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
408    /// and infinite upper bound. It means that the objective function
409    /// is unbounded on that arc, however, note that it could actually be
410    /// bounded over the feasible flows, but this algroithm cannot handle
411    /// these cases.
412    ///
413    /// \see ProblemType, Method
414    /// \see resetParams(), reset()
415    ProblemType run(Method method = CANCEL_AND_TIGHTEN) {
416      ProblemType pt = init();
417      if (pt != OPTIMAL) return pt;
418      start(method);
419      return OPTIMAL;
420    }
421
422    /// \brief Reset all the parameters that have been given before.
423    ///
424    /// This function resets all the paramaters that have been given
425    /// before using functions \ref lowerMap(), \ref upperMap(),
426    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
427    ///
428    /// It is useful for multiple \ref run() calls. Basically, all the given
429    /// parameters are kept for the next \ref run() call, unless
430    /// \ref resetParams() or \ref reset() is used.
431    /// If the underlying digraph was also modified after the construction
432    /// of the class or the last \ref reset() call, then the \ref reset()
433    /// function must be used, otherwise \ref resetParams() is sufficient.
434    ///
435    /// For example,
436    /// \code
437    ///   CycleCanceling<ListDigraph> cs(graph);
438    ///
439    ///   // First run
440    ///   cc.lowerMap(lower).upperMap(upper).costMap(cost)
441    ///     .supplyMap(sup).run();
442    ///
443    ///   // Run again with modified cost map (resetParams() is not called,
444    ///   // so only the cost map have to be set again)
445    ///   cost[e] += 100;
446    ///   cc.costMap(cost).run();
447    ///
448    ///   // Run again from scratch using resetParams()
449    ///   // (the lower bounds will be set to zero on all arcs)
450    ///   cc.resetParams();
451    ///   cc.upperMap(capacity).costMap(cost)
452    ///     .supplyMap(sup).run();
453    /// \endcode
454    ///
455    /// \return <tt>(*this)</tt>
456    ///
457    /// \see reset(), run()
458    CycleCanceling& resetParams() {
459      for (int i = 0; i != _res_node_num; ++i) {
460        _supply[i] = 0;
461      }
462      int limit = _first_out[_root];
463      for (int j = 0; j != limit; ++j) {
464        _lower[j] = 0;
465        _upper[j] = INF;
466        _cost[j] = _forward[j] ? 1 : -1;
467      }
468      for (int j = limit; j != _res_arc_num; ++j) {
469        _lower[j] = 0;
470        _upper[j] = INF;
471        _cost[j] = 0;
472        _cost[_reverse[j]] = 0;
473      }
474      _have_lower = false;
475      return *this;
476    }
477
478    /// \brief Reset the internal data structures and all the parameters
479    /// that have been given before.
480    ///
481    /// This function resets the internal data structures and all the
482    /// paramaters that have been given before using functions \ref lowerMap(),
483    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
484    ///
485    /// It is useful for multiple \ref run() calls. Basically, all the given
486    /// parameters are kept for the next \ref run() call, unless
487    /// \ref resetParams() or \ref reset() is used.
488    /// If the underlying digraph was also modified after the construction
489    /// of the class or the last \ref reset() call, then the \ref reset()
490    /// function must be used, otherwise \ref resetParams() is sufficient.
491    ///
492    /// See \ref resetParams() for examples.
493    ///
494    /// \return <tt>(*this)</tt>
495    ///
496    /// \see resetParams(), run()
497    CycleCanceling& reset() {
498      // Resize vectors
499      _node_num = countNodes(_graph);
500      _arc_num = countArcs(_graph);
501      _res_node_num = _node_num + 1;
502      _res_arc_num = 2 * (_arc_num + _node_num);
503      _root = _node_num;
504
505      _first_out.resize(_res_node_num + 1);
506      _forward.resize(_res_arc_num);
507      _source.resize(_res_arc_num);
508      _target.resize(_res_arc_num);
509      _reverse.resize(_res_arc_num);
510
511      _lower.resize(_res_arc_num);
512      _upper.resize(_res_arc_num);
513      _cost.resize(_res_arc_num);
514      _supply.resize(_res_node_num);
515
516      _res_cap.resize(_res_arc_num);
517      _pi.resize(_res_node_num);
518
519      _arc_vec.reserve(_res_arc_num);
520      _cost_vec.reserve(_res_arc_num);
521      _id_vec.reserve(_res_arc_num);
522
523      // Copy the graph
524      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
525      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
526        _node_id[n] = i;
527      }
528      i = 0;
529      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
530        _first_out[i] = j;
531        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
532          _arc_idf[a] = j;
533          _forward[j] = true;
534          _source[j] = i;
535          _target[j] = _node_id[_graph.runningNode(a)];
536        }
537        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
538          _arc_idb[a] = j;
539          _forward[j] = false;
540          _source[j] = i;
541          _target[j] = _node_id[_graph.runningNode(a)];
542        }
543        _forward[j] = false;
544        _source[j] = i;
545        _target[j] = _root;
546        _reverse[j] = k;
547        _forward[k] = true;
548        _source[k] = _root;
549        _target[k] = i;
550        _reverse[k] = j;
551        ++j; ++k;
552      }
553      _first_out[i] = j;
554      _first_out[_res_node_num] = k;
555      for (ArcIt a(_graph); a != INVALID; ++a) {
556        int fi = _arc_idf[a];
557        int bi = _arc_idb[a];
558        _reverse[fi] = bi;
559        _reverse[bi] = fi;
560      }
561
562      // Reset parameters
563      resetParams();
564      return *this;
565    }
566
567    /// @}
568
569    /// \name Query Functions
570    /// The results of the algorithm can be obtained using these
571    /// functions.\n
572    /// The \ref run() function must be called before using them.
573
574    /// @{
575
576    /// \brief Return the total cost of the found flow.
577    ///
578    /// This function returns the total cost of the found flow.
579    /// Its complexity is O(e).
580    ///
581    /// \note The return type of the function can be specified as a
582    /// template parameter. For example,
583    /// \code
584    ///   cc.totalCost<double>();
585    /// \endcode
586    /// It is useful if the total cost cannot be stored in the \c Cost
587    /// type of the algorithm, which is the default return type of the
588    /// function.
589    ///
590    /// \pre \ref run() must be called before using this function.
591    template <typename Number>
592    Number totalCost() const {
593      Number c = 0;
594      for (ArcIt a(_graph); a != INVALID; ++a) {
595        int i = _arc_idb[a];
596        c += static_cast<Number>(_res_cap[i]) *
597             (-static_cast<Number>(_cost[i]));
598      }
599      return c;
600    }
601
602#ifndef DOXYGEN
603    Cost totalCost() const {
604      return totalCost<Cost>();
605    }
606#endif
607
608    /// \brief Return the flow on the given arc.
609    ///
610    /// This function returns the flow on the given arc.
611    ///
612    /// \pre \ref run() must be called before using this function.
613    Value flow(const Arc& a) const {
614      return _res_cap[_arc_idb[a]];
615    }
616
617    /// \brief Copy the flow values (the primal solution) into the
618    /// given map.
619    ///
620    /// This function copies the flow value on each arc into the given
621    /// map. The \c Value type of the algorithm must be convertible to
622    /// the \c Value type of the map.
623    ///
624    /// \pre \ref run() must be called before using this function.
625    template <typename FlowMap>
626    void flowMap(FlowMap &map) const {
627      for (ArcIt a(_graph); a != INVALID; ++a) {
628        map.set(a, _res_cap[_arc_idb[a]]);
629      }
630    }
631
632    /// \brief Return the potential (dual value) of the given node.
633    ///
634    /// This function returns the potential (dual value) of the
635    /// given node.
636    ///
637    /// \pre \ref run() must be called before using this function.
638    Cost potential(const Node& n) const {
639      return static_cast<Cost>(_pi[_node_id[n]]);
640    }
641
642    /// \brief Copy the potential values (the dual solution) into the
643    /// given map.
644    ///
645    /// This function copies the potential (dual value) of each node
646    /// into the given map.
647    /// The \c Cost type of the algorithm must be convertible to the
648    /// \c Value type of the map.
649    ///
650    /// \pre \ref run() must be called before using this function.
651    template <typename PotentialMap>
652    void potentialMap(PotentialMap &map) const {
653      for (NodeIt n(_graph); n != INVALID; ++n) {
654        map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
655      }
656    }
657
658    /// @}
659
660  private:
661
662    // Initialize the algorithm
663    ProblemType init() {
664      if (_res_node_num <= 1) return INFEASIBLE;
665
666      // Check the sum of supply values
667      _sum_supply = 0;
668      for (int i = 0; i != _root; ++i) {
669        _sum_supply += _supply[i];
670      }
671      if (_sum_supply > 0) return INFEASIBLE;
672
673
674      // Initialize vectors
675      for (int i = 0; i != _res_node_num; ++i) {
676        _pi[i] = 0;
677      }
678      ValueVector excess(_supply);
679
680      // Remove infinite upper bounds and check negative arcs
681      const Value MAX = std::numeric_limits<Value>::max();
682      int last_out;
683      if (_have_lower) {
684        for (int i = 0; i != _root; ++i) {
685          last_out = _first_out[i+1];
686          for (int j = _first_out[i]; j != last_out; ++j) {
687            if (_forward[j]) {
688              Value c = _cost[j] < 0 ? _upper[j] : _lower[j];
689              if (c >= MAX) return UNBOUNDED;
690              excess[i] -= c;
691              excess[_target[j]] += c;
692            }
693          }
694        }
695      } else {
696        for (int i = 0; i != _root; ++i) {
697          last_out = _first_out[i+1];
698          for (int j = _first_out[i]; j != last_out; ++j) {
699            if (_forward[j] && _cost[j] < 0) {
700              Value c = _upper[j];
701              if (c >= MAX) return UNBOUNDED;
702              excess[i] -= c;
703              excess[_target[j]] += c;
704            }
705          }
706        }
707      }
708      Value ex, max_cap = 0;
709      for (int i = 0; i != _res_node_num; ++i) {
710        ex = excess[i];
711        if (ex < 0) max_cap -= ex;
712      }
713      for (int j = 0; j != _res_arc_num; ++j) {
714        if (_upper[j] >= MAX) _upper[j] = max_cap;
715      }
716
717      // Initialize maps for Circulation and remove non-zero lower bounds
718      ConstMap<Arc, Value> low(0);
719      typedef typename Digraph::template ArcMap<Value> ValueArcMap;
720      typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
721      ValueArcMap cap(_graph), flow(_graph);
722      ValueNodeMap sup(_graph);
723      for (NodeIt n(_graph); n != INVALID; ++n) {
724        sup[n] = _supply[_node_id[n]];
725      }
726      if (_have_lower) {
727        for (ArcIt a(_graph); a != INVALID; ++a) {
728          int j = _arc_idf[a];
729          Value c = _lower[j];
730          cap[a] = _upper[j] - c;
731          sup[_graph.source(a)] -= c;
732          sup[_graph.target(a)] += c;
733        }
734      } else {
735        for (ArcIt a(_graph); a != INVALID; ++a) {
736          cap[a] = _upper[_arc_idf[a]];
737        }
738      }
739
740      // Find a feasible flow using Circulation
741      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
742        circ(_graph, low, cap, sup);
743      if (!circ.flowMap(flow).run()) return INFEASIBLE;
744
745      // Set residual capacities and handle GEQ supply type
746      if (_sum_supply < 0) {
747        for (ArcIt a(_graph); a != INVALID; ++a) {
748          Value fa = flow[a];
749          _res_cap[_arc_idf[a]] = cap[a] - fa;
750          _res_cap[_arc_idb[a]] = fa;
751          sup[_graph.source(a)] -= fa;
752          sup[_graph.target(a)] += fa;
753        }
754        for (NodeIt n(_graph); n != INVALID; ++n) {
755          excess[_node_id[n]] = sup[n];
756        }
757        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
758          int u = _target[a];
759          int ra = _reverse[a];
760          _res_cap[a] = -_sum_supply + 1;
761          _res_cap[ra] = -excess[u];
762          _cost[a] = 0;
763          _cost[ra] = 0;
764        }
765      } else {
766        for (ArcIt a(_graph); a != INVALID; ++a) {
767          Value fa = flow[a];
768          _res_cap[_arc_idf[a]] = cap[a] - fa;
769          _res_cap[_arc_idb[a]] = fa;
770        }
771        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
772          int ra = _reverse[a];
773          _res_cap[a] = 1;
774          _res_cap[ra] = 0;
775          _cost[a] = 0;
776          _cost[ra] = 0;
777        }
778      }
779
780      return OPTIMAL;
781    }
782
783    // Build a StaticDigraph structure containing the current
784    // residual network
785    void buildResidualNetwork() {
786      _arc_vec.clear();
787      _cost_vec.clear();
788      _id_vec.clear();
789      for (int j = 0; j != _res_arc_num; ++j) {
790        if (_res_cap[j] > 0) {
791          _arc_vec.push_back(IntPair(_source[j], _target[j]));
792          _cost_vec.push_back(_cost[j]);
793          _id_vec.push_back(j);
794        }
795      }
796      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
797    }
798
799    // Execute the algorithm and transform the results
800    void start(Method method) {
801      // Execute the algorithm
802      switch (method) {
803        case SIMPLE_CYCLE_CANCELING:
804          startSimpleCycleCanceling();
805          break;
806        case MINIMUM_MEAN_CYCLE_CANCELING:
807          startMinMeanCycleCanceling();
808          break;
809        case CANCEL_AND_TIGHTEN:
810          startCancelAndTighten();
811          break;
812      }
813
814      // Compute node potentials
815      if (method != SIMPLE_CYCLE_CANCELING) {
816        buildResidualNetwork();
817        typename BellmanFord<StaticDigraph, CostArcMap>
818          ::template SetDistMap<CostNodeMap>::Create bf(_sgr, _cost_map);
819        bf.distMap(_pi_map);
820        bf.init(0);
821        bf.start();
822      }
823
824      // Handle non-zero lower bounds
825      if (_have_lower) {
826        int limit = _first_out[_root];
827        for (int j = 0; j != limit; ++j) {
828          if (!_forward[j]) _res_cap[j] += _lower[j];
829        }
830      }
831    }
832
833    // Execute the "Simple Cycle Canceling" method
834    void startSimpleCycleCanceling() {
835      // Constants for computing the iteration limits
836      const int BF_FIRST_LIMIT  = 2;
837      const double BF_LIMIT_FACTOR = 1.5;
838
839      typedef StaticVectorMap<StaticDigraph::Arc, Value> FilterMap;
840      typedef FilterArcs<StaticDigraph, FilterMap> ResDigraph;
841      typedef StaticVectorMap<StaticDigraph::Node, StaticDigraph::Arc> PredMap;
842      typedef typename BellmanFord<ResDigraph, CostArcMap>
843        ::template SetDistMap<CostNodeMap>
844        ::template SetPredMap<PredMap>::Create BF;
845
846      // Build the residual network
847      _arc_vec.clear();
848      _cost_vec.clear();
849      for (int j = 0; j != _res_arc_num; ++j) {
850        _arc_vec.push_back(IntPair(_source[j], _target[j]));
851        _cost_vec.push_back(_cost[j]);
852      }
853      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
854
855      FilterMap filter_map(_res_cap);
856      ResDigraph rgr(_sgr, filter_map);
857      std::vector<int> cycle;
858      std::vector<StaticDigraph::Arc> pred(_res_arc_num);
859      PredMap pred_map(pred);
860      BF bf(rgr, _cost_map);
861      bf.distMap(_pi_map).predMap(pred_map);
862
863      int length_bound = BF_FIRST_LIMIT;
864      bool optimal = false;
865      while (!optimal) {
866        bf.init(0);
867        int iter_num = 0;
868        bool cycle_found = false;
869        while (!cycle_found) {
870          // Perform some iterations of the Bellman-Ford algorithm
871          int curr_iter_num = iter_num + length_bound <= _node_num ?
872            length_bound : _node_num - iter_num;
873          iter_num += curr_iter_num;
874          int real_iter_num = curr_iter_num;
875          for (int i = 0; i < curr_iter_num; ++i) {
876            if (bf.processNextWeakRound()) {
877              real_iter_num = i;
878              break;
879            }
880          }
881          if (real_iter_num < curr_iter_num) {
882            // Optimal flow is found
883            optimal = true;
884            break;
885          } else {
886            // Search for node disjoint negative cycles
887            std::vector<int> state(_res_node_num, 0);
888            int id = 0;
889            for (int u = 0; u != _res_node_num; ++u) {
890              if (state[u] != 0) continue;
891              ++id;
892              int v = u;
893              for (; v != -1 && state[v] == 0; v = pred[v] == INVALID ?
894                   -1 : rgr.id(rgr.source(pred[v]))) {
895                state[v] = id;
896              }
897              if (v != -1 && state[v] == id) {
898                // A negative cycle is found
899                cycle_found = true;
900                cycle.clear();
901                StaticDigraph::Arc a = pred[v];
902                Value d, delta = _res_cap[rgr.id(a)];
903                cycle.push_back(rgr.id(a));
904                while (rgr.id(rgr.source(a)) != v) {
905                  a = pred_map[rgr.source(a)];
906                  d = _res_cap[rgr.id(a)];
907                  if (d < delta) delta = d;
908                  cycle.push_back(rgr.id(a));
909                }
910
911                // Augment along the cycle
912                for (int i = 0; i < int(cycle.size()); ++i) {
913                  int j = cycle[i];
914                  _res_cap[j] -= delta;
915                  _res_cap[_reverse[j]] += delta;
916                }
917              }
918            }
919          }
920
921          // Increase iteration limit if no cycle is found
922          if (!cycle_found) {
923            length_bound = static_cast<int>(length_bound * BF_LIMIT_FACTOR);
924          }
925        }
926      }
927    }
928
929    // Execute the "Minimum Mean Cycle Canceling" method
930    void startMinMeanCycleCanceling() {
931      typedef Path<StaticDigraph> SPath;
932      typedef typename SPath::ArcIt SPathArcIt;
933      typedef typename HowardMmc<StaticDigraph, CostArcMap>
934        ::template SetPath<SPath>::Create HwMmc;
935      typedef typename HartmannOrlinMmc<StaticDigraph, CostArcMap>
936        ::template SetPath<SPath>::Create HoMmc;
937
938      const double HW_ITER_LIMIT_FACTOR = 1.0;
939      const int HW_ITER_LIMIT_MIN_VALUE = 5;
940
941      const int hw_iter_limit =
942          std::max(static_cast<int>(HW_ITER_LIMIT_FACTOR * _node_num),
943                   HW_ITER_LIMIT_MIN_VALUE);
944
945      SPath cycle;
946      HwMmc hw_mmc(_sgr, _cost_map);
947      hw_mmc.cycle(cycle);
948      buildResidualNetwork();
949      while (true) {
950       
951        typename HwMmc::TerminationCause hw_tc =
952            hw_mmc.findCycleMean(hw_iter_limit);
953        if (hw_tc == HwMmc::ITERATION_LIMIT) {
954          // Howard's algorithm reached the iteration limit, start a
955          // strongly polynomial algorithm instead
956          HoMmc ho_mmc(_sgr, _cost_map);
957          ho_mmc.cycle(cycle);
958          // Find a minimum mean cycle (Hartmann-Orlin algorithm)
959          if (!(ho_mmc.findCycleMean() && ho_mmc.cycleCost() < 0)) break;
960          ho_mmc.findCycle();
961        } else {
962          // Find a minimum mean cycle (Howard algorithm)
963          if (!(hw_tc == HwMmc::OPTIMAL && hw_mmc.cycleCost() < 0)) break;
964          hw_mmc.findCycle();
965        }
966       
967        // Compute delta value
968        Value delta = INF;
969        for (SPathArcIt a(cycle); a != INVALID; ++a) {
970          Value d = _res_cap[_id_vec[_sgr.id(a)]];
971          if (d < delta) delta = d;
972        }
973
974        // Augment along the cycle
975        for (SPathArcIt a(cycle); a != INVALID; ++a) {
976          int j = _id_vec[_sgr.id(a)];
977          _res_cap[j] -= delta;
978          _res_cap[_reverse[j]] += delta;
979        }
980
981        // Rebuild the residual network
982        buildResidualNetwork();
983      }
984    }
985
986    // Execute the "Cancel-and-Tighten" method
987    void startCancelAndTighten() {
988      // Constants for the min mean cycle computations
989      const double LIMIT_FACTOR = 1.0;
990      const int MIN_LIMIT = 5;
991      const double HW_ITER_LIMIT_FACTOR = 1.0;
992      const int HW_ITER_LIMIT_MIN_VALUE = 5;
993
994      const int hw_iter_limit =
995          std::max(static_cast<int>(HW_ITER_LIMIT_FACTOR * _node_num),
996                   HW_ITER_LIMIT_MIN_VALUE);
997
998      // Contruct auxiliary data vectors
999      DoubleVector pi(_res_node_num, 0.0);
1000      IntVector level(_res_node_num);
1001      BoolVector reached(_res_node_num);
1002      BoolVector processed(_res_node_num);
1003      IntVector pred_node(_res_node_num);
1004      IntVector pred_arc(_res_node_num);
1005      std::vector<int> stack(_res_node_num);
1006      std::vector<int> proc_vector(_res_node_num);
1007
1008      // Initialize epsilon
1009      double epsilon = 0;
1010      for (int a = 0; a != _res_arc_num; ++a) {
1011        if (_res_cap[a] > 0 && -_cost[a] > epsilon)
1012          epsilon = -_cost[a];
1013      }
1014
1015      // Start phases
1016      Tolerance<double> tol;
1017      tol.epsilon(1e-6);
1018      int limit = int(LIMIT_FACTOR * std::sqrt(double(_res_node_num)));
1019      if (limit < MIN_LIMIT) limit = MIN_LIMIT;
1020      int iter = limit;
1021      while (epsilon * _res_node_num >= 1) {
1022        // Find and cancel cycles in the admissible network using DFS
1023        for (int u = 0; u != _res_node_num; ++u) {
1024          reached[u] = false;
1025          processed[u] = false;
1026        }
1027        int stack_head = -1;
1028        int proc_head = -1;
1029        for (int start = 0; start != _res_node_num; ++start) {
1030          if (reached[start]) continue;
1031
1032          // New start node
1033          reached[start] = true;
1034          pred_arc[start] = -1;
1035          pred_node[start] = -1;
1036
1037          // Find the first admissible outgoing arc
1038          double p = pi[start];
1039          int a = _first_out[start];
1040          int last_out = _first_out[start+1];
1041          for (; a != last_out && (_res_cap[a] == 0 ||
1042               !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
1043          if (a == last_out) {
1044            processed[start] = true;
1045            proc_vector[++proc_head] = start;
1046            continue;
1047          }
1048          stack[++stack_head] = a;
1049
1050          while (stack_head >= 0) {
1051            int sa = stack[stack_head];
1052            int u = _source[sa];
1053            int v = _target[sa];
1054
1055            if (!reached[v]) {
1056              // A new node is reached
1057              reached[v] = true;
1058              pred_node[v] = u;
1059              pred_arc[v] = sa;
1060              p = pi[v];
1061              a = _first_out[v];
1062              last_out = _first_out[v+1];
1063              for (; a != last_out && (_res_cap[a] == 0 ||
1064                   !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
1065              stack[++stack_head] = a == last_out ? -1 : a;
1066            } else {
1067              if (!processed[v]) {
1068                // A cycle is found
1069                int n, w = u;
1070                Value d, delta = _res_cap[sa];
1071                for (n = u; n != v; n = pred_node[n]) {
1072                  d = _res_cap[pred_arc[n]];
1073                  if (d <= delta) {
1074                    delta = d;
1075                    w = pred_node[n];
1076                  }
1077                }
1078
1079                // Augment along the cycle
1080                _res_cap[sa] -= delta;
1081                _res_cap[_reverse[sa]] += delta;
1082                for (n = u; n != v; n = pred_node[n]) {
1083                  int pa = pred_arc[n];
1084                  _res_cap[pa] -= delta;
1085                  _res_cap[_reverse[pa]] += delta;
1086                }
1087                for (n = u; stack_head > 0 && n != w; n = pred_node[n]) {
1088                  --stack_head;
1089                  reached[n] = false;
1090                }
1091                u = w;
1092              }
1093              v = u;
1094
1095              // Find the next admissible outgoing arc
1096              p = pi[v];
1097              a = stack[stack_head] + 1;
1098              last_out = _first_out[v+1];
1099              for (; a != last_out && (_res_cap[a] == 0 ||
1100                   !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
1101              stack[stack_head] = a == last_out ? -1 : a;
1102            }
1103
1104            while (stack_head >= 0 && stack[stack_head] == -1) {
1105              processed[v] = true;
1106              proc_vector[++proc_head] = v;
1107              if (--stack_head >= 0) {
1108                // Find the next admissible outgoing arc
1109                v = _source[stack[stack_head]];
1110                p = pi[v];
1111                a = stack[stack_head] + 1;
1112                last_out = _first_out[v+1];
1113                for (; a != last_out && (_res_cap[a] == 0 ||
1114                     !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
1115                stack[stack_head] = a == last_out ? -1 : a;
1116              }
1117            }
1118          }
1119        }
1120
1121        // Tighten potentials and epsilon
1122        if (--iter > 0) {
1123          for (int u = 0; u != _res_node_num; ++u) {
1124            level[u] = 0;
1125          }
1126          for (int i = proc_head; i > 0; --i) {
1127            int u = proc_vector[i];
1128            double p = pi[u];
1129            int l = level[u] + 1;
1130            int last_out = _first_out[u+1];
1131            for (int a = _first_out[u]; a != last_out; ++a) {
1132              int v = _target[a];
1133              if (_res_cap[a] > 0 && tol.negative(_cost[a] + p - pi[v]) &&
1134                  l > level[v]) level[v] = l;
1135            }
1136          }
1137
1138          // Modify potentials
1139          double q = std::numeric_limits<double>::max();
1140          for (int u = 0; u != _res_node_num; ++u) {
1141            int lu = level[u];
1142            double p, pu = pi[u];
1143            int last_out = _first_out[u+1];
1144            for (int a = _first_out[u]; a != last_out; ++a) {
1145              if (_res_cap[a] == 0) continue;
1146              int v = _target[a];
1147              int ld = lu - level[v];
1148              if (ld > 0) {
1149                p = (_cost[a] + pu - pi[v] + epsilon) / (ld + 1);
1150                if (p < q) q = p;
1151              }
1152            }
1153          }
1154          for (int u = 0; u != _res_node_num; ++u) {
1155            pi[u] -= q * level[u];
1156          }
1157
1158          // Modify epsilon
1159          epsilon = 0;
1160          for (int u = 0; u != _res_node_num; ++u) {
1161            double curr, pu = pi[u];
1162            int last_out = _first_out[u+1];
1163            for (int a = _first_out[u]; a != last_out; ++a) {
1164              if (_res_cap[a] == 0) continue;
1165              curr = _cost[a] + pu - pi[_target[a]];
1166              if (-curr > epsilon) epsilon = -curr;
1167            }
1168          }
1169        } else {
1170          typedef HowardMmc<StaticDigraph, CostArcMap> HwMmc;
1171          typedef HartmannOrlinMmc<StaticDigraph, CostArcMap> HoMmc;
1172          typedef typename BellmanFord<StaticDigraph, CostArcMap>
1173            ::template SetDistMap<CostNodeMap>::Create BF;
1174
1175          // Set epsilon to the minimum cycle mean
1176          Cost cycle_cost = 0;
1177          int cycle_size = 1;
1178          buildResidualNetwork();
1179          HwMmc hw_mmc(_sgr, _cost_map);
1180          if (hw_mmc.findCycleMean(hw_iter_limit) == HwMmc::ITERATION_LIMIT) {
1181            // Howard's algorithm reached the iteration limit, start a
1182            // strongly polynomial algorithm instead
1183            HoMmc ho_mmc(_sgr, _cost_map);
1184            ho_mmc.findCycleMean();
1185            epsilon = -ho_mmc.cycleMean();
1186            cycle_cost = ho_mmc.cycleCost();
1187            cycle_size = ho_mmc.cycleSize();
1188          } else {
1189            // Set epsilon
1190            epsilon = -hw_mmc.cycleMean();
1191            cycle_cost = hw_mmc.cycleCost();
1192            cycle_size = hw_mmc.cycleSize();
1193          }
1194
1195          // Compute feasible potentials for the current epsilon
1196          for (int i = 0; i != int(_cost_vec.size()); ++i) {
1197            _cost_vec[i] = cycle_size * _cost_vec[i] - cycle_cost;
1198          }
1199          BF bf(_sgr, _cost_map);
1200          bf.distMap(_pi_map);
1201          bf.init(0);
1202          bf.start();
1203          for (int u = 0; u != _res_node_num; ++u) {
1204            pi[u] = static_cast<double>(_pi[u]) / cycle_size;
1205          }
1206
1207          iter = limit;
1208        }
1209      }
1210    }
1211
1212  }; //class CycleCanceling
1213
1214  ///@}
1215
1216} //namespace lemon
1217
1218#endif //LEMON_CYCLE_CANCELING_H
Note: See TracBrowser for help on using the repository browser.