COIN-OR::LEMON - Graph Library

source: lemon/lemon/cycle_canceling.h @ 942:d3ea191c3412

Last change on this file since 942:d3ea191c3412 was 942:d3ea191c3412, checked in by Peter Kovacs <kpeter@…>, 15 years ago

Rename min mean cycle classes and their members (#179)
with respect to the possible introduction of min ratio
cycle algorithms in the future.

The renamed classes:

The renamed members:

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