1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/cycle_canceling.h Sun Aug 11 15:28:12 2013 +0200
1.3 @@ -0,0 +1,1170 @@
1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library.
1.7 + *
1.8 + * Copyright (C) 2003-2010
1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 + *
1.12 + * Permission to use, modify and distribute this software is granted
1.13 + * provided that this copyright notice appears in all copies. For
1.14 + * precise terms see the accompanying LICENSE file.
1.15 + *
1.16 + * This software is provided "AS IS" with no warranty of any kind,
1.17 + * express or implied, and with no claim as to its suitability for any
1.18 + * purpose.
1.19 + *
1.20 + */
1.21 +
1.22 +#ifndef LEMON_CYCLE_CANCELING_H
1.23 +#define LEMON_CYCLE_CANCELING_H
1.24 +
1.25 +/// \ingroup min_cost_flow_algs
1.26 +/// \file
1.27 +/// \brief Cycle-canceling algorithms for finding a minimum cost flow.
1.28 +
1.29 +#include <vector>
1.30 +#include <limits>
1.31 +
1.32 +#include <lemon/core.h>
1.33 +#include <lemon/maps.h>
1.34 +#include <lemon/path.h>
1.35 +#include <lemon/math.h>
1.36 +#include <lemon/static_graph.h>
1.37 +#include <lemon/adaptors.h>
1.38 +#include <lemon/circulation.h>
1.39 +#include <lemon/bellman_ford.h>
1.40 +#include <lemon/howard_mmc.h>
1.41 +
1.42 +namespace lemon {
1.43 +
1.44 + /// \addtogroup min_cost_flow_algs
1.45 + /// @{
1.46 +
1.47 + /// \brief Implementation of cycle-canceling algorithms for
1.48 + /// finding a \ref min_cost_flow "minimum cost flow".
1.49 + ///
1.50 + /// \ref CycleCanceling implements three different cycle-canceling
1.51 + /// algorithms for finding a \ref min_cost_flow "minimum cost flow"
1.52 + /// \ref amo93networkflows, \ref klein67primal,
1.53 + /// \ref goldberg89cyclecanceling.
1.54 + /// The most efficent one (both theoretically and practically)
1.55 + /// is the \ref CANCEL_AND_TIGHTEN "Cancel and Tighten" algorithm,
1.56 + /// thus it is the default method.
1.57 + /// It is strongly polynomial, but in practice, it is typically much
1.58 + /// slower than the scaling algorithms and NetworkSimplex.
1.59 + ///
1.60 + /// Most of the parameters of the problem (except for the digraph)
1.61 + /// can be given using separate functions, and the algorithm can be
1.62 + /// executed using the \ref run() function. If some parameters are not
1.63 + /// specified, then default values will be used.
1.64 + ///
1.65 + /// \tparam GR The digraph type the algorithm runs on.
1.66 + /// \tparam V The number type used for flow amounts, capacity bounds
1.67 + /// and supply values in the algorithm. By default, it is \c int.
1.68 + /// \tparam C The number type used for costs and potentials in the
1.69 + /// algorithm. By default, it is the same as \c V.
1.70 + ///
1.71 + /// \warning Both number types must be signed and all input data must
1.72 + /// be integer.
1.73 + /// \warning This algorithm does not support negative costs for such
1.74 + /// arcs that have infinite upper bound.
1.75 + ///
1.76 + /// \note For more information about the three available methods,
1.77 + /// see \ref Method.
1.78 +#ifdef DOXYGEN
1.79 + template <typename GR, typename V, typename C>
1.80 +#else
1.81 + template <typename GR, typename V = int, typename C = V>
1.82 +#endif
1.83 + class CycleCanceling
1.84 + {
1.85 + public:
1.86 +
1.87 + /// The type of the digraph
1.88 + typedef GR Digraph;
1.89 + /// The type of the flow amounts, capacity bounds and supply values
1.90 + typedef V Value;
1.91 + /// The type of the arc costs
1.92 + typedef C Cost;
1.93 +
1.94 + public:
1.95 +
1.96 + /// \brief Problem type constants for the \c run() function.
1.97 + ///
1.98 + /// Enum type containing the problem type constants that can be
1.99 + /// returned by the \ref run() function of the algorithm.
1.100 + enum ProblemType {
1.101 + /// The problem has no feasible solution (flow).
1.102 + INFEASIBLE,
1.103 + /// The problem has optimal solution (i.e. it is feasible and
1.104 + /// bounded), and the algorithm has found optimal flow and node
1.105 + /// potentials (primal and dual solutions).
1.106 + OPTIMAL,
1.107 + /// The digraph contains an arc of negative cost and infinite
1.108 + /// upper bound. It means that the objective function is unbounded
1.109 + /// on that arc, however, note that it could actually be bounded
1.110 + /// over the feasible flows, but this algroithm cannot handle
1.111 + /// these cases.
1.112 + UNBOUNDED
1.113 + };
1.114 +
1.115 + /// \brief Constants for selecting the used method.
1.116 + ///
1.117 + /// Enum type containing constants for selecting the used method
1.118 + /// for the \ref run() function.
1.119 + ///
1.120 + /// \ref CycleCanceling provides three different cycle-canceling
1.121 + /// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel and Tighten"
1.122 + /// is used, which proved to be the most efficient and the most robust
1.123 + /// on various test inputs.
1.124 + /// However, the other methods can be selected using the \ref run()
1.125 + /// function with the proper parameter.
1.126 + enum Method {
1.127 + /// A simple cycle-canceling method, which uses the
1.128 + /// \ref BellmanFord "Bellman-Ford" algorithm with limited iteration
1.129 + /// number for detecting negative cycles in the residual network.
1.130 + SIMPLE_CYCLE_CANCELING,
1.131 + /// The "Minimum Mean Cycle-Canceling" algorithm, which is a
1.132 + /// well-known strongly polynomial method
1.133 + /// \ref goldberg89cyclecanceling. It improves along a
1.134 + /// \ref min_mean_cycle "minimum mean cycle" in each iteration.
1.135 + /// Its running time complexity is O(n<sup>2</sup>m<sup>3</sup>log(n)).
1.136 + MINIMUM_MEAN_CYCLE_CANCELING,
1.137 + /// The "Cancel And Tighten" algorithm, which can be viewed as an
1.138 + /// improved version of the previous method
1.139 + /// \ref goldberg89cyclecanceling.
1.140 + /// It is faster both in theory and in practice, its running time
1.141 + /// complexity is O(n<sup>2</sup>m<sup>2</sup>log(n)).
1.142 + CANCEL_AND_TIGHTEN
1.143 + };
1.144 +
1.145 + private:
1.146 +
1.147 + TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.148 +
1.149 + typedef std::vector<int> IntVector;
1.150 + typedef std::vector<double> DoubleVector;
1.151 + typedef std::vector<Value> ValueVector;
1.152 + typedef std::vector<Cost> CostVector;
1.153 + typedef std::vector<char> BoolVector;
1.154 + // Note: vector<char> is used instead of vector<bool> for efficiency reasons
1.155 +
1.156 + private:
1.157 +
1.158 + template <typename KT, typename VT>
1.159 + class StaticVectorMap {
1.160 + public:
1.161 + typedef KT Key;
1.162 + typedef VT Value;
1.163 +
1.164 + StaticVectorMap(std::vector<Value>& v) : _v(v) {}
1.165 +
1.166 + const Value& operator[](const Key& key) const {
1.167 + return _v[StaticDigraph::id(key)];
1.168 + }
1.169 +
1.170 + Value& operator[](const Key& key) {
1.171 + return _v[StaticDigraph::id(key)];
1.172 + }
1.173 +
1.174 + void set(const Key& key, const Value& val) {
1.175 + _v[StaticDigraph::id(key)] = val;
1.176 + }
1.177 +
1.178 + private:
1.179 + std::vector<Value>& _v;
1.180 + };
1.181 +
1.182 + typedef StaticVectorMap<StaticDigraph::Node, Cost> CostNodeMap;
1.183 + typedef StaticVectorMap<StaticDigraph::Arc, Cost> CostArcMap;
1.184 +
1.185 + private:
1.186 +
1.187 +
1.188 + // Data related to the underlying digraph
1.189 + const GR &_graph;
1.190 + int _node_num;
1.191 + int _arc_num;
1.192 + int _res_node_num;
1.193 + int _res_arc_num;
1.194 + int _root;
1.195 +
1.196 + // Parameters of the problem
1.197 + bool _have_lower;
1.198 + Value _sum_supply;
1.199 +
1.200 + // Data structures for storing the digraph
1.201 + IntNodeMap _node_id;
1.202 + IntArcMap _arc_idf;
1.203 + IntArcMap _arc_idb;
1.204 + IntVector _first_out;
1.205 + BoolVector _forward;
1.206 + IntVector _source;
1.207 + IntVector _target;
1.208 + IntVector _reverse;
1.209 +
1.210 + // Node and arc data
1.211 + ValueVector _lower;
1.212 + ValueVector _upper;
1.213 + CostVector _cost;
1.214 + ValueVector _supply;
1.215 +
1.216 + ValueVector _res_cap;
1.217 + CostVector _pi;
1.218 +
1.219 + // Data for a StaticDigraph structure
1.220 + typedef std::pair<int, int> IntPair;
1.221 + StaticDigraph _sgr;
1.222 + std::vector<IntPair> _arc_vec;
1.223 + std::vector<Cost> _cost_vec;
1.224 + IntVector _id_vec;
1.225 + CostArcMap _cost_map;
1.226 + CostNodeMap _pi_map;
1.227 +
1.228 + public:
1.229 +
1.230 + /// \brief Constant for infinite upper bounds (capacities).
1.231 + ///
1.232 + /// Constant for infinite upper bounds (capacities).
1.233 + /// It is \c std::numeric_limits<Value>::infinity() if available,
1.234 + /// \c std::numeric_limits<Value>::max() otherwise.
1.235 + const Value INF;
1.236 +
1.237 + public:
1.238 +
1.239 + /// \brief Constructor.
1.240 + ///
1.241 + /// The constructor of the class.
1.242 + ///
1.243 + /// \param graph The digraph the algorithm runs on.
1.244 + CycleCanceling(const GR& graph) :
1.245 + _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
1.246 + _cost_map(_cost_vec), _pi_map(_pi),
1.247 + INF(std::numeric_limits<Value>::has_infinity ?
1.248 + std::numeric_limits<Value>::infinity() :
1.249 + std::numeric_limits<Value>::max())
1.250 + {
1.251 + // Check the number types
1.252 + LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
1.253 + "The flow type of CycleCanceling must be signed");
1.254 + LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
1.255 + "The cost type of CycleCanceling must be signed");
1.256 +
1.257 + // Reset data structures
1.258 + reset();
1.259 + }
1.260 +
1.261 + /// \name Parameters
1.262 + /// The parameters of the algorithm can be specified using these
1.263 + /// functions.
1.264 +
1.265 + /// @{
1.266 +
1.267 + /// \brief Set the lower bounds on the arcs.
1.268 + ///
1.269 + /// This function sets the lower bounds on the arcs.
1.270 + /// If it is not used before calling \ref run(), the lower bounds
1.271 + /// will be set to zero on all arcs.
1.272 + ///
1.273 + /// \param map An arc map storing the lower bounds.
1.274 + /// Its \c Value type must be convertible to the \c Value type
1.275 + /// of the algorithm.
1.276 + ///
1.277 + /// \return <tt>(*this)</tt>
1.278 + template <typename LowerMap>
1.279 + CycleCanceling& lowerMap(const LowerMap& map) {
1.280 + _have_lower = true;
1.281 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.282 + _lower[_arc_idf[a]] = map[a];
1.283 + _lower[_arc_idb[a]] = map[a];
1.284 + }
1.285 + return *this;
1.286 + }
1.287 +
1.288 + /// \brief Set the upper bounds (capacities) on the arcs.
1.289 + ///
1.290 + /// This function sets the upper bounds (capacities) on the arcs.
1.291 + /// If it is not used before calling \ref run(), the upper bounds
1.292 + /// will be set to \ref INF on all arcs (i.e. the flow value will be
1.293 + /// unbounded from above).
1.294 + ///
1.295 + /// \param map An arc map storing the upper bounds.
1.296 + /// Its \c Value type must be convertible to the \c Value type
1.297 + /// of the algorithm.
1.298 + ///
1.299 + /// \return <tt>(*this)</tt>
1.300 + template<typename UpperMap>
1.301 + CycleCanceling& upperMap(const UpperMap& map) {
1.302 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.303 + _upper[_arc_idf[a]] = map[a];
1.304 + }
1.305 + return *this;
1.306 + }
1.307 +
1.308 + /// \brief Set the costs of the arcs.
1.309 + ///
1.310 + /// This function sets the costs of the arcs.
1.311 + /// If it is not used before calling \ref run(), the costs
1.312 + /// will be set to \c 1 on all arcs.
1.313 + ///
1.314 + /// \param map An arc map storing the costs.
1.315 + /// Its \c Value type must be convertible to the \c Cost type
1.316 + /// of the algorithm.
1.317 + ///
1.318 + /// \return <tt>(*this)</tt>
1.319 + template<typename CostMap>
1.320 + CycleCanceling& costMap(const CostMap& map) {
1.321 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.322 + _cost[_arc_idf[a]] = map[a];
1.323 + _cost[_arc_idb[a]] = -map[a];
1.324 + }
1.325 + return *this;
1.326 + }
1.327 +
1.328 + /// \brief Set the supply values of the nodes.
1.329 + ///
1.330 + /// This function sets the supply values of the nodes.
1.331 + /// If neither this function nor \ref stSupply() is used before
1.332 + /// calling \ref run(), the supply of each node will be set to zero.
1.333 + ///
1.334 + /// \param map A node map storing the supply values.
1.335 + /// Its \c Value type must be convertible to the \c Value type
1.336 + /// of the algorithm.
1.337 + ///
1.338 + /// \return <tt>(*this)</tt>
1.339 + template<typename SupplyMap>
1.340 + CycleCanceling& supplyMap(const SupplyMap& map) {
1.341 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.342 + _supply[_node_id[n]] = map[n];
1.343 + }
1.344 + return *this;
1.345 + }
1.346 +
1.347 + /// \brief Set single source and target nodes and a supply value.
1.348 + ///
1.349 + /// This function sets a single source node and a single target node
1.350 + /// and the required flow value.
1.351 + /// If neither this function nor \ref supplyMap() is used before
1.352 + /// calling \ref run(), the supply of each node will be set to zero.
1.353 + ///
1.354 + /// Using this function has the same effect as using \ref supplyMap()
1.355 + /// with such a map in which \c k is assigned to \c s, \c -k is
1.356 + /// assigned to \c t and all other nodes have zero supply value.
1.357 + ///
1.358 + /// \param s The source node.
1.359 + /// \param t The target node.
1.360 + /// \param k The required amount of flow from node \c s to node \c t
1.361 + /// (i.e. the supply of \c s and the demand of \c t).
1.362 + ///
1.363 + /// \return <tt>(*this)</tt>
1.364 + CycleCanceling& stSupply(const Node& s, const Node& t, Value k) {
1.365 + for (int i = 0; i != _res_node_num; ++i) {
1.366 + _supply[i] = 0;
1.367 + }
1.368 + _supply[_node_id[s]] = k;
1.369 + _supply[_node_id[t]] = -k;
1.370 + return *this;
1.371 + }
1.372 +
1.373 + /// @}
1.374 +
1.375 + /// \name Execution control
1.376 + /// The algorithm can be executed using \ref run().
1.377 +
1.378 + /// @{
1.379 +
1.380 + /// \brief Run the algorithm.
1.381 + ///
1.382 + /// This function runs the algorithm.
1.383 + /// The paramters can be specified using functions \ref lowerMap(),
1.384 + /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
1.385 + /// For example,
1.386 + /// \code
1.387 + /// CycleCanceling<ListDigraph> cc(graph);
1.388 + /// cc.lowerMap(lower).upperMap(upper).costMap(cost)
1.389 + /// .supplyMap(sup).run();
1.390 + /// \endcode
1.391 + ///
1.392 + /// This function can be called more than once. All the given parameters
1.393 + /// are kept for the next call, unless \ref resetParams() or \ref reset()
1.394 + /// is used, thus only the modified parameters have to be set again.
1.395 + /// If the underlying digraph was also modified after the construction
1.396 + /// of the class (or the last \ref reset() call), then the \ref reset()
1.397 + /// function must be called.
1.398 + ///
1.399 + /// \param method The cycle-canceling method that will be used.
1.400 + /// For more information, see \ref Method.
1.401 + ///
1.402 + /// \return \c INFEASIBLE if no feasible flow exists,
1.403 + /// \n \c OPTIMAL if the problem has optimal solution
1.404 + /// (i.e. it is feasible and bounded), and the algorithm has found
1.405 + /// optimal flow and node potentials (primal and dual solutions),
1.406 + /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
1.407 + /// and infinite upper bound. It means that the objective function
1.408 + /// is unbounded on that arc, however, note that it could actually be
1.409 + /// bounded over the feasible flows, but this algroithm cannot handle
1.410 + /// these cases.
1.411 + ///
1.412 + /// \see ProblemType, Method
1.413 + /// \see resetParams(), reset()
1.414 + ProblemType run(Method method = CANCEL_AND_TIGHTEN) {
1.415 + ProblemType pt = init();
1.416 + if (pt != OPTIMAL) return pt;
1.417 + start(method);
1.418 + return OPTIMAL;
1.419 + }
1.420 +
1.421 + /// \brief Reset all the parameters that have been given before.
1.422 + ///
1.423 + /// This function resets all the paramaters that have been given
1.424 + /// before using functions \ref lowerMap(), \ref upperMap(),
1.425 + /// \ref costMap(), \ref supplyMap(), \ref stSupply().
1.426 + ///
1.427 + /// It is useful for multiple \ref run() calls. Basically, all the given
1.428 + /// parameters are kept for the next \ref run() call, unless
1.429 + /// \ref resetParams() or \ref reset() is used.
1.430 + /// If the underlying digraph was also modified after the construction
1.431 + /// of the class or the last \ref reset() call, then the \ref reset()
1.432 + /// function must be used, otherwise \ref resetParams() is sufficient.
1.433 + ///
1.434 + /// For example,
1.435 + /// \code
1.436 + /// CycleCanceling<ListDigraph> cs(graph);
1.437 + ///
1.438 + /// // First run
1.439 + /// cc.lowerMap(lower).upperMap(upper).costMap(cost)
1.440 + /// .supplyMap(sup).run();
1.441 + ///
1.442 + /// // Run again with modified cost map (resetParams() is not called,
1.443 + /// // so only the cost map have to be set again)
1.444 + /// cost[e] += 100;
1.445 + /// cc.costMap(cost).run();
1.446 + ///
1.447 + /// // Run again from scratch using resetParams()
1.448 + /// // (the lower bounds will be set to zero on all arcs)
1.449 + /// cc.resetParams();
1.450 + /// cc.upperMap(capacity).costMap(cost)
1.451 + /// .supplyMap(sup).run();
1.452 + /// \endcode
1.453 + ///
1.454 + /// \return <tt>(*this)</tt>
1.455 + ///
1.456 + /// \see reset(), run()
1.457 + CycleCanceling& resetParams() {
1.458 + for (int i = 0; i != _res_node_num; ++i) {
1.459 + _supply[i] = 0;
1.460 + }
1.461 + int limit = _first_out[_root];
1.462 + for (int j = 0; j != limit; ++j) {
1.463 + _lower[j] = 0;
1.464 + _upper[j] = INF;
1.465 + _cost[j] = _forward[j] ? 1 : -1;
1.466 + }
1.467 + for (int j = limit; j != _res_arc_num; ++j) {
1.468 + _lower[j] = 0;
1.469 + _upper[j] = INF;
1.470 + _cost[j] = 0;
1.471 + _cost[_reverse[j]] = 0;
1.472 + }
1.473 + _have_lower = false;
1.474 + return *this;
1.475 + }
1.476 +
1.477 + /// \brief Reset the internal data structures and all the parameters
1.478 + /// that have been given before.
1.479 + ///
1.480 + /// This function resets the internal data structures and all the
1.481 + /// paramaters that have been given before using functions \ref lowerMap(),
1.482 + /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
1.483 + ///
1.484 + /// It is useful for multiple \ref run() calls. Basically, all the given
1.485 + /// parameters are kept for the next \ref run() call, unless
1.486 + /// \ref resetParams() or \ref reset() is used.
1.487 + /// If the underlying digraph was also modified after the construction
1.488 + /// of the class or the last \ref reset() call, then the \ref reset()
1.489 + /// function must be used, otherwise \ref resetParams() is sufficient.
1.490 + ///
1.491 + /// See \ref resetParams() for examples.
1.492 + ///
1.493 + /// \return <tt>(*this)</tt>
1.494 + ///
1.495 + /// \see resetParams(), run()
1.496 + CycleCanceling& reset() {
1.497 + // Resize vectors
1.498 + _node_num = countNodes(_graph);
1.499 + _arc_num = countArcs(_graph);
1.500 + _res_node_num = _node_num + 1;
1.501 + _res_arc_num = 2 * (_arc_num + _node_num);
1.502 + _root = _node_num;
1.503 +
1.504 + _first_out.resize(_res_node_num + 1);
1.505 + _forward.resize(_res_arc_num);
1.506 + _source.resize(_res_arc_num);
1.507 + _target.resize(_res_arc_num);
1.508 + _reverse.resize(_res_arc_num);
1.509 +
1.510 + _lower.resize(_res_arc_num);
1.511 + _upper.resize(_res_arc_num);
1.512 + _cost.resize(_res_arc_num);
1.513 + _supply.resize(_res_node_num);
1.514 +
1.515 + _res_cap.resize(_res_arc_num);
1.516 + _pi.resize(_res_node_num);
1.517 +
1.518 + _arc_vec.reserve(_res_arc_num);
1.519 + _cost_vec.reserve(_res_arc_num);
1.520 + _id_vec.reserve(_res_arc_num);
1.521 +
1.522 + // Copy the graph
1.523 + int i = 0, j = 0, k = 2 * _arc_num + _node_num;
1.524 + for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.525 + _node_id[n] = i;
1.526 + }
1.527 + i = 0;
1.528 + for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.529 + _first_out[i] = j;
1.530 + for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
1.531 + _arc_idf[a] = j;
1.532 + _forward[j] = true;
1.533 + _source[j] = i;
1.534 + _target[j] = _node_id[_graph.runningNode(a)];
1.535 + }
1.536 + for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
1.537 + _arc_idb[a] = j;
1.538 + _forward[j] = false;
1.539 + _source[j] = i;
1.540 + _target[j] = _node_id[_graph.runningNode(a)];
1.541 + }
1.542 + _forward[j] = false;
1.543 + _source[j] = i;
1.544 + _target[j] = _root;
1.545 + _reverse[j] = k;
1.546 + _forward[k] = true;
1.547 + _source[k] = _root;
1.548 + _target[k] = i;
1.549 + _reverse[k] = j;
1.550 + ++j; ++k;
1.551 + }
1.552 + _first_out[i] = j;
1.553 + _first_out[_res_node_num] = k;
1.554 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.555 + int fi = _arc_idf[a];
1.556 + int bi = _arc_idb[a];
1.557 + _reverse[fi] = bi;
1.558 + _reverse[bi] = fi;
1.559 + }
1.560 +
1.561 + // Reset parameters
1.562 + resetParams();
1.563 + return *this;
1.564 + }
1.565 +
1.566 + /// @}
1.567 +
1.568 + /// \name Query Functions
1.569 + /// The results of the algorithm can be obtained using these
1.570 + /// functions.\n
1.571 + /// The \ref run() function must be called before using them.
1.572 +
1.573 + /// @{
1.574 +
1.575 + /// \brief Return the total cost of the found flow.
1.576 + ///
1.577 + /// This function returns the total cost of the found flow.
1.578 + /// Its complexity is O(e).
1.579 + ///
1.580 + /// \note The return type of the function can be specified as a
1.581 + /// template parameter. For example,
1.582 + /// \code
1.583 + /// cc.totalCost<double>();
1.584 + /// \endcode
1.585 + /// It is useful if the total cost cannot be stored in the \c Cost
1.586 + /// type of the algorithm, which is the default return type of the
1.587 + /// function.
1.588 + ///
1.589 + /// \pre \ref run() must be called before using this function.
1.590 + template <typename Number>
1.591 + Number totalCost() const {
1.592 + Number c = 0;
1.593 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.594 + int i = _arc_idb[a];
1.595 + c += static_cast<Number>(_res_cap[i]) *
1.596 + (-static_cast<Number>(_cost[i]));
1.597 + }
1.598 + return c;
1.599 + }
1.600 +
1.601 +#ifndef DOXYGEN
1.602 + Cost totalCost() const {
1.603 + return totalCost<Cost>();
1.604 + }
1.605 +#endif
1.606 +
1.607 + /// \brief Return the flow on the given arc.
1.608 + ///
1.609 + /// This function returns the flow on the given arc.
1.610 + ///
1.611 + /// \pre \ref run() must be called before using this function.
1.612 + Value flow(const Arc& a) const {
1.613 + return _res_cap[_arc_idb[a]];
1.614 + }
1.615 +
1.616 + /// \brief Return the flow map (the primal solution).
1.617 + ///
1.618 + /// This function copies the flow value on each arc into the given
1.619 + /// map. The \c Value type of the algorithm must be convertible to
1.620 + /// the \c Value type of the map.
1.621 + ///
1.622 + /// \pre \ref run() must be called before using this function.
1.623 + template <typename FlowMap>
1.624 + void flowMap(FlowMap &map) const {
1.625 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.626 + map.set(a, _res_cap[_arc_idb[a]]);
1.627 + }
1.628 + }
1.629 +
1.630 + /// \brief Return the potential (dual value) of the given node.
1.631 + ///
1.632 + /// This function returns the potential (dual value) of the
1.633 + /// given node.
1.634 + ///
1.635 + /// \pre \ref run() must be called before using this function.
1.636 + Cost potential(const Node& n) const {
1.637 + return static_cast<Cost>(_pi[_node_id[n]]);
1.638 + }
1.639 +
1.640 + /// \brief Return the potential map (the dual solution).
1.641 + ///
1.642 + /// This function copies the potential (dual value) of each node
1.643 + /// into the given map.
1.644 + /// The \c Cost type of the algorithm must be convertible to the
1.645 + /// \c Value type of the map.
1.646 + ///
1.647 + /// \pre \ref run() must be called before using this function.
1.648 + template <typename PotentialMap>
1.649 + void potentialMap(PotentialMap &map) const {
1.650 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.651 + map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
1.652 + }
1.653 + }
1.654 +
1.655 + /// @}
1.656 +
1.657 + private:
1.658 +
1.659 + // Initialize the algorithm
1.660 + ProblemType init() {
1.661 + if (_res_node_num <= 1) return INFEASIBLE;
1.662 +
1.663 + // Check the sum of supply values
1.664 + _sum_supply = 0;
1.665 + for (int i = 0; i != _root; ++i) {
1.666 + _sum_supply += _supply[i];
1.667 + }
1.668 + if (_sum_supply > 0) return INFEASIBLE;
1.669 +
1.670 +
1.671 + // Initialize vectors
1.672 + for (int i = 0; i != _res_node_num; ++i) {
1.673 + _pi[i] = 0;
1.674 + }
1.675 + ValueVector excess(_supply);
1.676 +
1.677 + // Remove infinite upper bounds and check negative arcs
1.678 + const Value MAX = std::numeric_limits<Value>::max();
1.679 + int last_out;
1.680 + if (_have_lower) {
1.681 + for (int i = 0; i != _root; ++i) {
1.682 + last_out = _first_out[i+1];
1.683 + for (int j = _first_out[i]; j != last_out; ++j) {
1.684 + if (_forward[j]) {
1.685 + Value c = _cost[j] < 0 ? _upper[j] : _lower[j];
1.686 + if (c >= MAX) return UNBOUNDED;
1.687 + excess[i] -= c;
1.688 + excess[_target[j]] += c;
1.689 + }
1.690 + }
1.691 + }
1.692 + } else {
1.693 + for (int i = 0; i != _root; ++i) {
1.694 + last_out = _first_out[i+1];
1.695 + for (int j = _first_out[i]; j != last_out; ++j) {
1.696 + if (_forward[j] && _cost[j] < 0) {
1.697 + Value c = _upper[j];
1.698 + if (c >= MAX) return UNBOUNDED;
1.699 + excess[i] -= c;
1.700 + excess[_target[j]] += c;
1.701 + }
1.702 + }
1.703 + }
1.704 + }
1.705 + Value ex, max_cap = 0;
1.706 + for (int i = 0; i != _res_node_num; ++i) {
1.707 + ex = excess[i];
1.708 + if (ex < 0) max_cap -= ex;
1.709 + }
1.710 + for (int j = 0; j != _res_arc_num; ++j) {
1.711 + if (_upper[j] >= MAX) _upper[j] = max_cap;
1.712 + }
1.713 +
1.714 + // Initialize maps for Circulation and remove non-zero lower bounds
1.715 + ConstMap<Arc, Value> low(0);
1.716 + typedef typename Digraph::template ArcMap<Value> ValueArcMap;
1.717 + typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
1.718 + ValueArcMap cap(_graph), flow(_graph);
1.719 + ValueNodeMap sup(_graph);
1.720 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.721 + sup[n] = _supply[_node_id[n]];
1.722 + }
1.723 + if (_have_lower) {
1.724 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.725 + int j = _arc_idf[a];
1.726 + Value c = _lower[j];
1.727 + cap[a] = _upper[j] - c;
1.728 + sup[_graph.source(a)] -= c;
1.729 + sup[_graph.target(a)] += c;
1.730 + }
1.731 + } else {
1.732 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.733 + cap[a] = _upper[_arc_idf[a]];
1.734 + }
1.735 + }
1.736 +
1.737 + // Find a feasible flow using Circulation
1.738 + Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
1.739 + circ(_graph, low, cap, sup);
1.740 + if (!circ.flowMap(flow).run()) return INFEASIBLE;
1.741 +
1.742 + // Set residual capacities and handle GEQ supply type
1.743 + if (_sum_supply < 0) {
1.744 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.745 + Value fa = flow[a];
1.746 + _res_cap[_arc_idf[a]] = cap[a] - fa;
1.747 + _res_cap[_arc_idb[a]] = fa;
1.748 + sup[_graph.source(a)] -= fa;
1.749 + sup[_graph.target(a)] += fa;
1.750 + }
1.751 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.752 + excess[_node_id[n]] = sup[n];
1.753 + }
1.754 + for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
1.755 + int u = _target[a];
1.756 + int ra = _reverse[a];
1.757 + _res_cap[a] = -_sum_supply + 1;
1.758 + _res_cap[ra] = -excess[u];
1.759 + _cost[a] = 0;
1.760 + _cost[ra] = 0;
1.761 + }
1.762 + } else {
1.763 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.764 + Value fa = flow[a];
1.765 + _res_cap[_arc_idf[a]] = cap[a] - fa;
1.766 + _res_cap[_arc_idb[a]] = fa;
1.767 + }
1.768 + for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
1.769 + int ra = _reverse[a];
1.770 + _res_cap[a] = 1;
1.771 + _res_cap[ra] = 0;
1.772 + _cost[a] = 0;
1.773 + _cost[ra] = 0;
1.774 + }
1.775 + }
1.776 +
1.777 + return OPTIMAL;
1.778 + }
1.779 +
1.780 + // Build a StaticDigraph structure containing the current
1.781 + // residual network
1.782 + void buildResidualNetwork() {
1.783 + _arc_vec.clear();
1.784 + _cost_vec.clear();
1.785 + _id_vec.clear();
1.786 + for (int j = 0; j != _res_arc_num; ++j) {
1.787 + if (_res_cap[j] > 0) {
1.788 + _arc_vec.push_back(IntPair(_source[j], _target[j]));
1.789 + _cost_vec.push_back(_cost[j]);
1.790 + _id_vec.push_back(j);
1.791 + }
1.792 + }
1.793 + _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1.794 + }
1.795 +
1.796 + // Execute the algorithm and transform the results
1.797 + void start(Method method) {
1.798 + // Execute the algorithm
1.799 + switch (method) {
1.800 + case SIMPLE_CYCLE_CANCELING:
1.801 + startSimpleCycleCanceling();
1.802 + break;
1.803 + case MINIMUM_MEAN_CYCLE_CANCELING:
1.804 + startMinMeanCycleCanceling();
1.805 + break;
1.806 + case CANCEL_AND_TIGHTEN:
1.807 + startCancelAndTighten();
1.808 + break;
1.809 + }
1.810 +
1.811 + // Compute node potentials
1.812 + if (method != SIMPLE_CYCLE_CANCELING) {
1.813 + buildResidualNetwork();
1.814 + typename BellmanFord<StaticDigraph, CostArcMap>
1.815 + ::template SetDistMap<CostNodeMap>::Create bf(_sgr, _cost_map);
1.816 + bf.distMap(_pi_map);
1.817 + bf.init(0);
1.818 + bf.start();
1.819 + }
1.820 +
1.821 + // Handle non-zero lower bounds
1.822 + if (_have_lower) {
1.823 + int limit = _first_out[_root];
1.824 + for (int j = 0; j != limit; ++j) {
1.825 + if (!_forward[j]) _res_cap[j] += _lower[j];
1.826 + }
1.827 + }
1.828 + }
1.829 +
1.830 + // Execute the "Simple Cycle Canceling" method
1.831 + void startSimpleCycleCanceling() {
1.832 + // Constants for computing the iteration limits
1.833 + const int BF_FIRST_LIMIT = 2;
1.834 + const double BF_LIMIT_FACTOR = 1.5;
1.835 +
1.836 + typedef StaticVectorMap<StaticDigraph::Arc, Value> FilterMap;
1.837 + typedef FilterArcs<StaticDigraph, FilterMap> ResDigraph;
1.838 + typedef StaticVectorMap<StaticDigraph::Node, StaticDigraph::Arc> PredMap;
1.839 + typedef typename BellmanFord<ResDigraph, CostArcMap>
1.840 + ::template SetDistMap<CostNodeMap>
1.841 + ::template SetPredMap<PredMap>::Create BF;
1.842 +
1.843 + // Build the residual network
1.844 + _arc_vec.clear();
1.845 + _cost_vec.clear();
1.846 + for (int j = 0; j != _res_arc_num; ++j) {
1.847 + _arc_vec.push_back(IntPair(_source[j], _target[j]));
1.848 + _cost_vec.push_back(_cost[j]);
1.849 + }
1.850 + _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1.851 +
1.852 + FilterMap filter_map(_res_cap);
1.853 + ResDigraph rgr(_sgr, filter_map);
1.854 + std::vector<int> cycle;
1.855 + std::vector<StaticDigraph::Arc> pred(_res_arc_num);
1.856 + PredMap pred_map(pred);
1.857 + BF bf(rgr, _cost_map);
1.858 + bf.distMap(_pi_map).predMap(pred_map);
1.859 +
1.860 + int length_bound = BF_FIRST_LIMIT;
1.861 + bool optimal = false;
1.862 + while (!optimal) {
1.863 + bf.init(0);
1.864 + int iter_num = 0;
1.865 + bool cycle_found = false;
1.866 + while (!cycle_found) {
1.867 + // Perform some iterations of the Bellman-Ford algorithm
1.868 + int curr_iter_num = iter_num + length_bound <= _node_num ?
1.869 + length_bound : _node_num - iter_num;
1.870 + iter_num += curr_iter_num;
1.871 + int real_iter_num = curr_iter_num;
1.872 + for (int i = 0; i < curr_iter_num; ++i) {
1.873 + if (bf.processNextWeakRound()) {
1.874 + real_iter_num = i;
1.875 + break;
1.876 + }
1.877 + }
1.878 + if (real_iter_num < curr_iter_num) {
1.879 + // Optimal flow is found
1.880 + optimal = true;
1.881 + break;
1.882 + } else {
1.883 + // Search for node disjoint negative cycles
1.884 + std::vector<int> state(_res_node_num, 0);
1.885 + int id = 0;
1.886 + for (int u = 0; u != _res_node_num; ++u) {
1.887 + if (state[u] != 0) continue;
1.888 + ++id;
1.889 + int v = u;
1.890 + for (; v != -1 && state[v] == 0; v = pred[v] == INVALID ?
1.891 + -1 : rgr.id(rgr.source(pred[v]))) {
1.892 + state[v] = id;
1.893 + }
1.894 + if (v != -1 && state[v] == id) {
1.895 + // A negative cycle is found
1.896 + cycle_found = true;
1.897 + cycle.clear();
1.898 + StaticDigraph::Arc a = pred[v];
1.899 + Value d, delta = _res_cap[rgr.id(a)];
1.900 + cycle.push_back(rgr.id(a));
1.901 + while (rgr.id(rgr.source(a)) != v) {
1.902 + a = pred_map[rgr.source(a)];
1.903 + d = _res_cap[rgr.id(a)];
1.904 + if (d < delta) delta = d;
1.905 + cycle.push_back(rgr.id(a));
1.906 + }
1.907 +
1.908 + // Augment along the cycle
1.909 + for (int i = 0; i < int(cycle.size()); ++i) {
1.910 + int j = cycle[i];
1.911 + _res_cap[j] -= delta;
1.912 + _res_cap[_reverse[j]] += delta;
1.913 + }
1.914 + }
1.915 + }
1.916 + }
1.917 +
1.918 + // Increase iteration limit if no cycle is found
1.919 + if (!cycle_found) {
1.920 + length_bound = static_cast<int>(length_bound * BF_LIMIT_FACTOR);
1.921 + }
1.922 + }
1.923 + }
1.924 + }
1.925 +
1.926 + // Execute the "Minimum Mean Cycle Canceling" method
1.927 + void startMinMeanCycleCanceling() {
1.928 + typedef SimplePath<StaticDigraph> SPath;
1.929 + typedef typename SPath::ArcIt SPathArcIt;
1.930 + typedef typename HowardMmc<StaticDigraph, CostArcMap>
1.931 + ::template SetPath<SPath>::Create MMC;
1.932 +
1.933 + SPath cycle;
1.934 + MMC mmc(_sgr, _cost_map);
1.935 + mmc.cycle(cycle);
1.936 + buildResidualNetwork();
1.937 + while (mmc.findCycleMean() && mmc.cycleCost() < 0) {
1.938 + // Find the cycle
1.939 + mmc.findCycle();
1.940 +
1.941 + // Compute delta value
1.942 + Value delta = INF;
1.943 + for (SPathArcIt a(cycle); a != INVALID; ++a) {
1.944 + Value d = _res_cap[_id_vec[_sgr.id(a)]];
1.945 + if (d < delta) delta = d;
1.946 + }
1.947 +
1.948 + // Augment along the cycle
1.949 + for (SPathArcIt a(cycle); a != INVALID; ++a) {
1.950 + int j = _id_vec[_sgr.id(a)];
1.951 + _res_cap[j] -= delta;
1.952 + _res_cap[_reverse[j]] += delta;
1.953 + }
1.954 +
1.955 + // Rebuild the residual network
1.956 + buildResidualNetwork();
1.957 + }
1.958 + }
1.959 +
1.960 + // Execute the "Cancel And Tighten" method
1.961 + void startCancelAndTighten() {
1.962 + // Constants for the min mean cycle computations
1.963 + const double LIMIT_FACTOR = 1.0;
1.964 + const int MIN_LIMIT = 5;
1.965 +
1.966 + // Contruct auxiliary data vectors
1.967 + DoubleVector pi(_res_node_num, 0.0);
1.968 + IntVector level(_res_node_num);
1.969 + BoolVector reached(_res_node_num);
1.970 + BoolVector processed(_res_node_num);
1.971 + IntVector pred_node(_res_node_num);
1.972 + IntVector pred_arc(_res_node_num);
1.973 + std::vector<int> stack(_res_node_num);
1.974 + std::vector<int> proc_vector(_res_node_num);
1.975 +
1.976 + // Initialize epsilon
1.977 + double epsilon = 0;
1.978 + for (int a = 0; a != _res_arc_num; ++a) {
1.979 + if (_res_cap[a] > 0 && -_cost[a] > epsilon)
1.980 + epsilon = -_cost[a];
1.981 + }
1.982 +
1.983 + // Start phases
1.984 + Tolerance<double> tol;
1.985 + tol.epsilon(1e-6);
1.986 + int limit = int(LIMIT_FACTOR * std::sqrt(double(_res_node_num)));
1.987 + if (limit < MIN_LIMIT) limit = MIN_LIMIT;
1.988 + int iter = limit;
1.989 + while (epsilon * _res_node_num >= 1) {
1.990 + // Find and cancel cycles in the admissible network using DFS
1.991 + for (int u = 0; u != _res_node_num; ++u) {
1.992 + reached[u] = false;
1.993 + processed[u] = false;
1.994 + }
1.995 + int stack_head = -1;
1.996 + int proc_head = -1;
1.997 + for (int start = 0; start != _res_node_num; ++start) {
1.998 + if (reached[start]) continue;
1.999 +
1.1000 + // New start node
1.1001 + reached[start] = true;
1.1002 + pred_arc[start] = -1;
1.1003 + pred_node[start] = -1;
1.1004 +
1.1005 + // Find the first admissible outgoing arc
1.1006 + double p = pi[start];
1.1007 + int a = _first_out[start];
1.1008 + int last_out = _first_out[start+1];
1.1009 + for (; a != last_out && (_res_cap[a] == 0 ||
1.1010 + !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
1.1011 + if (a == last_out) {
1.1012 + processed[start] = true;
1.1013 + proc_vector[++proc_head] = start;
1.1014 + continue;
1.1015 + }
1.1016 + stack[++stack_head] = a;
1.1017 +
1.1018 + while (stack_head >= 0) {
1.1019 + int sa = stack[stack_head];
1.1020 + int u = _source[sa];
1.1021 + int v = _target[sa];
1.1022 +
1.1023 + if (!reached[v]) {
1.1024 + // A new node is reached
1.1025 + reached[v] = true;
1.1026 + pred_node[v] = u;
1.1027 + pred_arc[v] = sa;
1.1028 + p = pi[v];
1.1029 + a = _first_out[v];
1.1030 + last_out = _first_out[v+1];
1.1031 + for (; a != last_out && (_res_cap[a] == 0 ||
1.1032 + !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
1.1033 + stack[++stack_head] = a == last_out ? -1 : a;
1.1034 + } else {
1.1035 + if (!processed[v]) {
1.1036 + // A cycle is found
1.1037 + int n, w = u;
1.1038 + Value d, delta = _res_cap[sa];
1.1039 + for (n = u; n != v; n = pred_node[n]) {
1.1040 + d = _res_cap[pred_arc[n]];
1.1041 + if (d <= delta) {
1.1042 + delta = d;
1.1043 + w = pred_node[n];
1.1044 + }
1.1045 + }
1.1046 +
1.1047 + // Augment along the cycle
1.1048 + _res_cap[sa] -= delta;
1.1049 + _res_cap[_reverse[sa]] += delta;
1.1050 + for (n = u; n != v; n = pred_node[n]) {
1.1051 + int pa = pred_arc[n];
1.1052 + _res_cap[pa] -= delta;
1.1053 + _res_cap[_reverse[pa]] += delta;
1.1054 + }
1.1055 + for (n = u; stack_head > 0 && n != w; n = pred_node[n]) {
1.1056 + --stack_head;
1.1057 + reached[n] = false;
1.1058 + }
1.1059 + u = w;
1.1060 + }
1.1061 + v = u;
1.1062 +
1.1063 + // Find the next admissible outgoing arc
1.1064 + p = pi[v];
1.1065 + a = stack[stack_head] + 1;
1.1066 + last_out = _first_out[v+1];
1.1067 + for (; a != last_out && (_res_cap[a] == 0 ||
1.1068 + !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
1.1069 + stack[stack_head] = a == last_out ? -1 : a;
1.1070 + }
1.1071 +
1.1072 + while (stack_head >= 0 && stack[stack_head] == -1) {
1.1073 + processed[v] = true;
1.1074 + proc_vector[++proc_head] = v;
1.1075 + if (--stack_head >= 0) {
1.1076 + // Find the next admissible outgoing arc
1.1077 + v = _source[stack[stack_head]];
1.1078 + p = pi[v];
1.1079 + a = stack[stack_head] + 1;
1.1080 + last_out = _first_out[v+1];
1.1081 + for (; a != last_out && (_res_cap[a] == 0 ||
1.1082 + !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
1.1083 + stack[stack_head] = a == last_out ? -1 : a;
1.1084 + }
1.1085 + }
1.1086 + }
1.1087 + }
1.1088 +
1.1089 + // Tighten potentials and epsilon
1.1090 + if (--iter > 0) {
1.1091 + for (int u = 0; u != _res_node_num; ++u) {
1.1092 + level[u] = 0;
1.1093 + }
1.1094 + for (int i = proc_head; i > 0; --i) {
1.1095 + int u = proc_vector[i];
1.1096 + double p = pi[u];
1.1097 + int l = level[u] + 1;
1.1098 + int last_out = _first_out[u+1];
1.1099 + for (int a = _first_out[u]; a != last_out; ++a) {
1.1100 + int v = _target[a];
1.1101 + if (_res_cap[a] > 0 && tol.negative(_cost[a] + p - pi[v]) &&
1.1102 + l > level[v]) level[v] = l;
1.1103 + }
1.1104 + }
1.1105 +
1.1106 + // Modify potentials
1.1107 + double q = std::numeric_limits<double>::max();
1.1108 + for (int u = 0; u != _res_node_num; ++u) {
1.1109 + int lu = level[u];
1.1110 + double p, pu = pi[u];
1.1111 + int last_out = _first_out[u+1];
1.1112 + for (int a = _first_out[u]; a != last_out; ++a) {
1.1113 + if (_res_cap[a] == 0) continue;
1.1114 + int v = _target[a];
1.1115 + int ld = lu - level[v];
1.1116 + if (ld > 0) {
1.1117 + p = (_cost[a] + pu - pi[v] + epsilon) / (ld + 1);
1.1118 + if (p < q) q = p;
1.1119 + }
1.1120 + }
1.1121 + }
1.1122 + for (int u = 0; u != _res_node_num; ++u) {
1.1123 + pi[u] -= q * level[u];
1.1124 + }
1.1125 +
1.1126 + // Modify epsilon
1.1127 + epsilon = 0;
1.1128 + for (int u = 0; u != _res_node_num; ++u) {
1.1129 + double curr, pu = pi[u];
1.1130 + int last_out = _first_out[u+1];
1.1131 + for (int a = _first_out[u]; a != last_out; ++a) {
1.1132 + if (_res_cap[a] == 0) continue;
1.1133 + curr = _cost[a] + pu - pi[_target[a]];
1.1134 + if (-curr > epsilon) epsilon = -curr;
1.1135 + }
1.1136 + }
1.1137 + } else {
1.1138 + typedef HowardMmc<StaticDigraph, CostArcMap> MMC;
1.1139 + typedef typename BellmanFord<StaticDigraph, CostArcMap>
1.1140 + ::template SetDistMap<CostNodeMap>::Create BF;
1.1141 +
1.1142 + // Set epsilon to the minimum cycle mean
1.1143 + buildResidualNetwork();
1.1144 + MMC mmc(_sgr, _cost_map);
1.1145 + mmc.findCycleMean();
1.1146 + epsilon = -mmc.cycleMean();
1.1147 + Cost cycle_cost = mmc.cycleCost();
1.1148 + int cycle_size = mmc.cycleSize();
1.1149 +
1.1150 + // Compute feasible potentials for the current epsilon
1.1151 + for (int i = 0; i != int(_cost_vec.size()); ++i) {
1.1152 + _cost_vec[i] = cycle_size * _cost_vec[i] - cycle_cost;
1.1153 + }
1.1154 + BF bf(_sgr, _cost_map);
1.1155 + bf.distMap(_pi_map);
1.1156 + bf.init(0);
1.1157 + bf.start();
1.1158 + for (int u = 0; u != _res_node_num; ++u) {
1.1159 + pi[u] = static_cast<double>(_pi[u]) / cycle_size;
1.1160 + }
1.1161 +
1.1162 + iter = limit;
1.1163 + }
1.1164 + }
1.1165 + }
1.1166 +
1.1167 + }; //class CycleCanceling
1.1168 +
1.1169 + ///@}
1.1170 +
1.1171 +} //namespace lemon
1.1172 +
1.1173 +#endif //LEMON_CYCLE_CANCELING_H