COIN-OR::LEMON - Graph Library

Changeset 815:aef153f430e1 in lemon-main


Ignore:
Timestamp:
11/13/09 00:10:33 (15 years ago)
Author:
Peter Kovacs <kpeter@…>
Branch:
default
Phase:
public
Rebase:
62373361373236626665663832376334666366353631373365333236323466343130303133633939
Message:

Entirely rework cycle canceling algorithms (#180)

  • Move the cycle canceling algorithms (CycleCanceling?, CancelAndTighten?) into one class (CycleCanceling?).
  • Add a Method parameter to the run() function to be able to select the used cycle canceling method.
  • Use the new interface similarly to NetworkSimplex?.
  • Rework the implementations using an efficient internal structure for handling the residual network. This improvement made the codes much faster.
  • Handle GEQ supply type (LEQ is not supported).
  • Handle infinite upper bounds.
  • Handle negative costs (for arcs of finite upper bound).
  • Extend the documentation.
Location:
lemon
Files:
1 deleted
2 edited

Legend:

Unmodified
Added
Removed
  • lemon/Makefile.am

    r814 r815  
    6363        lemon/binom_heap.h \
    6464        lemon/bucket_heap.h \
    65         lemon/cancel_and_tighten.h \
    6665        lemon/capacity_scaling.h \
    6766        lemon/cbc.h \
  • lemon/cycle_canceling.h

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