Port cycle canceling algorithms from SVN -r3524 (#180)
authorPeter Kovacs <kpeter@inf.elte.hu>
Fri, 13 Nov 2009 00:09:35 +0100 (2009-11-12)
changeset 8140643a9c2c3ae
parent 813 25804ef35064
child 815 aef153f430e1
Port cycle canceling algorithms from SVN -r3524 (#180)
lemon/Makefile.am
lemon/cancel_and_tighten.h
lemon/cycle_canceling.h
     1.1 --- a/lemon/Makefile.am	Thu Nov 12 23:52:51 2009 +0100
     1.2 +++ b/lemon/Makefile.am	Fri Nov 13 00:09:35 2009 +0100
     1.3 @@ -62,6 +62,7 @@
     1.4  	lemon/bin_heap.h \
     1.5  	lemon/binom_heap.h \
     1.6  	lemon/bucket_heap.h \
     1.7 +	lemon/cancel_and_tighten.h \
     1.8  	lemon/capacity_scaling.h \
     1.9  	lemon/cbc.h \
    1.10  	lemon/circulation.h \
    1.11 @@ -73,6 +74,7 @@
    1.12  	lemon/cost_scaling.h \
    1.13  	lemon/counter.h \
    1.14  	lemon/cplex.h \
    1.15 +	lemon/cycle_canceling.h \
    1.16  	lemon/dfs.h \
    1.17  	lemon/dijkstra.h \
    1.18  	lemon/dim2.h \
     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/lemon/cancel_and_tighten.h	Fri Nov 13 00:09:35 2009 +0100
     2.3 @@ -0,0 +1,797 @@
     2.4 +/* -*- C++ -*-
     2.5 + *
     2.6 + * This file is a part of LEMON, a generic C++ optimization library
     2.7 + *
     2.8 + * Copyright (C) 2003-2008
     2.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    2.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    2.11 + *
    2.12 + * Permission to use, modify and distribute this software is granted
    2.13 + * provided that this copyright notice appears in all copies. For
    2.14 + * precise terms see the accompanying LICENSE file.
    2.15 + *
    2.16 + * This software is provided "AS IS" with no warranty of any kind,
    2.17 + * express or implied, and with no claim as to its suitability for any
    2.18 + * purpose.
    2.19 + *
    2.20 + */
    2.21 +
    2.22 +#ifndef LEMON_CANCEL_AND_TIGHTEN_H
    2.23 +#define LEMON_CANCEL_AND_TIGHTEN_H
    2.24 +
    2.25 +/// \ingroup min_cost_flow
    2.26 +///
    2.27 +/// \file
    2.28 +/// \brief Cancel and Tighten algorithm for finding a minimum cost flow.
    2.29 +
    2.30 +#include <vector>
    2.31 +
    2.32 +#include <lemon/circulation.h>
    2.33 +#include <lemon/bellman_ford.h>
    2.34 +#include <lemon/howard.h>
    2.35 +#include <lemon/adaptors.h>
    2.36 +#include <lemon/tolerance.h>
    2.37 +#include <lemon/math.h>
    2.38 +
    2.39 +#include <lemon/static_graph.h>
    2.40 +
    2.41 +namespace lemon {
    2.42 +
    2.43 +  /// \addtogroup min_cost_flow
    2.44 +  /// @{
    2.45 +
    2.46 +  /// \brief Implementation of the Cancel and Tighten algorithm for
    2.47 +  /// finding a minimum cost flow.
    2.48 +  ///
    2.49 +  /// \ref CancelAndTighten implements the Cancel and Tighten algorithm for
    2.50 +  /// finding a minimum cost flow.
    2.51 +  ///
    2.52 +  /// \tparam Digraph The digraph type the algorithm runs on.
    2.53 +  /// \tparam LowerMap The type of the lower bound map.
    2.54 +  /// \tparam CapacityMap The type of the capacity (upper bound) map.
    2.55 +  /// \tparam CostMap The type of the cost (length) map.
    2.56 +  /// \tparam SupplyMap The type of the supply map.
    2.57 +  ///
    2.58 +  /// \warning
    2.59 +  /// - Arc capacities and costs should be \e non-negative \e integers.
    2.60 +  /// - Supply values should be \e signed \e integers.
    2.61 +  /// - The value types of the maps should be convertible to each other.
    2.62 +  /// - \c CostMap::Value must be signed type.
    2.63 +  ///
    2.64 +  /// \author Peter Kovacs
    2.65 +  template < typename Digraph,
    2.66 +             typename LowerMap = typename Digraph::template ArcMap<int>,
    2.67 +             typename CapacityMap = typename Digraph::template ArcMap<int>,
    2.68 +             typename CostMap = typename Digraph::template ArcMap<int>,
    2.69 +             typename SupplyMap = typename Digraph::template NodeMap<int> >
    2.70 +  class CancelAndTighten
    2.71 +  {
    2.72 +    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
    2.73 +
    2.74 +    typedef typename CapacityMap::Value Capacity;
    2.75 +    typedef typename CostMap::Value Cost;
    2.76 +    typedef typename SupplyMap::Value Supply;
    2.77 +    typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
    2.78 +    typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
    2.79 +
    2.80 +    typedef ResidualDigraph< const Digraph,
    2.81 +      CapacityArcMap, CapacityArcMap > ResDigraph;
    2.82 +
    2.83 +  public:
    2.84 +
    2.85 +    /// The type of the flow map.
    2.86 +    typedef typename Digraph::template ArcMap<Capacity> FlowMap;
    2.87 +    /// The type of the potential map.
    2.88 +    typedef typename Digraph::template NodeMap<Cost> PotentialMap;
    2.89 +
    2.90 +  private:
    2.91 +
    2.92 +    /// \brief Map adaptor class for handling residual arc costs.
    2.93 +    ///
    2.94 +    /// Map adaptor class for handling residual arc costs.
    2.95 +    class ResidualCostMap : public MapBase<typename ResDigraph::Arc, Cost>
    2.96 +    {
    2.97 +      typedef typename ResDigraph::Arc Arc;
    2.98 +      
    2.99 +    private:
   2.100 +
   2.101 +      const CostMap &_cost_map;
   2.102 +
   2.103 +    public:
   2.104 +
   2.105 +      ///\e
   2.106 +      ResidualCostMap(const CostMap &cost_map) : _cost_map(cost_map) {}
   2.107 +
   2.108 +      ///\e
   2.109 +      Cost operator[](const Arc &e) const {
   2.110 +        return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
   2.111 +      }
   2.112 +
   2.113 +    }; //class ResidualCostMap
   2.114 +
   2.115 +    /// \brief Map adaptor class for handling reduced arc costs.
   2.116 +    ///
   2.117 +    /// Map adaptor class for handling reduced arc costs.
   2.118 +    class ReducedCostMap : public MapBase<Arc, Cost>
   2.119 +    {
   2.120 +    private:
   2.121 +
   2.122 +      const Digraph &_gr;
   2.123 +      const CostMap &_cost_map;
   2.124 +      const PotentialMap &_pot_map;
   2.125 +
   2.126 +    public:
   2.127 +
   2.128 +      ///\e
   2.129 +      ReducedCostMap( const Digraph &gr,
   2.130 +                      const CostMap &cost_map,
   2.131 +                      const PotentialMap &pot_map ) :
   2.132 +        _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
   2.133 +
   2.134 +      ///\e
   2.135 +      inline Cost operator[](const Arc &e) const {
   2.136 +        return _cost_map[e] + _pot_map[_gr.source(e)]
   2.137 +                            - _pot_map[_gr.target(e)];
   2.138 +      }
   2.139 +
   2.140 +    }; //class ReducedCostMap
   2.141 +
   2.142 +    struct BFOperationTraits {
   2.143 +      static double zero() { return 0; }
   2.144 +
   2.145 +      static double infinity() {
   2.146 +        return std::numeric_limits<double>::infinity();
   2.147 +      }
   2.148 +
   2.149 +      static double plus(const double& left, const double& right) {
   2.150 +        return left + right;
   2.151 +      }
   2.152 +
   2.153 +      static bool less(const double& left, const double& right) {
   2.154 +        return left + 1e-6 < right;
   2.155 +      }
   2.156 +    }; // class BFOperationTraits
   2.157 +
   2.158 +  private:
   2.159 +
   2.160 +    // The digraph the algorithm runs on
   2.161 +    const Digraph &_graph;
   2.162 +    // The original lower bound map
   2.163 +    const LowerMap *_lower;
   2.164 +    // The modified capacity map
   2.165 +    CapacityArcMap _capacity;
   2.166 +    // The original cost map
   2.167 +    const CostMap &_cost;
   2.168 +    // The modified supply map
   2.169 +    SupplyNodeMap _supply;
   2.170 +    bool _valid_supply;
   2.171 +
   2.172 +    // Arc map of the current flow
   2.173 +    FlowMap *_flow;
   2.174 +    bool _local_flow;
   2.175 +    // Node map of the current potentials
   2.176 +    PotentialMap *_potential;
   2.177 +    bool _local_potential;
   2.178 +
   2.179 +    // The residual digraph
   2.180 +    ResDigraph *_res_graph;
   2.181 +    // The residual cost map
   2.182 +    ResidualCostMap _res_cost;
   2.183 +
   2.184 +  public:
   2.185 +
   2.186 +    /// \brief General constructor (with lower bounds).
   2.187 +    ///
   2.188 +    /// General constructor (with lower bounds).
   2.189 +    ///
   2.190 +    /// \param digraph The digraph the algorithm runs on.
   2.191 +    /// \param lower The lower bounds of the arcs.
   2.192 +    /// \param capacity The capacities (upper bounds) of the arcs.
   2.193 +    /// \param cost The cost (length) values of the arcs.
   2.194 +    /// \param supply The supply values of the nodes (signed).
   2.195 +    CancelAndTighten( const Digraph &digraph,
   2.196 +                      const LowerMap &lower,
   2.197 +                      const CapacityMap &capacity,
   2.198 +                      const CostMap &cost,
   2.199 +                      const SupplyMap &supply ) :
   2.200 +      _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost),
   2.201 +      _supply(digraph), _flow(NULL), _local_flow(false),
   2.202 +      _potential(NULL), _local_potential(false),
   2.203 +      _res_graph(NULL), _res_cost(_cost)
   2.204 +    {
   2.205 +      // Check the sum of supply values
   2.206 +      Supply sum = 0;
   2.207 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   2.208 +        _supply[n] = supply[n];
   2.209 +        sum += _supply[n];
   2.210 +      }
   2.211 +      _valid_supply = sum == 0;
   2.212 +
   2.213 +      // Remove non-zero lower bounds
   2.214 +      for (ArcIt e(_graph); e != INVALID; ++e) {
   2.215 +        _capacity[e] = capacity[e];
   2.216 +        if (lower[e] != 0) {
   2.217 +          _capacity[e] -= lower[e];
   2.218 +          _supply[_graph.source(e)] -= lower[e];
   2.219 +          _supply[_graph.target(e)] += lower[e];
   2.220 +        }
   2.221 +      }
   2.222 +    }
   2.223 +/*
   2.224 +    /// \brief General constructor (without lower bounds).
   2.225 +    ///
   2.226 +    /// General constructor (without lower bounds).
   2.227 +    ///
   2.228 +    /// \param digraph The digraph the algorithm runs on.
   2.229 +    /// \param capacity The capacities (upper bounds) of the arcs.
   2.230 +    /// \param cost The cost (length) values of the arcs.
   2.231 +    /// \param supply The supply values of the nodes (signed).
   2.232 +    CancelAndTighten( const Digraph &digraph,
   2.233 +                      const CapacityMap &capacity,
   2.234 +                      const CostMap &cost,
   2.235 +                      const SupplyMap &supply ) :
   2.236 +      _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
   2.237 +      _supply(supply), _flow(NULL), _local_flow(false),
   2.238 +      _potential(NULL), _local_potential(false),
   2.239 +      _res_graph(NULL), _res_cost(_cost)
   2.240 +    {
   2.241 +      // Check the sum of supply values
   2.242 +      Supply sum = 0;
   2.243 +      for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
   2.244 +      _valid_supply = sum == 0;
   2.245 +    }
   2.246 +
   2.247 +    /// \brief Simple constructor (with lower bounds).
   2.248 +    ///
   2.249 +    /// Simple constructor (with lower bounds).
   2.250 +    ///
   2.251 +    /// \param digraph The digraph the algorithm runs on.
   2.252 +    /// \param lower The lower bounds of the arcs.
   2.253 +    /// \param capacity The capacities (upper bounds) of the arcs.
   2.254 +    /// \param cost The cost (length) values of the arcs.
   2.255 +    /// \param s The source node.
   2.256 +    /// \param t The target node.
   2.257 +    /// \param flow_value The required amount of flow from node \c s
   2.258 +    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
   2.259 +    CancelAndTighten( const Digraph &digraph,
   2.260 +                      const LowerMap &lower,
   2.261 +                      const CapacityMap &capacity,
   2.262 +                      const CostMap &cost,
   2.263 +                      Node s, Node t,
   2.264 +                      Supply flow_value ) :
   2.265 +      _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost),
   2.266 +      _supply(digraph, 0), _flow(NULL), _local_flow(false),
   2.267 +      _potential(NULL), _local_potential(false),
   2.268 +      _res_graph(NULL), _res_cost(_cost)
   2.269 +    {
   2.270 +      // Remove non-zero lower bounds
   2.271 +      _supply[s] =  flow_value;
   2.272 +      _supply[t] = -flow_value;
   2.273 +      for (ArcIt e(_graph); e != INVALID; ++e) {
   2.274 +        if (lower[e] != 0) {
   2.275 +          _capacity[e] -= lower[e];
   2.276 +          _supply[_graph.source(e)] -= lower[e];
   2.277 +          _supply[_graph.target(e)] += lower[e];
   2.278 +        }
   2.279 +      }
   2.280 +      _valid_supply = true;
   2.281 +    }
   2.282 +
   2.283 +    /// \brief Simple constructor (without lower bounds).
   2.284 +    ///
   2.285 +    /// Simple constructor (without lower bounds).
   2.286 +    ///
   2.287 +    /// \param digraph The digraph the algorithm runs on.
   2.288 +    /// \param capacity The capacities (upper bounds) of the arcs.
   2.289 +    /// \param cost The cost (length) values of the arcs.
   2.290 +    /// \param s The source node.
   2.291 +    /// \param t The target node.
   2.292 +    /// \param flow_value The required amount of flow from node \c s
   2.293 +    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
   2.294 +    CancelAndTighten( const Digraph &digraph,
   2.295 +                      const CapacityMap &capacity,
   2.296 +                      const CostMap &cost,
   2.297 +                      Node s, Node t,
   2.298 +                      Supply flow_value ) :
   2.299 +      _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
   2.300 +      _supply(digraph, 0), _flow(NULL), _local_flow(false),
   2.301 +      _potential(NULL), _local_potential(false),
   2.302 +      _res_graph(NULL), _res_cost(_cost)
   2.303 +    {
   2.304 +      _supply[s] =  flow_value;
   2.305 +      _supply[t] = -flow_value;
   2.306 +      _valid_supply = true;
   2.307 +    }
   2.308 +*/
   2.309 +    /// Destructor.
   2.310 +    ~CancelAndTighten() {
   2.311 +      if (_local_flow) delete _flow;
   2.312 +      if (_local_potential) delete _potential;
   2.313 +      delete _res_graph;
   2.314 +    }
   2.315 +
   2.316 +    /// \brief Set the flow map.
   2.317 +    ///
   2.318 +    /// Set the flow map.
   2.319 +    ///
   2.320 +    /// \return \c (*this)
   2.321 +    CancelAndTighten& flowMap(FlowMap &map) {
   2.322 +      if (_local_flow) {
   2.323 +        delete _flow;
   2.324 +        _local_flow = false;
   2.325 +      }
   2.326 +      _flow = &map;
   2.327 +      return *this;
   2.328 +    }
   2.329 +
   2.330 +    /// \brief Set the potential map.
   2.331 +    ///
   2.332 +    /// Set the potential map.
   2.333 +    ///
   2.334 +    /// \return \c (*this)
   2.335 +    CancelAndTighten& potentialMap(PotentialMap &map) {
   2.336 +      if (_local_potential) {
   2.337 +        delete _potential;
   2.338 +        _local_potential = false;
   2.339 +      }
   2.340 +      _potential = &map;
   2.341 +      return *this;
   2.342 +    }
   2.343 +
   2.344 +    /// \name Execution control
   2.345 +
   2.346 +    /// @{
   2.347 +
   2.348 +    /// \brief Run the algorithm.
   2.349 +    ///
   2.350 +    /// Run the algorithm.
   2.351 +    ///
   2.352 +    /// \return \c true if a feasible flow can be found.
   2.353 +    bool run() {
   2.354 +      return init() && start();
   2.355 +    }
   2.356 +
   2.357 +    /// @}
   2.358 +
   2.359 +    /// \name Query Functions
   2.360 +    /// The result of the algorithm can be obtained using these
   2.361 +    /// functions.\n
   2.362 +    /// \ref lemon::CancelAndTighten::run() "run()" must be called before
   2.363 +    /// using them.
   2.364 +
   2.365 +    /// @{
   2.366 +
   2.367 +    /// \brief Return a const reference to the arc map storing the
   2.368 +    /// found flow.
   2.369 +    ///
   2.370 +    /// Return a const reference to the arc map storing the found flow.
   2.371 +    ///
   2.372 +    /// \pre \ref run() must be called before using this function.
   2.373 +    const FlowMap& flowMap() const {
   2.374 +      return *_flow;
   2.375 +    }
   2.376 +
   2.377 +    /// \brief Return a const reference to the node map storing the
   2.378 +    /// found potentials (the dual solution).
   2.379 +    ///
   2.380 +    /// Return a const reference to the node map storing the found
   2.381 +    /// potentials (the dual solution).
   2.382 +    ///
   2.383 +    /// \pre \ref run() must be called before using this function.
   2.384 +    const PotentialMap& potentialMap() const {
   2.385 +      return *_potential;
   2.386 +    }
   2.387 +
   2.388 +    /// \brief Return the flow on the given arc.
   2.389 +    ///
   2.390 +    /// Return the flow on the given arc.
   2.391 +    ///
   2.392 +    /// \pre \ref run() must be called before using this function.
   2.393 +    Capacity flow(const Arc& arc) const {
   2.394 +      return (*_flow)[arc];
   2.395 +    }
   2.396 +
   2.397 +    /// \brief Return the potential of the given node.
   2.398 +    ///
   2.399 +    /// Return the potential of the given node.
   2.400 +    ///
   2.401 +    /// \pre \ref run() must be called before using this function.
   2.402 +    Cost potential(const Node& node) const {
   2.403 +      return (*_potential)[node];
   2.404 +    }
   2.405 +
   2.406 +    /// \brief Return the total cost of the found flow.
   2.407 +    ///
   2.408 +    /// Return the total cost of the found flow. The complexity of the
   2.409 +    /// function is \f$ O(e) \f$.
   2.410 +    ///
   2.411 +    /// \pre \ref run() must be called before using this function.
   2.412 +    Cost totalCost() const {
   2.413 +      Cost c = 0;
   2.414 +      for (ArcIt e(_graph); e != INVALID; ++e)
   2.415 +        c += (*_flow)[e] * _cost[e];
   2.416 +      return c;
   2.417 +    }
   2.418 +
   2.419 +    /// @}
   2.420 +
   2.421 +  private:
   2.422 +
   2.423 +    /// Initialize the algorithm.
   2.424 +    bool init() {
   2.425 +      if (!_valid_supply) return false;
   2.426 +
   2.427 +      // Initialize flow and potential maps
   2.428 +      if (!_flow) {
   2.429 +        _flow = new FlowMap(_graph);
   2.430 +        _local_flow = true;
   2.431 +      }
   2.432 +      if (!_potential) {
   2.433 +        _potential = new PotentialMap(_graph);
   2.434 +        _local_potential = true;
   2.435 +      }
   2.436 +
   2.437 +      _res_graph = new ResDigraph(_graph, _capacity, *_flow);
   2.438 +
   2.439 +      // Find a feasible flow using Circulation
   2.440 +      Circulation< Digraph, ConstMap<Arc, Capacity>,
   2.441 +                   CapacityArcMap, SupplyMap >
   2.442 +        circulation( _graph, constMap<Arc>(Capacity(0)),
   2.443 +                     _capacity, _supply );
   2.444 +      return circulation.flowMap(*_flow).run();
   2.445 +    }
   2.446 +
   2.447 +    bool start() {
   2.448 +      const double LIMIT_FACTOR = 0.01;
   2.449 +      const int MIN_LIMIT = 3;
   2.450 +
   2.451 +      typedef typename Digraph::template NodeMap<double> FloatPotentialMap;
   2.452 +      typedef typename Digraph::template NodeMap<int> LevelMap;
   2.453 +      typedef typename Digraph::template NodeMap<bool> BoolNodeMap;
   2.454 +      typedef typename Digraph::template NodeMap<Node> PredNodeMap;
   2.455 +      typedef typename Digraph::template NodeMap<Arc> PredArcMap;
   2.456 +      typedef typename ResDigraph::template ArcMap<double> ResShiftCostMap;
   2.457 +      FloatPotentialMap pi(_graph);
   2.458 +      LevelMap level(_graph);
   2.459 +      BoolNodeMap reached(_graph);
   2.460 +      BoolNodeMap processed(_graph);
   2.461 +      PredNodeMap pred_node(_graph);
   2.462 +      PredArcMap pred_arc(_graph);
   2.463 +      int node_num = countNodes(_graph);
   2.464 +      typedef std::pair<Arc, bool> pair;
   2.465 +      std::vector<pair> stack(node_num);
   2.466 +      std::vector<Node> proc_vector(node_num);
   2.467 +      ResShiftCostMap shift_cost(*_res_graph);
   2.468 +
   2.469 +      Tolerance<double> tol;
   2.470 +      tol.epsilon(1e-6);
   2.471 +
   2.472 +      Timer t1, t2, t3;
   2.473 +      t1.reset();
   2.474 +      t2.reset();
   2.475 +      t3.reset();
   2.476 +
   2.477 +      // Initialize epsilon and the node potentials
   2.478 +      double epsilon = 0;
   2.479 +      for (ArcIt e(_graph); e != INVALID; ++e) {
   2.480 +        if (_capacity[e] - (*_flow)[e] > 0 && _cost[e] < -epsilon)
   2.481 +          epsilon = -_cost[e];
   2.482 +        else if ((*_flow)[e] > 0 && _cost[e] > epsilon)
   2.483 +          epsilon = _cost[e];
   2.484 +      }
   2.485 +      for (NodeIt v(_graph); v != INVALID; ++v) {
   2.486 +        pi[v] = 0;
   2.487 +      }
   2.488 +
   2.489 +      // Start phases
   2.490 +      int limit = int(LIMIT_FACTOR * node_num);
   2.491 +      if (limit < MIN_LIMIT) limit = MIN_LIMIT;
   2.492 +      int iter = limit;
   2.493 +      while (epsilon * node_num >= 1) {
   2.494 +        t1.start();
   2.495 +        // Find and cancel cycles in the admissible digraph using DFS
   2.496 +        for (NodeIt n(_graph); n != INVALID; ++n) {
   2.497 +          reached[n] = false;
   2.498 +          processed[n] = false;
   2.499 +        }
   2.500 +        int stack_head = -1;
   2.501 +        int proc_head = -1;
   2.502 +
   2.503 +        for (NodeIt start(_graph); start != INVALID; ++start) {
   2.504 +          if (reached[start]) continue;
   2.505 +
   2.506 +          // New start node
   2.507 +          reached[start] = true;
   2.508 +          pred_arc[start] = INVALID;
   2.509 +          pred_node[start] = INVALID;
   2.510 +
   2.511 +          // Find the first admissible residual outgoing arc
   2.512 +          double p = pi[start];
   2.513 +          Arc e;
   2.514 +          _graph.firstOut(e, start);
   2.515 +          while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
   2.516 +                  !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
   2.517 +            _graph.nextOut(e);
   2.518 +          if (e != INVALID) {
   2.519 +            stack[++stack_head] = pair(e, true);
   2.520 +            goto next_step_1;
   2.521 +          }
   2.522 +          _graph.firstIn(e, start);
   2.523 +          while ( e != INVALID && ((*_flow)[e] == 0 ||
   2.524 +                  !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
   2.525 +            _graph.nextIn(e);
   2.526 +          if (e != INVALID) {
   2.527 +            stack[++stack_head] = pair(e, false);
   2.528 +            goto next_step_1;
   2.529 +          }
   2.530 +          processed[start] = true;
   2.531 +          proc_vector[++proc_head] = start;
   2.532 +          continue;
   2.533 +        next_step_1:
   2.534 +
   2.535 +          while (stack_head >= 0) {
   2.536 +            Arc se = stack[stack_head].first;
   2.537 +            bool sf = stack[stack_head].second;
   2.538 +            Node u, v;
   2.539 +            if (sf) {
   2.540 +              u = _graph.source(se);
   2.541 +              v = _graph.target(se);
   2.542 +            } else {
   2.543 +              u = _graph.target(se);
   2.544 +              v = _graph.source(se);
   2.545 +            }
   2.546 +
   2.547 +            if (!reached[v]) {
   2.548 +              // A new node is reached
   2.549 +              reached[v] = true;
   2.550 +              pred_node[v] = u;
   2.551 +              pred_arc[v] = se;
   2.552 +              // Find the first admissible residual outgoing arc
   2.553 +              double p = pi[v];
   2.554 +              Arc e;
   2.555 +              _graph.firstOut(e, v);
   2.556 +              while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
   2.557 +                      !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
   2.558 +                _graph.nextOut(e);
   2.559 +              if (e != INVALID) {
   2.560 +                stack[++stack_head] = pair(e, true);
   2.561 +                goto next_step_2;
   2.562 +              }
   2.563 +              _graph.firstIn(e, v);
   2.564 +              while ( e != INVALID && ((*_flow)[e] == 0 ||
   2.565 +                      !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
   2.566 +                _graph.nextIn(e);
   2.567 +              stack[++stack_head] = pair(e, false);
   2.568 +            next_step_2: ;
   2.569 +            } else {
   2.570 +              if (!processed[v]) {
   2.571 +                // A cycle is found
   2.572 +                Node n, w = u;
   2.573 +                Capacity d, delta = sf ? _capacity[se] - (*_flow)[se] :
   2.574 +                                         (*_flow)[se];
   2.575 +                for (n = u; n != v; n = pred_node[n]) {
   2.576 +                  d = _graph.target(pred_arc[n]) == n ?
   2.577 +                      _capacity[pred_arc[n]] - (*_flow)[pred_arc[n]] :
   2.578 +                      (*_flow)[pred_arc[n]];
   2.579 +                  if (d <= delta) {
   2.580 +                    delta = d;
   2.581 +                    w = pred_node[n];
   2.582 +                  }
   2.583 +                }
   2.584 +
   2.585 +/*
   2.586 +                std::cout << "CYCLE FOUND: ";
   2.587 +                if (sf)
   2.588 +                  std::cout << _cost[se] + pi[_graph.source(se)] - pi[_graph.target(se)];
   2.589 +                else
   2.590 +                  std::cout << _graph.id(se) << ":" << -(_cost[se] + pi[_graph.source(se)] - pi[_graph.target(se)]);
   2.591 +                for (n = u; n != v; n = pred_node[n]) {
   2.592 +                  if (_graph.target(pred_arc[n]) == n)
   2.593 +                    std::cout << " " << _cost[pred_arc[n]] + pi[_graph.source(pred_arc[n])] - pi[_graph.target(pred_arc[n])];
   2.594 +                  else
   2.595 +                    std::cout << " " << -(_cost[pred_arc[n]] + pi[_graph.source(pred_arc[n])] - pi[_graph.target(pred_arc[n])]);
   2.596 +                }
   2.597 +                std::cout << "\n";
   2.598 +*/
   2.599 +                // Augment along the cycle
   2.600 +                (*_flow)[se] = sf ? (*_flow)[se] + delta :
   2.601 +                                    (*_flow)[se] - delta;
   2.602 +                for (n = u; n != v; n = pred_node[n]) {
   2.603 +                  if (_graph.target(pred_arc[n]) == n)
   2.604 +                    (*_flow)[pred_arc[n]] += delta;
   2.605 +                  else
   2.606 +                    (*_flow)[pred_arc[n]] -= delta;
   2.607 +                }
   2.608 +                for (n = u; stack_head > 0 && n != w; n = pred_node[n]) {
   2.609 +                  --stack_head;
   2.610 +                  reached[n] = false;
   2.611 +                }
   2.612 +                u = w;
   2.613 +              }
   2.614 +              v = u;
   2.615 +
   2.616 +              // Find the next admissible residual outgoing arc
   2.617 +              double p = pi[v];
   2.618 +              Arc e = stack[stack_head].first;
   2.619 +              if (!stack[stack_head].second) {
   2.620 +                _graph.nextIn(e);
   2.621 +                goto in_arc_3;
   2.622 +              }
   2.623 +              _graph.nextOut(e);
   2.624 +              while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
   2.625 +                      !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
   2.626 +                _graph.nextOut(e);
   2.627 +              if (e != INVALID) {
   2.628 +                stack[stack_head] = pair(e, true);
   2.629 +                goto next_step_3;
   2.630 +              }
   2.631 +              _graph.firstIn(e, v);
   2.632 +            in_arc_3:
   2.633 +              while ( e != INVALID && ((*_flow)[e] == 0 ||
   2.634 +                      !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
   2.635 +                _graph.nextIn(e);
   2.636 +              stack[stack_head] = pair(e, false);
   2.637 +            next_step_3: ;
   2.638 +            }
   2.639 +
   2.640 +            while (stack_head >= 0 && stack[stack_head].first == INVALID) {
   2.641 +              processed[v] = true;
   2.642 +              proc_vector[++proc_head] = v;
   2.643 +              if (--stack_head >= 0) {
   2.644 +                v = stack[stack_head].second ?
   2.645 +                    _graph.source(stack[stack_head].first) :
   2.646 +                    _graph.target(stack[stack_head].first);
   2.647 +                // Find the next admissible residual outgoing arc
   2.648 +                double p = pi[v];
   2.649 +                Arc e = stack[stack_head].first;
   2.650 +                if (!stack[stack_head].second) {
   2.651 +                  _graph.nextIn(e);
   2.652 +                  goto in_arc_4;
   2.653 +                }
   2.654 +                _graph.nextOut(e);
   2.655 +                while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
   2.656 +                        !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
   2.657 +                  _graph.nextOut(e);
   2.658 +                if (e != INVALID) {
   2.659 +                  stack[stack_head] = pair(e, true);
   2.660 +                  goto next_step_4;
   2.661 +                }
   2.662 +                _graph.firstIn(e, v);
   2.663 +              in_arc_4:
   2.664 +                while ( e != INVALID && ((*_flow)[e] == 0 ||
   2.665 +                        !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
   2.666 +                  _graph.nextIn(e);
   2.667 +                stack[stack_head] = pair(e, false);
   2.668 +              next_step_4: ;
   2.669 +              }
   2.670 +            }
   2.671 +          }
   2.672 +        }
   2.673 +        t1.stop();
   2.674 +
   2.675 +        // Tighten potentials and epsilon
   2.676 +        if (--iter > 0) {
   2.677 +          // Compute levels
   2.678 +          t2.start();
   2.679 +          for (int i = proc_head; i >= 0; --i) {
   2.680 +            Node v = proc_vector[i];
   2.681 +            double p = pi[v];
   2.682 +            int l = 0;
   2.683 +            for (InArcIt e(_graph, v); e != INVALID; ++e) {
   2.684 +              Node u = _graph.source(e);
   2.685 +              if ( _capacity[e] - (*_flow)[e] > 0 &&
   2.686 +                   tol.negative(_cost[e] + pi[u] - p) &&
   2.687 +                   level[u] + 1 > l ) l = level[u] + 1;
   2.688 +            }
   2.689 +            for (OutArcIt e(_graph, v); e != INVALID; ++e) {
   2.690 +              Node u = _graph.target(e);
   2.691 +              if ( (*_flow)[e] > 0 &&
   2.692 +                   tol.negative(-_cost[e] + pi[u] - p) &&
   2.693 +                   level[u] + 1 > l ) l = level[u] + 1;
   2.694 +            }
   2.695 +            level[v] = l;
   2.696 +          }
   2.697 +
   2.698 +          // Modify potentials
   2.699 +          double p, q = -1;
   2.700 +          for (ArcIt e(_graph); e != INVALID; ++e) {
   2.701 +            Node u = _graph.source(e);
   2.702 +            Node v = _graph.target(e);
   2.703 +            if (_capacity[e] - (*_flow)[e] > 0 && level[u] - level[v] > 0) {
   2.704 +              p = (_cost[e] + pi[u] - pi[v] + epsilon) /
   2.705 +                  (level[u] - level[v] + 1);
   2.706 +              if (q < 0 || p < q) q = p;
   2.707 +            }
   2.708 +            else if ((*_flow)[e] > 0 && level[v] - level[u] > 0) {
   2.709 +              p = (-_cost[e] - pi[u] + pi[v] + epsilon) /
   2.710 +                  (level[v] - level[u] + 1);
   2.711 +              if (q < 0 || p < q) q = p;
   2.712 +            }
   2.713 +          }
   2.714 +          for (NodeIt v(_graph); v != INVALID; ++v) {
   2.715 +            pi[v] -= q * level[v];
   2.716 +          }
   2.717 +
   2.718 +          // Modify epsilon
   2.719 +          epsilon = 0;
   2.720 +          for (ArcIt e(_graph); e != INVALID; ++e) {
   2.721 +            double curr = _cost[e] + pi[_graph.source(e)]
   2.722 +                                   - pi[_graph.target(e)];
   2.723 +            if (_capacity[e] - (*_flow)[e] > 0 && curr < -epsilon)
   2.724 +              epsilon = -curr;
   2.725 +            else if ((*_flow)[e] > 0 && curr > epsilon)
   2.726 +              epsilon = curr;
   2.727 +          }
   2.728 +          t2.stop();
   2.729 +        } else {
   2.730 +          // Set epsilon to the minimum cycle mean
   2.731 +          t3.start();
   2.732 +
   2.733 +/**/
   2.734 +          StaticDigraph static_graph;
   2.735 +          typename ResDigraph::template NodeMap<typename StaticDigraph::Node> node_ref(*_res_graph);
   2.736 +          typename ResDigraph::template ArcMap<typename StaticDigraph::Arc> arc_ref(*_res_graph);
   2.737 +          static_graph.build(*_res_graph, node_ref, arc_ref);
   2.738 +          typename StaticDigraph::template NodeMap<double> static_pi(static_graph);
   2.739 +          typename StaticDigraph::template ArcMap<double> static_cost(static_graph);
   2.740 +
   2.741 +          for (typename ResDigraph::ArcIt e(*_res_graph); e != INVALID; ++e)
   2.742 +            static_cost[arc_ref[e]] = _res_cost[e];
   2.743 +
   2.744 +          Howard<StaticDigraph, typename StaticDigraph::template ArcMap<double> >
   2.745 +            mmc(static_graph, static_cost);
   2.746 +          mmc.findMinMean();
   2.747 +          epsilon = -mmc.cycleMean();
   2.748 +/**/
   2.749 +
   2.750 +/*
   2.751 +          Howard<ResDigraph, ResidualCostMap> mmc(*_res_graph, _res_cost);
   2.752 +          mmc.findMinMean();
   2.753 +          epsilon = -mmc.cycleMean();
   2.754 +*/
   2.755 +
   2.756 +          // Compute feasible potentials for the current epsilon
   2.757 +          for (typename StaticDigraph::ArcIt e(static_graph); e != INVALID; ++e)
   2.758 +            static_cost[e] += epsilon;
   2.759 +          typename BellmanFord<StaticDigraph, typename StaticDigraph::template ArcMap<double> >::
   2.760 +            template SetDistMap<typename StaticDigraph::template NodeMap<double> >::
   2.761 +            template SetOperationTraits<BFOperationTraits>::Create
   2.762 +              bf(static_graph, static_cost);
   2.763 +          bf.distMap(static_pi).init(0);
   2.764 +          bf.start();
   2.765 +          for (NodeIt n(_graph); n != INVALID; ++n)
   2.766 +            pi[n] = static_pi[node_ref[n]];
   2.767 +          
   2.768 +/*
   2.769 +          for (typename ResDigraph::ArcIt e(*_res_graph); e != INVALID; ++e)
   2.770 +            shift_cost[e] = _res_cost[e] + epsilon;
   2.771 +          typename BellmanFord<ResDigraph, ResShiftCostMap>::
   2.772 +            template SetDistMap<FloatPotentialMap>::
   2.773 +            template SetOperationTraits<BFOperationTraits>::Create
   2.774 +              bf(*_res_graph, shift_cost);
   2.775 +          bf.distMap(pi).init(0);
   2.776 +          bf.start();
   2.777 +*/
   2.778 +
   2.779 +          iter = limit;
   2.780 +          t3.stop();
   2.781 +        }
   2.782 +      }
   2.783 +
   2.784 +//      std::cout << t1.realTime() << " " << t2.realTime() << " " << t3.realTime() << "\n";
   2.785 +
   2.786 +      // Handle non-zero lower bounds
   2.787 +      if (_lower) {
   2.788 +        for (ArcIt e(_graph); e != INVALID; ++e)
   2.789 +          (*_flow)[e] += (*_lower)[e];
   2.790 +      }
   2.791 +      return true;
   2.792 +    }
   2.793 +
   2.794 +  }; //class CancelAndTighten
   2.795 +
   2.796 +  ///@}
   2.797 +
   2.798 +} //namespace lemon
   2.799 +
   2.800 +#endif //LEMON_CANCEL_AND_TIGHTEN_H
     3.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     3.2 +++ b/lemon/cycle_canceling.h	Fri Nov 13 00:09:35 2009 +0100
     3.3 @@ -0,0 +1,559 @@
     3.4 +/* -*- C++ -*-
     3.5 + *
     3.6 + * This file is a part of LEMON, a generic C++ optimization library
     3.7 + *
     3.8 + * Copyright (C) 2003-2008
     3.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    3.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    3.11 + *
    3.12 + * Permission to use, modify and distribute this software is granted
    3.13 + * provided that this copyright notice appears in all copies. For
    3.14 + * precise terms see the accompanying LICENSE file.
    3.15 + *
    3.16 + * This software is provided "AS IS" with no warranty of any kind,
    3.17 + * express or implied, and with no claim as to its suitability for any
    3.18 + * purpose.
    3.19 + *
    3.20 + */
    3.21 +
    3.22 +#ifndef LEMON_CYCLE_CANCELING_H
    3.23 +#define LEMON_CYCLE_CANCELING_H
    3.24 +
    3.25 +/// \ingroup min_cost_flow
    3.26 +///
    3.27 +/// \file
    3.28 +/// \brief Cycle-canceling algorithm for finding a minimum cost flow.
    3.29 +
    3.30 +#include <vector>
    3.31 +#include <lemon/adaptors.h>
    3.32 +#include <lemon/path.h>
    3.33 +
    3.34 +#include <lemon/circulation.h>
    3.35 +#include <lemon/bellman_ford.h>
    3.36 +#include <lemon/howard.h>
    3.37 +
    3.38 +namespace lemon {
    3.39 +
    3.40 +  /// \addtogroup min_cost_flow
    3.41 +  /// @{
    3.42 +
    3.43 +  /// \brief Implementation of a cycle-canceling algorithm for
    3.44 +  /// finding a minimum cost flow.
    3.45 +  ///
    3.46 +  /// \ref CycleCanceling implements a cycle-canceling algorithm for
    3.47 +  /// finding a minimum cost flow.
    3.48 +  ///
    3.49 +  /// \tparam Digraph The digraph type the algorithm runs on.
    3.50 +  /// \tparam LowerMap The type of the lower bound map.
    3.51 +  /// \tparam CapacityMap The type of the capacity (upper bound) map.
    3.52 +  /// \tparam CostMap The type of the cost (length) map.
    3.53 +  /// \tparam SupplyMap The type of the supply map.
    3.54 +  ///
    3.55 +  /// \warning
    3.56 +  /// - Arc capacities and costs should be \e non-negative \e integers.
    3.57 +  /// - Supply values should be \e signed \e integers.
    3.58 +  /// - The value types of the maps should be convertible to each other.
    3.59 +  /// - \c CostMap::Value must be signed type.
    3.60 +  ///
    3.61 +  /// \note By default the \ref BellmanFord "Bellman-Ford" algorithm is
    3.62 +  /// used for negative cycle detection with limited iteration number.
    3.63 +  /// However \ref CycleCanceling also provides the "Minimum Mean
    3.64 +  /// Cycle-Canceling" algorithm, which is \e strongly \e polynomial,
    3.65 +  /// but rather slower in practice.
    3.66 +  /// To use this version of the algorithm, call \ref run() with \c true
    3.67 +  /// parameter.
    3.68 +  ///
    3.69 +  /// \author Peter Kovacs
    3.70 +  template < typename Digraph,
    3.71 +             typename LowerMap = typename Digraph::template ArcMap<int>,
    3.72 +             typename CapacityMap = typename Digraph::template ArcMap<int>,
    3.73 +             typename CostMap = typename Digraph::template ArcMap<int>,
    3.74 +             typename SupplyMap = typename Digraph::template NodeMap<int> >
    3.75 +  class CycleCanceling
    3.76 +  {
    3.77 +    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
    3.78 +
    3.79 +    typedef typename CapacityMap::Value Capacity;
    3.80 +    typedef typename CostMap::Value Cost;
    3.81 +    typedef typename SupplyMap::Value Supply;
    3.82 +    typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
    3.83 +    typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
    3.84 +
    3.85 +    typedef ResidualDigraph< const Digraph,
    3.86 +      CapacityArcMap, CapacityArcMap > ResDigraph;
    3.87 +    typedef typename ResDigraph::Node ResNode;
    3.88 +    typedef typename ResDigraph::NodeIt ResNodeIt;
    3.89 +    typedef typename ResDigraph::Arc ResArc;
    3.90 +    typedef typename ResDigraph::ArcIt ResArcIt;
    3.91 +
    3.92 +  public:
    3.93 +
    3.94 +    /// The type of the flow map.
    3.95 +    typedef typename Digraph::template ArcMap<Capacity> FlowMap;
    3.96 +    /// The type of the potential map.
    3.97 +    typedef typename Digraph::template NodeMap<Cost> PotentialMap;
    3.98 +
    3.99 +  private:
   3.100 +
   3.101 +    /// \brief Map adaptor class for handling residual arc costs.
   3.102 +    ///
   3.103 +    /// Map adaptor class for handling residual arc costs.
   3.104 +    class ResidualCostMap : public MapBase<ResArc, Cost>
   3.105 +    {
   3.106 +    private:
   3.107 +
   3.108 +      const CostMap &_cost_map;
   3.109 +
   3.110 +    public:
   3.111 +
   3.112 +      ///\e
   3.113 +      ResidualCostMap(const CostMap &cost_map) : _cost_map(cost_map) {}
   3.114 +
   3.115 +      ///\e
   3.116 +      Cost operator[](const ResArc &e) const {
   3.117 +        return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
   3.118 +      }
   3.119 +
   3.120 +    }; //class ResidualCostMap
   3.121 +
   3.122 +  private:
   3.123 +
   3.124 +    // The maximum number of iterations for the first execution of the
   3.125 +    // Bellman-Ford algorithm. It should be at least 2.
   3.126 +    static const int BF_FIRST_LIMIT  = 2;
   3.127 +    // The iteration limit for the Bellman-Ford algorithm is multiplied
   3.128 +    // by BF_LIMIT_FACTOR/100 in every round.
   3.129 +    static const int BF_LIMIT_FACTOR = 150;
   3.130 +
   3.131 +  private:
   3.132 +
   3.133 +    // The digraph the algorithm runs on
   3.134 +    const Digraph &_graph;
   3.135 +    // The original lower bound map
   3.136 +    const LowerMap *_lower;
   3.137 +    // The modified capacity map
   3.138 +    CapacityArcMap _capacity;
   3.139 +    // The original cost map
   3.140 +    const CostMap &_cost;
   3.141 +    // The modified supply map
   3.142 +    SupplyNodeMap _supply;
   3.143 +    bool _valid_supply;
   3.144 +
   3.145 +    // Arc map of the current flow
   3.146 +    FlowMap *_flow;
   3.147 +    bool _local_flow;
   3.148 +    // Node map of the current potentials
   3.149 +    PotentialMap *_potential;
   3.150 +    bool _local_potential;
   3.151 +
   3.152 +    // The residual digraph
   3.153 +    ResDigraph *_res_graph;
   3.154 +    // The residual cost map
   3.155 +    ResidualCostMap _res_cost;
   3.156 +
   3.157 +  public:
   3.158 +
   3.159 +    /// \brief General constructor (with lower bounds).
   3.160 +    ///
   3.161 +    /// General constructor (with lower bounds).
   3.162 +    ///
   3.163 +    /// \param digraph The digraph the algorithm runs on.
   3.164 +    /// \param lower The lower bounds of the arcs.
   3.165 +    /// \param capacity The capacities (upper bounds) of the arcs.
   3.166 +    /// \param cost The cost (length) values of the arcs.
   3.167 +    /// \param supply The supply values of the nodes (signed).
   3.168 +    CycleCanceling( const Digraph &digraph,
   3.169 +                    const LowerMap &lower,
   3.170 +                    const CapacityMap &capacity,
   3.171 +                    const CostMap &cost,
   3.172 +                    const SupplyMap &supply ) :
   3.173 +      _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost),
   3.174 +      _supply(digraph), _flow(NULL), _local_flow(false),
   3.175 +      _potential(NULL), _local_potential(false),
   3.176 +      _res_graph(NULL), _res_cost(_cost)
   3.177 +    {
   3.178 +      // Check the sum of supply values
   3.179 +      Supply sum = 0;
   3.180 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   3.181 +        _supply[n] = supply[n];
   3.182 +        sum += _supply[n];
   3.183 +      }
   3.184 +      _valid_supply = sum == 0;
   3.185 +
   3.186 +      // Remove non-zero lower bounds
   3.187 +      for (ArcIt e(_graph); e != INVALID; ++e) {
   3.188 +        _capacity[e] = capacity[e];
   3.189 +        if (lower[e] != 0) {
   3.190 +          _capacity[e] -= lower[e];
   3.191 +          _supply[_graph.source(e)] -= lower[e];
   3.192 +          _supply[_graph.target(e)] += lower[e];
   3.193 +        }
   3.194 +      }
   3.195 +    }
   3.196 +/*
   3.197 +    /// \brief General constructor (without lower bounds).
   3.198 +    ///
   3.199 +    /// General constructor (without lower bounds).
   3.200 +    ///
   3.201 +    /// \param digraph The digraph the algorithm runs on.
   3.202 +    /// \param capacity The capacities (upper bounds) of the arcs.
   3.203 +    /// \param cost The cost (length) values of the arcs.
   3.204 +    /// \param supply The supply values of the nodes (signed).
   3.205 +    CycleCanceling( const Digraph &digraph,
   3.206 +                    const CapacityMap &capacity,
   3.207 +                    const CostMap &cost,
   3.208 +                    const SupplyMap &supply ) :
   3.209 +      _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
   3.210 +      _supply(supply), _flow(NULL), _local_flow(false),
   3.211 +      _potential(NULL), _local_potential(false), _res_graph(NULL),
   3.212 +      _res_cost(_cost)
   3.213 +    {
   3.214 +      // Check the sum of supply values
   3.215 +      Supply sum = 0;
   3.216 +      for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
   3.217 +      _valid_supply = sum == 0;
   3.218 +    }
   3.219 +
   3.220 +    /// \brief Simple constructor (with lower bounds).
   3.221 +    ///
   3.222 +    /// Simple constructor (with lower bounds).
   3.223 +    ///
   3.224 +    /// \param digraph The digraph the algorithm runs on.
   3.225 +    /// \param lower The lower bounds of the arcs.
   3.226 +    /// \param capacity The capacities (upper bounds) of the arcs.
   3.227 +    /// \param cost The cost (length) values of the arcs.
   3.228 +    /// \param s The source node.
   3.229 +    /// \param t The target node.
   3.230 +    /// \param flow_value The required amount of flow from node \c s
   3.231 +    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
   3.232 +    CycleCanceling( const Digraph &digraph,
   3.233 +                    const LowerMap &lower,
   3.234 +                    const CapacityMap &capacity,
   3.235 +                    const CostMap &cost,
   3.236 +                    Node s, Node t,
   3.237 +                    Supply flow_value ) :
   3.238 +      _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost),
   3.239 +      _supply(digraph, 0), _flow(NULL), _local_flow(false),
   3.240 +      _potential(NULL), _local_potential(false), _res_graph(NULL),
   3.241 +      _res_cost(_cost)
   3.242 +    {
   3.243 +      // Remove non-zero lower bounds
   3.244 +      _supply[s] =  flow_value;
   3.245 +      _supply[t] = -flow_value;
   3.246 +      for (ArcIt e(_graph); e != INVALID; ++e) {
   3.247 +        if (lower[e] != 0) {
   3.248 +          _capacity[e] -= lower[e];
   3.249 +          _supply[_graph.source(e)] -= lower[e];
   3.250 +          _supply[_graph.target(e)] += lower[e];
   3.251 +        }
   3.252 +      }
   3.253 +      _valid_supply = true;
   3.254 +    }
   3.255 +
   3.256 +    /// \brief Simple constructor (without lower bounds).
   3.257 +    ///
   3.258 +    /// Simple constructor (without lower bounds).
   3.259 +    ///
   3.260 +    /// \param digraph The digraph the algorithm runs on.
   3.261 +    /// \param capacity The capacities (upper bounds) of the arcs.
   3.262 +    /// \param cost The cost (length) values of the arcs.
   3.263 +    /// \param s The source node.
   3.264 +    /// \param t The target node.
   3.265 +    /// \param flow_value The required amount of flow from node \c s
   3.266 +    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
   3.267 +    CycleCanceling( const Digraph &digraph,
   3.268 +                    const CapacityMap &capacity,
   3.269 +                    const CostMap &cost,
   3.270 +                    Node s, Node t,
   3.271 +                    Supply flow_value ) :
   3.272 +      _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
   3.273 +      _supply(digraph, 0), _flow(NULL), _local_flow(false),
   3.274 +      _potential(NULL), _local_potential(false), _res_graph(NULL),
   3.275 +      _res_cost(_cost)
   3.276 +    {
   3.277 +      _supply[s] =  flow_value;
   3.278 +      _supply[t] = -flow_value;
   3.279 +      _valid_supply = true;
   3.280 +    }
   3.281 +*/
   3.282 +    /// Destructor.
   3.283 +    ~CycleCanceling() {
   3.284 +      if (_local_flow) delete _flow;
   3.285 +      if (_local_potential) delete _potential;
   3.286 +      delete _res_graph;
   3.287 +    }
   3.288 +
   3.289 +    /// \brief Set the flow map.
   3.290 +    ///
   3.291 +    /// Set the flow map.
   3.292 +    ///
   3.293 +    /// \return \c (*this)
   3.294 +    CycleCanceling& flowMap(FlowMap &map) {
   3.295 +      if (_local_flow) {
   3.296 +        delete _flow;
   3.297 +        _local_flow = false;
   3.298 +      }
   3.299 +      _flow = &map;
   3.300 +      return *this;
   3.301 +    }
   3.302 +
   3.303 +    /// \brief Set the potential map.
   3.304 +    ///
   3.305 +    /// Set the potential map.
   3.306 +    ///
   3.307 +    /// \return \c (*this)
   3.308 +    CycleCanceling& potentialMap(PotentialMap &map) {
   3.309 +      if (_local_potential) {
   3.310 +        delete _potential;
   3.311 +        _local_potential = false;
   3.312 +      }
   3.313 +      _potential = &map;
   3.314 +      return *this;
   3.315 +    }
   3.316 +
   3.317 +    /// \name Execution control
   3.318 +
   3.319 +    /// @{
   3.320 +
   3.321 +    /// \brief Run the algorithm.
   3.322 +    ///
   3.323 +    /// Run the algorithm.
   3.324 +    ///
   3.325 +    /// \param min_mean_cc Set this parameter to \c true to run the
   3.326 +    /// "Minimum Mean Cycle-Canceling" algorithm, which is strongly
   3.327 +    /// polynomial, but rather slower in practice.
   3.328 +    ///
   3.329 +    /// \return \c true if a feasible flow can be found.
   3.330 +    bool run(bool min_mean_cc = false) {
   3.331 +      return init() && start(min_mean_cc);
   3.332 +    }
   3.333 +
   3.334 +    /// @}
   3.335 +
   3.336 +    /// \name Query Functions
   3.337 +    /// The result of the algorithm can be obtained using these
   3.338 +    /// functions.\n
   3.339 +    /// \ref lemon::CycleCanceling::run() "run()" must be called before
   3.340 +    /// using them.
   3.341 +
   3.342 +    /// @{
   3.343 +
   3.344 +    /// \brief Return a const reference to the arc map storing the
   3.345 +    /// found flow.
   3.346 +    ///
   3.347 +    /// Return a const reference to the arc map storing the found flow.
   3.348 +    ///
   3.349 +    /// \pre \ref run() must be called before using this function.
   3.350 +    const FlowMap& flowMap() const {
   3.351 +      return *_flow;
   3.352 +    }
   3.353 +
   3.354 +    /// \brief Return a const reference to the node map storing the
   3.355 +    /// found potentials (the dual solution).
   3.356 +    ///
   3.357 +    /// Return a const reference to the node map storing the found
   3.358 +    /// potentials (the dual solution).
   3.359 +    ///
   3.360 +    /// \pre \ref run() must be called before using this function.
   3.361 +    const PotentialMap& potentialMap() const {
   3.362 +      return *_potential;
   3.363 +    }
   3.364 +
   3.365 +    /// \brief Return the flow on the given arc.
   3.366 +    ///
   3.367 +    /// Return the flow on the given arc.
   3.368 +    ///
   3.369 +    /// \pre \ref run() must be called before using this function.
   3.370 +    Capacity flow(const Arc& arc) const {
   3.371 +      return (*_flow)[arc];
   3.372 +    }
   3.373 +
   3.374 +    /// \brief Return the potential of the given node.
   3.375 +    ///
   3.376 +    /// Return the potential of the given node.
   3.377 +    ///
   3.378 +    /// \pre \ref run() must be called before using this function.
   3.379 +    Cost potential(const Node& node) const {
   3.380 +      return (*_potential)[node];
   3.381 +    }
   3.382 +
   3.383 +    /// \brief Return the total cost of the found flow.
   3.384 +    ///
   3.385 +    /// Return the total cost of the found flow. The complexity of the
   3.386 +    /// function is \f$ O(e) \f$.
   3.387 +    ///
   3.388 +    /// \pre \ref run() must be called before using this function.
   3.389 +    Cost totalCost() const {
   3.390 +      Cost c = 0;
   3.391 +      for (ArcIt e(_graph); e != INVALID; ++e)
   3.392 +        c += (*_flow)[e] * _cost[e];
   3.393 +      return c;
   3.394 +    }
   3.395 +
   3.396 +    /// @}
   3.397 +
   3.398 +  private:
   3.399 +
   3.400 +    /// Initialize the algorithm.
   3.401 +    bool init() {
   3.402 +      if (!_valid_supply) return false;
   3.403 +
   3.404 +      // Initializing flow and potential maps
   3.405 +      if (!_flow) {
   3.406 +        _flow = new FlowMap(_graph);
   3.407 +        _local_flow = true;
   3.408 +      }
   3.409 +      if (!_potential) {
   3.410 +        _potential = new PotentialMap(_graph);
   3.411 +        _local_potential = true;
   3.412 +      }
   3.413 +
   3.414 +      _res_graph = new ResDigraph(_graph, _capacity, *_flow);
   3.415 +
   3.416 +      // Finding a feasible flow using Circulation
   3.417 +      Circulation< Digraph, ConstMap<Arc, Capacity>, CapacityArcMap,
   3.418 +                   SupplyMap >
   3.419 +        circulation( _graph, constMap<Arc>(Capacity(0)), _capacity,
   3.420 +                     _supply );
   3.421 +      return circulation.flowMap(*_flow).run();
   3.422 +    }
   3.423 +
   3.424 +    bool start(bool min_mean_cc) {
   3.425 +      if (min_mean_cc)
   3.426 +        startMinMean();
   3.427 +      else
   3.428 +        start();
   3.429 +
   3.430 +      // Handling non-zero lower bounds
   3.431 +      if (_lower) {
   3.432 +        for (ArcIt e(_graph); e != INVALID; ++e)
   3.433 +          (*_flow)[e] += (*_lower)[e];
   3.434 +      }
   3.435 +      return true;
   3.436 +    }
   3.437 +
   3.438 +    /// \brief Execute the algorithm using \ref BellmanFord.
   3.439 +    ///
   3.440 +    /// Execute the algorithm using the \ref BellmanFord
   3.441 +    /// "Bellman-Ford" algorithm for negative cycle detection with
   3.442 +    /// successively larger limit for the number of iterations.
   3.443 +    void start() {
   3.444 +      typename BellmanFord<ResDigraph, ResidualCostMap>::PredMap pred(*_res_graph);
   3.445 +      typename ResDigraph::template NodeMap<int> visited(*_res_graph);
   3.446 +      std::vector<ResArc> cycle;
   3.447 +      int node_num = countNodes(_graph);
   3.448 +
   3.449 +      int length_bound = BF_FIRST_LIMIT;
   3.450 +      bool optimal = false;
   3.451 +      while (!optimal) {
   3.452 +        BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost);
   3.453 +        bf.predMap(pred);
   3.454 +        bf.init(0);
   3.455 +        int iter_num = 0;
   3.456 +        bool cycle_found = false;
   3.457 +        while (!cycle_found) {
   3.458 +          int curr_iter_num = iter_num + length_bound <= node_num ?
   3.459 +                              length_bound : node_num - iter_num;
   3.460 +          iter_num += curr_iter_num;
   3.461 +          int real_iter_num = curr_iter_num;
   3.462 +          for (int i = 0; i < curr_iter_num; ++i) {
   3.463 +            if (bf.processNextWeakRound()) {
   3.464 +              real_iter_num = i;
   3.465 +              break;
   3.466 +            }
   3.467 +          }
   3.468 +          if (real_iter_num < curr_iter_num) {
   3.469 +            // Optimal flow is found
   3.470 +            optimal = true;
   3.471 +            // Setting node potentials
   3.472 +            for (NodeIt n(_graph); n != INVALID; ++n)
   3.473 +              (*_potential)[n] = bf.dist(n);
   3.474 +            break;
   3.475 +          } else {
   3.476 +            // Searching for node disjoint negative cycles
   3.477 +            for (ResNodeIt n(*_res_graph); n != INVALID; ++n)
   3.478 +              visited[n] = 0;
   3.479 +            int id = 0;
   3.480 +            for (ResNodeIt n(*_res_graph); n != INVALID; ++n) {
   3.481 +              if (visited[n] > 0) continue;
   3.482 +              visited[n] = ++id;
   3.483 +              ResNode u = pred[n] == INVALID ?
   3.484 +                          INVALID : _res_graph->source(pred[n]);
   3.485 +              while (u != INVALID && visited[u] == 0) {
   3.486 +                visited[u] = id;
   3.487 +                u = pred[u] == INVALID ?
   3.488 +                    INVALID : _res_graph->source(pred[u]);
   3.489 +              }
   3.490 +              if (u != INVALID && visited[u] == id) {
   3.491 +                // Finding the negative cycle
   3.492 +                cycle_found = true;
   3.493 +                cycle.clear();
   3.494 +                ResArc e = pred[u];
   3.495 +                cycle.push_back(e);
   3.496 +                Capacity d = _res_graph->residualCapacity(e);
   3.497 +                while (_res_graph->source(e) != u) {
   3.498 +                  cycle.push_back(e = pred[_res_graph->source(e)]);
   3.499 +                  if (_res_graph->residualCapacity(e) < d)
   3.500 +                    d = _res_graph->residualCapacity(e);
   3.501 +                }
   3.502 +
   3.503 +                // Augmenting along the cycle
   3.504 +                for (int i = 0; i < int(cycle.size()); ++i)
   3.505 +                  _res_graph->augment(cycle[i], d);
   3.506 +              }
   3.507 +            }
   3.508 +          }
   3.509 +
   3.510 +          if (!cycle_found)
   3.511 +            length_bound = length_bound * BF_LIMIT_FACTOR / 100;
   3.512 +        }
   3.513 +      }
   3.514 +    }
   3.515 +
   3.516 +    /// \brief Execute the algorithm using \ref Howard.
   3.517 +    ///
   3.518 +    /// Execute the algorithm using \ref Howard for negative
   3.519 +    /// cycle detection.
   3.520 +    void startMinMean() {
   3.521 +      typedef Path<ResDigraph> ResPath;
   3.522 +      Howard<ResDigraph, ResidualCostMap> mmc(*_res_graph, _res_cost);
   3.523 +      ResPath cycle;
   3.524 +
   3.525 +      mmc.cycle(cycle);
   3.526 +      if (mmc.findMinMean()) {
   3.527 +        while (mmc.cycleLength() < 0) {
   3.528 +          // Finding the cycle
   3.529 +          mmc.findCycle();
   3.530 +
   3.531 +          // Finding the largest flow amount that can be augmented
   3.532 +          // along the cycle
   3.533 +          Capacity delta = 0;
   3.534 +          for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e) {
   3.535 +            if (delta == 0 || _res_graph->residualCapacity(e) < delta)
   3.536 +              delta = _res_graph->residualCapacity(e);
   3.537 +          }
   3.538 +
   3.539 +          // Augmenting along the cycle
   3.540 +          for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e)
   3.541 +            _res_graph->augment(e, delta);
   3.542 +
   3.543 +          // Finding the minimum cycle mean for the modified residual
   3.544 +          // digraph
   3.545 +          if (!mmc.findMinMean()) break;
   3.546 +        }
   3.547 +      }
   3.548 +
   3.549 +      // Computing node potentials
   3.550 +      BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost);
   3.551 +      bf.init(0); bf.start();
   3.552 +      for (NodeIt n(_graph); n != INVALID; ++n)
   3.553 +        (*_potential)[n] = bf.dist(n);
   3.554 +    }
   3.555 +
   3.556 +  }; //class CycleCanceling
   3.557 +
   3.558 +  ///@}
   3.559 +
   3.560 +} //namespace lemon
   3.561 +
   3.562 +#endif //LEMON_CYCLE_CANCELING_H