1.1 --- a/lemon/cycle_canceling.h Fri Nov 13 00:09:35 2009 +0100
1.2 +++ b/lemon/cycle_canceling.h Fri Nov 13 00:10:33 2009 +0100
1.3 @@ -19,441 +19,817 @@
1.4 #ifndef LEMON_CYCLE_CANCELING_H
1.5 #define LEMON_CYCLE_CANCELING_H
1.6
1.7 -/// \ingroup min_cost_flow
1.8 -///
1.9 +/// \ingroup min_cost_flow_algs
1.10 /// \file
1.11 -/// \brief Cycle-canceling algorithm for finding a minimum cost flow.
1.12 +/// \brief Cycle-canceling algorithms for finding a minimum cost flow.
1.13
1.14 #include <vector>
1.15 +#include <limits>
1.16 +
1.17 +#include <lemon/core.h>
1.18 +#include <lemon/maps.h>
1.19 +#include <lemon/path.h>
1.20 +#include <lemon/math.h>
1.21 +#include <lemon/static_graph.h>
1.22 #include <lemon/adaptors.h>
1.23 -#include <lemon/path.h>
1.24 -
1.25 #include <lemon/circulation.h>
1.26 #include <lemon/bellman_ford.h>
1.27 #include <lemon/howard.h>
1.28
1.29 namespace lemon {
1.30
1.31 - /// \addtogroup min_cost_flow
1.32 + /// \addtogroup min_cost_flow_algs
1.33 /// @{
1.34
1.35 - /// \brief Implementation of a cycle-canceling algorithm for
1.36 - /// finding a minimum cost flow.
1.37 + /// \brief Implementation of cycle-canceling algorithms for
1.38 + /// finding a \ref min_cost_flow "minimum cost flow".
1.39 ///
1.40 - /// \ref CycleCanceling implements a cycle-canceling algorithm for
1.41 - /// finding a minimum cost flow.
1.42 + /// \ref CycleCanceling implements three different cycle-canceling
1.43 + /// algorithms for finding a \ref min_cost_flow "minimum cost flow".
1.44 + /// The most efficent one (both theoretically and practically)
1.45 + /// is the \ref CANCEL_AND_TIGHTEN "Cancel and Tighten" algorithm,
1.46 + /// thus it is the default method.
1.47 + /// It is strongly polynomial, but in practice, it is typically much
1.48 + /// slower than the scaling algorithms and NetworkSimplex.
1.49 ///
1.50 - /// \tparam Digraph The digraph type the algorithm runs on.
1.51 - /// \tparam LowerMap The type of the lower bound map.
1.52 - /// \tparam CapacityMap The type of the capacity (upper bound) map.
1.53 - /// \tparam CostMap The type of the cost (length) map.
1.54 - /// \tparam SupplyMap The type of the supply map.
1.55 + /// Most of the parameters of the problem (except for the digraph)
1.56 + /// can be given using separate functions, and the algorithm can be
1.57 + /// executed using the \ref run() function. If some parameters are not
1.58 + /// specified, then default values will be used.
1.59 ///
1.60 - /// \warning
1.61 - /// - Arc capacities and costs should be \e non-negative \e integers.
1.62 - /// - Supply values should be \e signed \e integers.
1.63 - /// - The value types of the maps should be convertible to each other.
1.64 - /// - \c CostMap::Value must be signed type.
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 - /// \note By default the \ref BellmanFord "Bellman-Ford" algorithm is
1.72 - /// used for negative cycle detection with limited iteration number.
1.73 - /// However \ref CycleCanceling also provides the "Minimum Mean
1.74 - /// Cycle-Canceling" algorithm, which is \e strongly \e polynomial,
1.75 - /// but rather slower in practice.
1.76 - /// To use this version of the algorithm, call \ref run() with \c true
1.77 - /// parameter.
1.78 + /// \warning Both number types must be signed and all input data must
1.79 + /// be integer.
1.80 + /// \warning This algorithm does not support negative costs for such
1.81 + /// arcs that have infinite upper bound.
1.82 ///
1.83 - /// \author Peter Kovacs
1.84 - template < typename Digraph,
1.85 - typename LowerMap = typename Digraph::template ArcMap<int>,
1.86 - typename CapacityMap = typename Digraph::template ArcMap<int>,
1.87 - typename CostMap = typename Digraph::template ArcMap<int>,
1.88 - typename SupplyMap = typename Digraph::template NodeMap<int> >
1.89 + /// \note For more information about the three available methods,
1.90 + /// see \ref Method.
1.91 +#ifdef DOXYGEN
1.92 + template <typename GR, typename V, typename C>
1.93 +#else
1.94 + template <typename GR, typename V = int, typename C = V>
1.95 +#endif
1.96 class CycleCanceling
1.97 {
1.98 - TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
1.99 + public:
1.100
1.101 - typedef typename CapacityMap::Value Capacity;
1.102 - typedef typename CostMap::Value Cost;
1.103 - typedef typename SupplyMap::Value Supply;
1.104 - typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
1.105 - typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
1.106 -
1.107 - typedef ResidualDigraph< const Digraph,
1.108 - CapacityArcMap, CapacityArcMap > ResDigraph;
1.109 - typedef typename ResDigraph::Node ResNode;
1.110 - typedef typename ResDigraph::NodeIt ResNodeIt;
1.111 - typedef typename ResDigraph::Arc ResArc;
1.112 - typedef typename ResDigraph::ArcIt ResArcIt;
1.113 + /// The type of the digraph
1.114 + typedef GR Digraph;
1.115 + /// The type of the flow amounts, capacity bounds and supply values
1.116 + typedef V Value;
1.117 + /// The type of the arc costs
1.118 + typedef C Cost;
1.119
1.120 public:
1.121
1.122 - /// The type of the flow map.
1.123 - typedef typename Digraph::template ArcMap<Capacity> FlowMap;
1.124 - /// The type of the potential map.
1.125 - typedef typename Digraph::template NodeMap<Cost> PotentialMap;
1.126 + /// \brief Problem type constants for the \c run() function.
1.127 + ///
1.128 + /// Enum type containing the problem type constants that can be
1.129 + /// returned by the \ref run() function of the algorithm.
1.130 + enum ProblemType {
1.131 + /// The problem has no feasible solution (flow).
1.132 + INFEASIBLE,
1.133 + /// The problem has optimal solution (i.e. it is feasible and
1.134 + /// bounded), and the algorithm has found optimal flow and node
1.135 + /// potentials (primal and dual solutions).
1.136 + OPTIMAL,
1.137 + /// The digraph contains an arc of negative cost and infinite
1.138 + /// upper bound. It means that the objective function is unbounded
1.139 + /// on that arc, however, note that it could actually be bounded
1.140 + /// over the feasible flows, but this algroithm cannot handle
1.141 + /// these cases.
1.142 + UNBOUNDED
1.143 + };
1.144 +
1.145 + /// \brief Constants for selecting the used method.
1.146 + ///
1.147 + /// Enum type containing constants for selecting the used method
1.148 + /// for the \ref run() function.
1.149 + ///
1.150 + /// \ref CycleCanceling provides three different cycle-canceling
1.151 + /// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel and Tighten"
1.152 + /// is used, which proved to be the most efficient and the most robust
1.153 + /// on various test inputs.
1.154 + /// However, the other methods can be selected using the \ref run()
1.155 + /// function with the proper parameter.
1.156 + enum Method {
1.157 + /// A simple cycle-canceling method, which uses the
1.158 + /// \ref BellmanFord "Bellman-Ford" algorithm with limited iteration
1.159 + /// number for detecting negative cycles in the residual network.
1.160 + SIMPLE_CYCLE_CANCELING,
1.161 + /// The "Minimum Mean Cycle-Canceling" algorithm, which is a
1.162 + /// well-known strongly polynomial method. It improves along a
1.163 + /// \ref min_mean_cycle "minimum mean cycle" in each iteration.
1.164 + /// Its running time complexity is O(n<sup>2</sup>m<sup>3</sup>log(n)).
1.165 + MINIMUM_MEAN_CYCLE_CANCELING,
1.166 + /// The "Cancel And Tighten" algorithm, which can be viewed as an
1.167 + /// improved version of the previous method.
1.168 + /// It is faster both in theory and in practice, its running time
1.169 + /// complexity is O(n<sup>2</sup>m<sup>2</sup>log(n)).
1.170 + CANCEL_AND_TIGHTEN
1.171 + };
1.172
1.173 private:
1.174
1.175 - /// \brief Map adaptor class for handling residual arc costs.
1.176 - ///
1.177 - /// Map adaptor class for handling residual arc costs.
1.178 - class ResidualCostMap : public MapBase<ResArc, Cost>
1.179 - {
1.180 - private:
1.181 + TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.182 +
1.183 + typedef std::vector<int> IntVector;
1.184 + typedef std::vector<char> CharVector;
1.185 + typedef std::vector<double> DoubleVector;
1.186 + typedef std::vector<Value> ValueVector;
1.187 + typedef std::vector<Cost> CostVector;
1.188
1.189 - const CostMap &_cost_map;
1.190 -
1.191 + private:
1.192 +
1.193 + template <typename KT, typename VT>
1.194 + class VectorMap {
1.195 public:
1.196 -
1.197 - ///\e
1.198 - ResidualCostMap(const CostMap &cost_map) : _cost_map(cost_map) {}
1.199 -
1.200 - ///\e
1.201 - Cost operator[](const ResArc &e) const {
1.202 - return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
1.203 + typedef KT Key;
1.204 + typedef VT Value;
1.205 +
1.206 + VectorMap(std::vector<Value>& v) : _v(v) {}
1.207 +
1.208 + const Value& operator[](const Key& key) const {
1.209 + return _v[StaticDigraph::id(key)];
1.210 }
1.211
1.212 - }; //class ResidualCostMap
1.213 + Value& operator[](const Key& key) {
1.214 + return _v[StaticDigraph::id(key)];
1.215 + }
1.216 +
1.217 + void set(const Key& key, const Value& val) {
1.218 + _v[StaticDigraph::id(key)] = val;
1.219 + }
1.220 +
1.221 + private:
1.222 + std::vector<Value>& _v;
1.223 + };
1.224 +
1.225 + typedef VectorMap<StaticDigraph::Node, Cost> CostNodeMap;
1.226 + typedef VectorMap<StaticDigraph::Arc, Cost> CostArcMap;
1.227
1.228 private:
1.229
1.230 - // The maximum number of iterations for the first execution of the
1.231 - // Bellman-Ford algorithm. It should be at least 2.
1.232 - static const int BF_FIRST_LIMIT = 2;
1.233 - // The iteration limit for the Bellman-Ford algorithm is multiplied
1.234 - // by BF_LIMIT_FACTOR/100 in every round.
1.235 - static const int BF_LIMIT_FACTOR = 150;
1.236
1.237 - private:
1.238 + // Data related to the underlying digraph
1.239 + const GR &_graph;
1.240 + int _node_num;
1.241 + int _arc_num;
1.242 + int _res_node_num;
1.243 + int _res_arc_num;
1.244 + int _root;
1.245
1.246 - // The digraph the algorithm runs on
1.247 - const Digraph &_graph;
1.248 - // The original lower bound map
1.249 - const LowerMap *_lower;
1.250 - // The modified capacity map
1.251 - CapacityArcMap _capacity;
1.252 - // The original cost map
1.253 - const CostMap &_cost;
1.254 - // The modified supply map
1.255 - SupplyNodeMap _supply;
1.256 - bool _valid_supply;
1.257 + // Parameters of the problem
1.258 + bool _have_lower;
1.259 + Value _sum_supply;
1.260
1.261 - // Arc map of the current flow
1.262 - FlowMap *_flow;
1.263 - bool _local_flow;
1.264 - // Node map of the current potentials
1.265 - PotentialMap *_potential;
1.266 - bool _local_potential;
1.267 + // Data structures for storing the digraph
1.268 + IntNodeMap _node_id;
1.269 + IntArcMap _arc_idf;
1.270 + IntArcMap _arc_idb;
1.271 + IntVector _first_out;
1.272 + CharVector _forward;
1.273 + IntVector _source;
1.274 + IntVector _target;
1.275 + IntVector _reverse;
1.276
1.277 - // The residual digraph
1.278 - ResDigraph *_res_graph;
1.279 - // The residual cost map
1.280 - ResidualCostMap _res_cost;
1.281 + // Node and arc data
1.282 + ValueVector _lower;
1.283 + ValueVector _upper;
1.284 + CostVector _cost;
1.285 + ValueVector _supply;
1.286 +
1.287 + ValueVector _res_cap;
1.288 + CostVector _pi;
1.289 +
1.290 + // Data for a StaticDigraph structure
1.291 + typedef std::pair<int, int> IntPair;
1.292 + StaticDigraph _sgr;
1.293 + std::vector<IntPair> _arc_vec;
1.294 + std::vector<Cost> _cost_vec;
1.295 + IntVector _id_vec;
1.296 + CostArcMap _cost_map;
1.297 + CostNodeMap _pi_map;
1.298 +
1.299 + public:
1.300 +
1.301 + /// \brief Constant for infinite upper bounds (capacities).
1.302 + ///
1.303 + /// Constant for infinite upper bounds (capacities).
1.304 + /// It is \c std::numeric_limits<Value>::infinity() if available,
1.305 + /// \c std::numeric_limits<Value>::max() otherwise.
1.306 + const Value INF;
1.307
1.308 public:
1.309
1.310 - /// \brief General constructor (with lower bounds).
1.311 + /// \brief Constructor.
1.312 ///
1.313 - /// General constructor (with lower bounds).
1.314 + /// The constructor of the class.
1.315 ///
1.316 - /// \param digraph The digraph the algorithm runs on.
1.317 - /// \param lower The lower bounds of the arcs.
1.318 - /// \param capacity The capacities (upper bounds) of the arcs.
1.319 - /// \param cost The cost (length) values of the arcs.
1.320 - /// \param supply The supply values of the nodes (signed).
1.321 - CycleCanceling( const Digraph &digraph,
1.322 - const LowerMap &lower,
1.323 - const CapacityMap &capacity,
1.324 - const CostMap &cost,
1.325 - const SupplyMap &supply ) :
1.326 - _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost),
1.327 - _supply(digraph), _flow(NULL), _local_flow(false),
1.328 - _potential(NULL), _local_potential(false),
1.329 - _res_graph(NULL), _res_cost(_cost)
1.330 + /// \param graph The digraph the algorithm runs on.
1.331 + CycleCanceling(const GR& graph) :
1.332 + _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
1.333 + _cost_map(_cost_vec), _pi_map(_pi),
1.334 + INF(std::numeric_limits<Value>::has_infinity ?
1.335 + std::numeric_limits<Value>::infinity() :
1.336 + std::numeric_limits<Value>::max())
1.337 {
1.338 - // Check the sum of supply values
1.339 - Supply sum = 0;
1.340 - for (NodeIt n(_graph); n != INVALID; ++n) {
1.341 - _supply[n] = supply[n];
1.342 - sum += _supply[n];
1.343 + // Check the number types
1.344 + LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
1.345 + "The flow type of CycleCanceling must be signed");
1.346 + LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
1.347 + "The cost type of CycleCanceling must be signed");
1.348 +
1.349 + // Resize vectors
1.350 + _node_num = countNodes(_graph);
1.351 + _arc_num = countArcs(_graph);
1.352 + _res_node_num = _node_num + 1;
1.353 + _res_arc_num = 2 * (_arc_num + _node_num);
1.354 + _root = _node_num;
1.355 +
1.356 + _first_out.resize(_res_node_num + 1);
1.357 + _forward.resize(_res_arc_num);
1.358 + _source.resize(_res_arc_num);
1.359 + _target.resize(_res_arc_num);
1.360 + _reverse.resize(_res_arc_num);
1.361 +
1.362 + _lower.resize(_res_arc_num);
1.363 + _upper.resize(_res_arc_num);
1.364 + _cost.resize(_res_arc_num);
1.365 + _supply.resize(_res_node_num);
1.366 +
1.367 + _res_cap.resize(_res_arc_num);
1.368 + _pi.resize(_res_node_num);
1.369 +
1.370 + _arc_vec.reserve(_res_arc_num);
1.371 + _cost_vec.reserve(_res_arc_num);
1.372 + _id_vec.reserve(_res_arc_num);
1.373 +
1.374 + // Copy the graph
1.375 + int i = 0, j = 0, k = 2 * _arc_num + _node_num;
1.376 + for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.377 + _node_id[n] = i;
1.378 }
1.379 - _valid_supply = sum == 0;
1.380 -
1.381 - // Remove non-zero lower bounds
1.382 - for (ArcIt e(_graph); e != INVALID; ++e) {
1.383 - _capacity[e] = capacity[e];
1.384 - if (lower[e] != 0) {
1.385 - _capacity[e] -= lower[e];
1.386 - _supply[_graph.source(e)] -= lower[e];
1.387 - _supply[_graph.target(e)] += lower[e];
1.388 + i = 0;
1.389 + for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.390 + _first_out[i] = j;
1.391 + for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
1.392 + _arc_idf[a] = j;
1.393 + _forward[j] = true;
1.394 + _source[j] = i;
1.395 + _target[j] = _node_id[_graph.runningNode(a)];
1.396 }
1.397 + for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
1.398 + _arc_idb[a] = j;
1.399 + _forward[j] = false;
1.400 + _source[j] = i;
1.401 + _target[j] = _node_id[_graph.runningNode(a)];
1.402 + }
1.403 + _forward[j] = false;
1.404 + _source[j] = i;
1.405 + _target[j] = _root;
1.406 + _reverse[j] = k;
1.407 + _forward[k] = true;
1.408 + _source[k] = _root;
1.409 + _target[k] = i;
1.410 + _reverse[k] = j;
1.411 + ++j; ++k;
1.412 }
1.413 - }
1.414 -/*
1.415 - /// \brief General constructor (without lower bounds).
1.416 - ///
1.417 - /// General constructor (without lower bounds).
1.418 - ///
1.419 - /// \param digraph The digraph the algorithm runs on.
1.420 - /// \param capacity The capacities (upper bounds) of the arcs.
1.421 - /// \param cost The cost (length) values of the arcs.
1.422 - /// \param supply The supply values of the nodes (signed).
1.423 - CycleCanceling( const Digraph &digraph,
1.424 - const CapacityMap &capacity,
1.425 - const CostMap &cost,
1.426 - const SupplyMap &supply ) :
1.427 - _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
1.428 - _supply(supply), _flow(NULL), _local_flow(false),
1.429 - _potential(NULL), _local_potential(false), _res_graph(NULL),
1.430 - _res_cost(_cost)
1.431 - {
1.432 - // Check the sum of supply values
1.433 - Supply sum = 0;
1.434 - for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
1.435 - _valid_supply = sum == 0;
1.436 + _first_out[i] = j;
1.437 + _first_out[_res_node_num] = k;
1.438 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.439 + int fi = _arc_idf[a];
1.440 + int bi = _arc_idb[a];
1.441 + _reverse[fi] = bi;
1.442 + _reverse[bi] = fi;
1.443 + }
1.444 +
1.445 + // Reset parameters
1.446 + reset();
1.447 }
1.448
1.449 - /// \brief Simple constructor (with lower bounds).
1.450 + /// \name Parameters
1.451 + /// The parameters of the algorithm can be specified using these
1.452 + /// functions.
1.453 +
1.454 + /// @{
1.455 +
1.456 + /// \brief Set the lower bounds on the arcs.
1.457 ///
1.458 - /// Simple constructor (with lower bounds).
1.459 + /// This function sets the lower bounds on the arcs.
1.460 + /// If it is not used before calling \ref run(), the lower bounds
1.461 + /// will be set to zero on all arcs.
1.462 ///
1.463 - /// \param digraph The digraph the algorithm runs on.
1.464 - /// \param lower The lower bounds of the arcs.
1.465 - /// \param capacity The capacities (upper bounds) of the arcs.
1.466 - /// \param cost The cost (length) values of the arcs.
1.467 - /// \param s The source node.
1.468 - /// \param t The target node.
1.469 - /// \param flow_value The required amount of flow from node \c s
1.470 - /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.471 - CycleCanceling( const Digraph &digraph,
1.472 - const LowerMap &lower,
1.473 - const CapacityMap &capacity,
1.474 - const CostMap &cost,
1.475 - Node s, Node t,
1.476 - Supply flow_value ) :
1.477 - _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost),
1.478 - _supply(digraph, 0), _flow(NULL), _local_flow(false),
1.479 - _potential(NULL), _local_potential(false), _res_graph(NULL),
1.480 - _res_cost(_cost)
1.481 - {
1.482 - // Remove non-zero lower bounds
1.483 - _supply[s] = flow_value;
1.484 - _supply[t] = -flow_value;
1.485 - for (ArcIt e(_graph); e != INVALID; ++e) {
1.486 - if (lower[e] != 0) {
1.487 - _capacity[e] -= lower[e];
1.488 - _supply[_graph.source(e)] -= lower[e];
1.489 - _supply[_graph.target(e)] += lower[e];
1.490 - }
1.491 + /// \param map An arc map storing the lower bounds.
1.492 + /// Its \c Value type must be convertible to the \c Value type
1.493 + /// of the algorithm.
1.494 + ///
1.495 + /// \return <tt>(*this)</tt>
1.496 + template <typename LowerMap>
1.497 + CycleCanceling& lowerMap(const LowerMap& map) {
1.498 + _have_lower = true;
1.499 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.500 + _lower[_arc_idf[a]] = map[a];
1.501 + _lower[_arc_idb[a]] = map[a];
1.502 }
1.503 - _valid_supply = true;
1.504 - }
1.505 -
1.506 - /// \brief Simple constructor (without lower bounds).
1.507 - ///
1.508 - /// Simple constructor (without lower bounds).
1.509 - ///
1.510 - /// \param digraph The digraph the algorithm runs on.
1.511 - /// \param capacity The capacities (upper bounds) of the arcs.
1.512 - /// \param cost The cost (length) values of the arcs.
1.513 - /// \param s The source node.
1.514 - /// \param t The target node.
1.515 - /// \param flow_value The required amount of flow from node \c s
1.516 - /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.517 - CycleCanceling( const Digraph &digraph,
1.518 - const CapacityMap &capacity,
1.519 - const CostMap &cost,
1.520 - Node s, Node t,
1.521 - Supply flow_value ) :
1.522 - _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
1.523 - _supply(digraph, 0), _flow(NULL), _local_flow(false),
1.524 - _potential(NULL), _local_potential(false), _res_graph(NULL),
1.525 - _res_cost(_cost)
1.526 - {
1.527 - _supply[s] = flow_value;
1.528 - _supply[t] = -flow_value;
1.529 - _valid_supply = true;
1.530 - }
1.531 -*/
1.532 - /// Destructor.
1.533 - ~CycleCanceling() {
1.534 - if (_local_flow) delete _flow;
1.535 - if (_local_potential) delete _potential;
1.536 - delete _res_graph;
1.537 - }
1.538 -
1.539 - /// \brief Set the flow map.
1.540 - ///
1.541 - /// Set the flow map.
1.542 - ///
1.543 - /// \return \c (*this)
1.544 - CycleCanceling& flowMap(FlowMap &map) {
1.545 - if (_local_flow) {
1.546 - delete _flow;
1.547 - _local_flow = false;
1.548 - }
1.549 - _flow = ↦
1.550 return *this;
1.551 }
1.552
1.553 - /// \brief Set the potential map.
1.554 + /// \brief Set the upper bounds (capacities) on the arcs.
1.555 ///
1.556 - /// Set the potential map.
1.557 + /// This function sets the upper bounds (capacities) on the arcs.
1.558 + /// If it is not used before calling \ref run(), the upper bounds
1.559 + /// will be set to \ref INF on all arcs (i.e. the flow value will be
1.560 + /// unbounded from above).
1.561 ///
1.562 - /// \return \c (*this)
1.563 - CycleCanceling& potentialMap(PotentialMap &map) {
1.564 - if (_local_potential) {
1.565 - delete _potential;
1.566 - _local_potential = false;
1.567 + /// \param map An arc map storing the upper bounds.
1.568 + /// Its \c Value type must be convertible to the \c Value type
1.569 + /// of the algorithm.
1.570 + ///
1.571 + /// \return <tt>(*this)</tt>
1.572 + template<typename UpperMap>
1.573 + CycleCanceling& upperMap(const UpperMap& map) {
1.574 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.575 + _upper[_arc_idf[a]] = map[a];
1.576 }
1.577 - _potential = ↦
1.578 return *this;
1.579 }
1.580
1.581 + /// \brief Set the costs of the arcs.
1.582 + ///
1.583 + /// This function sets the costs of the arcs.
1.584 + /// If it is not used before calling \ref run(), the costs
1.585 + /// will be set to \c 1 on all arcs.
1.586 + ///
1.587 + /// \param map An arc map storing the costs.
1.588 + /// Its \c Value type must be convertible to the \c Cost type
1.589 + /// of the algorithm.
1.590 + ///
1.591 + /// \return <tt>(*this)</tt>
1.592 + template<typename CostMap>
1.593 + CycleCanceling& costMap(const CostMap& map) {
1.594 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.595 + _cost[_arc_idf[a]] = map[a];
1.596 + _cost[_arc_idb[a]] = -map[a];
1.597 + }
1.598 + return *this;
1.599 + }
1.600 +
1.601 + /// \brief Set the supply values of the nodes.
1.602 + ///
1.603 + /// This function sets the supply values of the nodes.
1.604 + /// If neither this function nor \ref stSupply() is used before
1.605 + /// calling \ref run(), the supply of each node will be set to zero.
1.606 + ///
1.607 + /// \param map A node map storing the supply values.
1.608 + /// Its \c Value type must be convertible to the \c Value type
1.609 + /// of the algorithm.
1.610 + ///
1.611 + /// \return <tt>(*this)</tt>
1.612 + template<typename SupplyMap>
1.613 + CycleCanceling& supplyMap(const SupplyMap& map) {
1.614 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.615 + _supply[_node_id[n]] = map[n];
1.616 + }
1.617 + return *this;
1.618 + }
1.619 +
1.620 + /// \brief Set single source and target nodes and a supply value.
1.621 + ///
1.622 + /// This function sets a single source node and a single target node
1.623 + /// and the required flow value.
1.624 + /// If neither this function nor \ref supplyMap() is used before
1.625 + /// calling \ref run(), the supply of each node will be set to zero.
1.626 + ///
1.627 + /// Using this function has the same effect as using \ref supplyMap()
1.628 + /// with such a map in which \c k is assigned to \c s, \c -k is
1.629 + /// assigned to \c t and all other nodes have zero supply value.
1.630 + ///
1.631 + /// \param s The source node.
1.632 + /// \param t The target node.
1.633 + /// \param k The required amount of flow from node \c s to node \c t
1.634 + /// (i.e. the supply of \c s and the demand of \c t).
1.635 + ///
1.636 + /// \return <tt>(*this)</tt>
1.637 + CycleCanceling& stSupply(const Node& s, const Node& t, Value k) {
1.638 + for (int i = 0; i != _res_node_num; ++i) {
1.639 + _supply[i] = 0;
1.640 + }
1.641 + _supply[_node_id[s]] = k;
1.642 + _supply[_node_id[t]] = -k;
1.643 + return *this;
1.644 + }
1.645 +
1.646 + /// @}
1.647 +
1.648 /// \name Execution control
1.649 + /// The algorithm can be executed using \ref run().
1.650
1.651 /// @{
1.652
1.653 /// \brief Run the algorithm.
1.654 ///
1.655 - /// Run the algorithm.
1.656 + /// This function runs the algorithm.
1.657 + /// The paramters can be specified using functions \ref lowerMap(),
1.658 + /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
1.659 + /// For example,
1.660 + /// \code
1.661 + /// CycleCanceling<ListDigraph> cc(graph);
1.662 + /// cc.lowerMap(lower).upperMap(upper).costMap(cost)
1.663 + /// .supplyMap(sup).run();
1.664 + /// \endcode
1.665 ///
1.666 - /// \param min_mean_cc Set this parameter to \c true to run the
1.667 - /// "Minimum Mean Cycle-Canceling" algorithm, which is strongly
1.668 - /// polynomial, but rather slower in practice.
1.669 + /// This function can be called more than once. All the parameters
1.670 + /// that have been given are kept for the next call, unless
1.671 + /// \ref reset() is called, thus only the modified parameters
1.672 + /// have to be set again. See \ref reset() for examples.
1.673 + /// However, the underlying digraph must not be modified after this
1.674 + /// class have been constructed, since it copies and extends the graph.
1.675 ///
1.676 - /// \return \c true if a feasible flow can be found.
1.677 - bool run(bool min_mean_cc = false) {
1.678 - return init() && start(min_mean_cc);
1.679 + /// \param method The cycle-canceling method that will be used.
1.680 + /// For more information, see \ref Method.
1.681 + ///
1.682 + /// \return \c INFEASIBLE if no feasible flow exists,
1.683 + /// \n \c OPTIMAL if the problem has optimal solution
1.684 + /// (i.e. it is feasible and bounded), and the algorithm has found
1.685 + /// optimal flow and node potentials (primal and dual solutions),
1.686 + /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
1.687 + /// and infinite upper bound. It means that the objective function
1.688 + /// is unbounded on that arc, however, note that it could actually be
1.689 + /// bounded over the feasible flows, but this algroithm cannot handle
1.690 + /// these cases.
1.691 + ///
1.692 + /// \see ProblemType, Method
1.693 + ProblemType run(Method method = CANCEL_AND_TIGHTEN) {
1.694 + ProblemType pt = init();
1.695 + if (pt != OPTIMAL) return pt;
1.696 + start(method);
1.697 + return OPTIMAL;
1.698 + }
1.699 +
1.700 + /// \brief Reset all the parameters that have been given before.
1.701 + ///
1.702 + /// This function resets all the paramaters that have been given
1.703 + /// before using functions \ref lowerMap(), \ref upperMap(),
1.704 + /// \ref costMap(), \ref supplyMap(), \ref stSupply().
1.705 + ///
1.706 + /// It is useful for multiple run() calls. If this function is not
1.707 + /// used, all the parameters given before are kept for the next
1.708 + /// \ref run() call.
1.709 + /// However, the underlying digraph must not be modified after this
1.710 + /// class have been constructed, since it copies and extends the graph.
1.711 + ///
1.712 + /// For example,
1.713 + /// \code
1.714 + /// CycleCanceling<ListDigraph> cs(graph);
1.715 + ///
1.716 + /// // First run
1.717 + /// cc.lowerMap(lower).upperMap(upper).costMap(cost)
1.718 + /// .supplyMap(sup).run();
1.719 + ///
1.720 + /// // Run again with modified cost map (reset() is not called,
1.721 + /// // so only the cost map have to be set again)
1.722 + /// cost[e] += 100;
1.723 + /// cc.costMap(cost).run();
1.724 + ///
1.725 + /// // Run again from scratch using reset()
1.726 + /// // (the lower bounds will be set to zero on all arcs)
1.727 + /// cc.reset();
1.728 + /// cc.upperMap(capacity).costMap(cost)
1.729 + /// .supplyMap(sup).run();
1.730 + /// \endcode
1.731 + ///
1.732 + /// \return <tt>(*this)</tt>
1.733 + CycleCanceling& reset() {
1.734 + for (int i = 0; i != _res_node_num; ++i) {
1.735 + _supply[i] = 0;
1.736 + }
1.737 + int limit = _first_out[_root];
1.738 + for (int j = 0; j != limit; ++j) {
1.739 + _lower[j] = 0;
1.740 + _upper[j] = INF;
1.741 + _cost[j] = _forward[j] ? 1 : -1;
1.742 + }
1.743 + for (int j = limit; j != _res_arc_num; ++j) {
1.744 + _lower[j] = 0;
1.745 + _upper[j] = INF;
1.746 + _cost[j] = 0;
1.747 + _cost[_reverse[j]] = 0;
1.748 + }
1.749 + _have_lower = false;
1.750 + return *this;
1.751 }
1.752
1.753 /// @}
1.754
1.755 /// \name Query Functions
1.756 - /// The result of the algorithm can be obtained using these
1.757 + /// The results of the algorithm can be obtained using these
1.758 /// functions.\n
1.759 - /// \ref lemon::CycleCanceling::run() "run()" must be called before
1.760 - /// using them.
1.761 + /// The \ref run() function must be called before using them.
1.762
1.763 /// @{
1.764
1.765 - /// \brief Return a const reference to the arc map storing the
1.766 - /// found flow.
1.767 + /// \brief Return the total cost of the found flow.
1.768 ///
1.769 - /// Return a const reference to the arc map storing the found flow.
1.770 + /// This function returns the total cost of the found flow.
1.771 + /// Its complexity is O(e).
1.772 + ///
1.773 + /// \note The return type of the function can be specified as a
1.774 + /// template parameter. For example,
1.775 + /// \code
1.776 + /// cc.totalCost<double>();
1.777 + /// \endcode
1.778 + /// It is useful if the total cost cannot be stored in the \c Cost
1.779 + /// type of the algorithm, which is the default return type of the
1.780 + /// function.
1.781 ///
1.782 /// \pre \ref run() must be called before using this function.
1.783 - const FlowMap& flowMap() const {
1.784 - return *_flow;
1.785 + template <typename Number>
1.786 + Number totalCost() const {
1.787 + Number c = 0;
1.788 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.789 + int i = _arc_idb[a];
1.790 + c += static_cast<Number>(_res_cap[i]) *
1.791 + (-static_cast<Number>(_cost[i]));
1.792 + }
1.793 + return c;
1.794 }
1.795
1.796 - /// \brief Return a const reference to the node map storing the
1.797 - /// found potentials (the dual solution).
1.798 - ///
1.799 - /// Return a const reference to the node map storing the found
1.800 - /// potentials (the dual solution).
1.801 - ///
1.802 - /// \pre \ref run() must be called before using this function.
1.803 - const PotentialMap& potentialMap() const {
1.804 - return *_potential;
1.805 +#ifndef DOXYGEN
1.806 + Cost totalCost() const {
1.807 + return totalCost<Cost>();
1.808 }
1.809 +#endif
1.810
1.811 /// \brief Return the flow on the given arc.
1.812 ///
1.813 - /// Return the flow on the given arc.
1.814 + /// This function returns the flow on the given arc.
1.815 ///
1.816 /// \pre \ref run() must be called before using this function.
1.817 - Capacity flow(const Arc& arc) const {
1.818 - return (*_flow)[arc];
1.819 + Value flow(const Arc& a) const {
1.820 + return _res_cap[_arc_idb[a]];
1.821 }
1.822
1.823 - /// \brief Return the potential of the given node.
1.824 + /// \brief Return the flow map (the primal solution).
1.825 ///
1.826 - /// Return the potential of the given node.
1.827 + /// This function copies the flow value on each arc into the given
1.828 + /// map. The \c Value type of the algorithm must be convertible to
1.829 + /// the \c Value type of the map.
1.830 ///
1.831 /// \pre \ref run() must be called before using this function.
1.832 - Cost potential(const Node& node) const {
1.833 - return (*_potential)[node];
1.834 + template <typename FlowMap>
1.835 + void flowMap(FlowMap &map) const {
1.836 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.837 + map.set(a, _res_cap[_arc_idb[a]]);
1.838 + }
1.839 }
1.840
1.841 - /// \brief Return the total cost of the found flow.
1.842 + /// \brief Return the potential (dual value) of the given node.
1.843 ///
1.844 - /// Return the total cost of the found flow. The complexity of the
1.845 - /// function is \f$ O(e) \f$.
1.846 + /// This function returns the potential (dual value) of the
1.847 + /// given node.
1.848 ///
1.849 /// \pre \ref run() must be called before using this function.
1.850 - Cost totalCost() const {
1.851 - Cost c = 0;
1.852 - for (ArcIt e(_graph); e != INVALID; ++e)
1.853 - c += (*_flow)[e] * _cost[e];
1.854 - return c;
1.855 + Cost potential(const Node& n) const {
1.856 + return static_cast<Cost>(_pi[_node_id[n]]);
1.857 + }
1.858 +
1.859 + /// \brief Return the potential map (the dual solution).
1.860 + ///
1.861 + /// This function copies the potential (dual value) of each node
1.862 + /// into the given map.
1.863 + /// The \c Cost type of the algorithm must be convertible to the
1.864 + /// \c Value type of the map.
1.865 + ///
1.866 + /// \pre \ref run() must be called before using this function.
1.867 + template <typename PotentialMap>
1.868 + void potentialMap(PotentialMap &map) const {
1.869 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.870 + map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
1.871 + }
1.872 }
1.873
1.874 /// @}
1.875
1.876 private:
1.877
1.878 - /// Initialize the algorithm.
1.879 - bool init() {
1.880 - if (!_valid_supply) return false;
1.881 + // Initialize the algorithm
1.882 + ProblemType init() {
1.883 + if (_res_node_num <= 1) return INFEASIBLE;
1.884
1.885 - // Initializing flow and potential maps
1.886 - if (!_flow) {
1.887 - _flow = new FlowMap(_graph);
1.888 - _local_flow = true;
1.889 + // Check the sum of supply values
1.890 + _sum_supply = 0;
1.891 + for (int i = 0; i != _root; ++i) {
1.892 + _sum_supply += _supply[i];
1.893 }
1.894 - if (!_potential) {
1.895 - _potential = new PotentialMap(_graph);
1.896 - _local_potential = true;
1.897 + if (_sum_supply > 0) return INFEASIBLE;
1.898 +
1.899 +
1.900 + // Initialize vectors
1.901 + for (int i = 0; i != _res_node_num; ++i) {
1.902 + _pi[i] = 0;
1.903 + }
1.904 + ValueVector excess(_supply);
1.905 +
1.906 + // Remove infinite upper bounds and check negative arcs
1.907 + const Value MAX = std::numeric_limits<Value>::max();
1.908 + int last_out;
1.909 + if (_have_lower) {
1.910 + for (int i = 0; i != _root; ++i) {
1.911 + last_out = _first_out[i+1];
1.912 + for (int j = _first_out[i]; j != last_out; ++j) {
1.913 + if (_forward[j]) {
1.914 + Value c = _cost[j] < 0 ? _upper[j] : _lower[j];
1.915 + if (c >= MAX) return UNBOUNDED;
1.916 + excess[i] -= c;
1.917 + excess[_target[j]] += c;
1.918 + }
1.919 + }
1.920 + }
1.921 + } else {
1.922 + for (int i = 0; i != _root; ++i) {
1.923 + last_out = _first_out[i+1];
1.924 + for (int j = _first_out[i]; j != last_out; ++j) {
1.925 + if (_forward[j] && _cost[j] < 0) {
1.926 + Value c = _upper[j];
1.927 + if (c >= MAX) return UNBOUNDED;
1.928 + excess[i] -= c;
1.929 + excess[_target[j]] += c;
1.930 + }
1.931 + }
1.932 + }
1.933 + }
1.934 + Value ex, max_cap = 0;
1.935 + for (int i = 0; i != _res_node_num; ++i) {
1.936 + ex = excess[i];
1.937 + if (ex < 0) max_cap -= ex;
1.938 + }
1.939 + for (int j = 0; j != _res_arc_num; ++j) {
1.940 + if (_upper[j] >= MAX) _upper[j] = max_cap;
1.941 }
1.942
1.943 - _res_graph = new ResDigraph(_graph, _capacity, *_flow);
1.944 + // Initialize maps for Circulation and remove non-zero lower bounds
1.945 + ConstMap<Arc, Value> low(0);
1.946 + typedef typename Digraph::template ArcMap<Value> ValueArcMap;
1.947 + typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
1.948 + ValueArcMap cap(_graph), flow(_graph);
1.949 + ValueNodeMap sup(_graph);
1.950 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.951 + sup[n] = _supply[_node_id[n]];
1.952 + }
1.953 + if (_have_lower) {
1.954 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.955 + int j = _arc_idf[a];
1.956 + Value c = _lower[j];
1.957 + cap[a] = _upper[j] - c;
1.958 + sup[_graph.source(a)] -= c;
1.959 + sup[_graph.target(a)] += c;
1.960 + }
1.961 + } else {
1.962 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.963 + cap[a] = _upper[_arc_idf[a]];
1.964 + }
1.965 + }
1.966
1.967 - // Finding a feasible flow using Circulation
1.968 - Circulation< Digraph, ConstMap<Arc, Capacity>, CapacityArcMap,
1.969 - SupplyMap >
1.970 - circulation( _graph, constMap<Arc>(Capacity(0)), _capacity,
1.971 - _supply );
1.972 - return circulation.flowMap(*_flow).run();
1.973 + // Find a feasible flow using Circulation
1.974 + Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
1.975 + circ(_graph, low, cap, sup);
1.976 + if (!circ.flowMap(flow).run()) return INFEASIBLE;
1.977 +
1.978 + // Set residual capacities and handle GEQ supply type
1.979 + if (_sum_supply < 0) {
1.980 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.981 + Value fa = flow[a];
1.982 + _res_cap[_arc_idf[a]] = cap[a] - fa;
1.983 + _res_cap[_arc_idb[a]] = fa;
1.984 + sup[_graph.source(a)] -= fa;
1.985 + sup[_graph.target(a)] += fa;
1.986 + }
1.987 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.988 + excess[_node_id[n]] = sup[n];
1.989 + }
1.990 + for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
1.991 + int u = _target[a];
1.992 + int ra = _reverse[a];
1.993 + _res_cap[a] = -_sum_supply + 1;
1.994 + _res_cap[ra] = -excess[u];
1.995 + _cost[a] = 0;
1.996 + _cost[ra] = 0;
1.997 + }
1.998 + } else {
1.999 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.1000 + Value fa = flow[a];
1.1001 + _res_cap[_arc_idf[a]] = cap[a] - fa;
1.1002 + _res_cap[_arc_idb[a]] = fa;
1.1003 + }
1.1004 + for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
1.1005 + int ra = _reverse[a];
1.1006 + _res_cap[a] = 1;
1.1007 + _res_cap[ra] = 0;
1.1008 + _cost[a] = 0;
1.1009 + _cost[ra] = 0;
1.1010 + }
1.1011 + }
1.1012 +
1.1013 + return OPTIMAL;
1.1014 + }
1.1015 +
1.1016 + // Build a StaticDigraph structure containing the current
1.1017 + // residual network
1.1018 + void buildResidualNetwork() {
1.1019 + _arc_vec.clear();
1.1020 + _cost_vec.clear();
1.1021 + _id_vec.clear();
1.1022 + for (int j = 0; j != _res_arc_num; ++j) {
1.1023 + if (_res_cap[j] > 0) {
1.1024 + _arc_vec.push_back(IntPair(_source[j], _target[j]));
1.1025 + _cost_vec.push_back(_cost[j]);
1.1026 + _id_vec.push_back(j);
1.1027 + }
1.1028 + }
1.1029 + _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1.1030 }
1.1031
1.1032 - bool start(bool min_mean_cc) {
1.1033 - if (min_mean_cc)
1.1034 - startMinMean();
1.1035 - else
1.1036 - start();
1.1037 + // Execute the algorithm and transform the results
1.1038 + void start(Method method) {
1.1039 + // Execute the algorithm
1.1040 + switch (method) {
1.1041 + case SIMPLE_CYCLE_CANCELING:
1.1042 + startSimpleCycleCanceling();
1.1043 + break;
1.1044 + case MINIMUM_MEAN_CYCLE_CANCELING:
1.1045 + startMinMeanCycleCanceling();
1.1046 + break;
1.1047 + case CANCEL_AND_TIGHTEN:
1.1048 + startCancelAndTighten();
1.1049 + break;
1.1050 + }
1.1051
1.1052 - // Handling non-zero lower bounds
1.1053 - if (_lower) {
1.1054 - for (ArcIt e(_graph); e != INVALID; ++e)
1.1055 - (*_flow)[e] += (*_lower)[e];
1.1056 + // Compute node potentials
1.1057 + if (method != SIMPLE_CYCLE_CANCELING) {
1.1058 + buildResidualNetwork();
1.1059 + typename BellmanFord<StaticDigraph, CostArcMap>
1.1060 + ::template SetDistMap<CostNodeMap>::Create bf(_sgr, _cost_map);
1.1061 + bf.distMap(_pi_map);
1.1062 + bf.init(0);
1.1063 + bf.start();
1.1064 }
1.1065 - return true;
1.1066 +
1.1067 + // Handle non-zero lower bounds
1.1068 + if (_have_lower) {
1.1069 + int limit = _first_out[_root];
1.1070 + for (int j = 0; j != limit; ++j) {
1.1071 + if (!_forward[j]) _res_cap[j] += _lower[j];
1.1072 + }
1.1073 + }
1.1074 }
1.1075
1.1076 - /// \brief Execute the algorithm using \ref BellmanFord.
1.1077 - ///
1.1078 - /// Execute the algorithm using the \ref BellmanFord
1.1079 - /// "Bellman-Ford" algorithm for negative cycle detection with
1.1080 - /// successively larger limit for the number of iterations.
1.1081 - void start() {
1.1082 - typename BellmanFord<ResDigraph, ResidualCostMap>::PredMap pred(*_res_graph);
1.1083 - typename ResDigraph::template NodeMap<int> visited(*_res_graph);
1.1084 - std::vector<ResArc> cycle;
1.1085 - int node_num = countNodes(_graph);
1.1086 + // Execute the "Simple Cycle Canceling" method
1.1087 + void startSimpleCycleCanceling() {
1.1088 + // Constants for computing the iteration limits
1.1089 + const int BF_FIRST_LIMIT = 2;
1.1090 + const double BF_LIMIT_FACTOR = 1.5;
1.1091 +
1.1092 + typedef VectorMap<StaticDigraph::Arc, Value> FilterMap;
1.1093 + typedef FilterArcs<StaticDigraph, FilterMap> ResDigraph;
1.1094 + typedef VectorMap<StaticDigraph::Node, StaticDigraph::Arc> PredMap;
1.1095 + typedef typename BellmanFord<ResDigraph, CostArcMap>
1.1096 + ::template SetDistMap<CostNodeMap>
1.1097 + ::template SetPredMap<PredMap>::Create BF;
1.1098 +
1.1099 + // Build the residual network
1.1100 + _arc_vec.clear();
1.1101 + _cost_vec.clear();
1.1102 + for (int j = 0; j != _res_arc_num; ++j) {
1.1103 + _arc_vec.push_back(IntPair(_source[j], _target[j]));
1.1104 + _cost_vec.push_back(_cost[j]);
1.1105 + }
1.1106 + _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1.1107 +
1.1108 + FilterMap filter_map(_res_cap);
1.1109 + ResDigraph rgr(_sgr, filter_map);
1.1110 + std::vector<int> cycle;
1.1111 + std::vector<StaticDigraph::Arc> pred(_res_arc_num);
1.1112 + PredMap pred_map(pred);
1.1113 + BF bf(rgr, _cost_map);
1.1114 + bf.distMap(_pi_map).predMap(pred_map);
1.1115
1.1116 int length_bound = BF_FIRST_LIMIT;
1.1117 bool optimal = false;
1.1118 while (!optimal) {
1.1119 - BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost);
1.1120 - bf.predMap(pred);
1.1121 bf.init(0);
1.1122 int iter_num = 0;
1.1123 bool cycle_found = false;
1.1124 while (!cycle_found) {
1.1125 - int curr_iter_num = iter_num + length_bound <= node_num ?
1.1126 - length_bound : node_num - iter_num;
1.1127 + // Perform some iterations of the Bellman-Ford algorithm
1.1128 + int curr_iter_num = iter_num + length_bound <= _node_num ?
1.1129 + length_bound : _node_num - iter_num;
1.1130 iter_num += curr_iter_num;
1.1131 int real_iter_num = curr_iter_num;
1.1132 for (int i = 0; i < curr_iter_num; ++i) {
1.1133 @@ -465,89 +841,290 @@
1.1134 if (real_iter_num < curr_iter_num) {
1.1135 // Optimal flow is found
1.1136 optimal = true;
1.1137 - // Setting node potentials
1.1138 - for (NodeIt n(_graph); n != INVALID; ++n)
1.1139 - (*_potential)[n] = bf.dist(n);
1.1140 break;
1.1141 } else {
1.1142 - // Searching for node disjoint negative cycles
1.1143 - for (ResNodeIt n(*_res_graph); n != INVALID; ++n)
1.1144 - visited[n] = 0;
1.1145 + // Search for node disjoint negative cycles
1.1146 + std::vector<int> state(_res_node_num, 0);
1.1147 int id = 0;
1.1148 - for (ResNodeIt n(*_res_graph); n != INVALID; ++n) {
1.1149 - if (visited[n] > 0) continue;
1.1150 - visited[n] = ++id;
1.1151 - ResNode u = pred[n] == INVALID ?
1.1152 - INVALID : _res_graph->source(pred[n]);
1.1153 - while (u != INVALID && visited[u] == 0) {
1.1154 - visited[u] = id;
1.1155 - u = pred[u] == INVALID ?
1.1156 - INVALID : _res_graph->source(pred[u]);
1.1157 + for (int u = 0; u != _res_node_num; ++u) {
1.1158 + if (state[u] != 0) continue;
1.1159 + ++id;
1.1160 + int v = u;
1.1161 + for (; v != -1 && state[v] == 0; v = pred[v] == INVALID ?
1.1162 + -1 : rgr.id(rgr.source(pred[v]))) {
1.1163 + state[v] = id;
1.1164 }
1.1165 - if (u != INVALID && visited[u] == id) {
1.1166 - // Finding the negative cycle
1.1167 + if (v != -1 && state[v] == id) {
1.1168 + // A negative cycle is found
1.1169 cycle_found = true;
1.1170 cycle.clear();
1.1171 - ResArc e = pred[u];
1.1172 - cycle.push_back(e);
1.1173 - Capacity d = _res_graph->residualCapacity(e);
1.1174 - while (_res_graph->source(e) != u) {
1.1175 - cycle.push_back(e = pred[_res_graph->source(e)]);
1.1176 - if (_res_graph->residualCapacity(e) < d)
1.1177 - d = _res_graph->residualCapacity(e);
1.1178 + StaticDigraph::Arc a = pred[v];
1.1179 + Value d, delta = _res_cap[rgr.id(a)];
1.1180 + cycle.push_back(rgr.id(a));
1.1181 + while (rgr.id(rgr.source(a)) != v) {
1.1182 + a = pred_map[rgr.source(a)];
1.1183 + d = _res_cap[rgr.id(a)];
1.1184 + if (d < delta) delta = d;
1.1185 + cycle.push_back(rgr.id(a));
1.1186 }
1.1187
1.1188 - // Augmenting along the cycle
1.1189 - for (int i = 0; i < int(cycle.size()); ++i)
1.1190 - _res_graph->augment(cycle[i], d);
1.1191 + // Augment along the cycle
1.1192 + for (int i = 0; i < int(cycle.size()); ++i) {
1.1193 + int j = cycle[i];
1.1194 + _res_cap[j] -= delta;
1.1195 + _res_cap[_reverse[j]] += delta;
1.1196 + }
1.1197 }
1.1198 }
1.1199 }
1.1200
1.1201 - if (!cycle_found)
1.1202 - length_bound = length_bound * BF_LIMIT_FACTOR / 100;
1.1203 + // Increase iteration limit if no cycle is found
1.1204 + if (!cycle_found) {
1.1205 + length_bound = static_cast<int>(length_bound * BF_LIMIT_FACTOR);
1.1206 + }
1.1207 }
1.1208 }
1.1209 }
1.1210
1.1211 - /// \brief Execute the algorithm using \ref Howard.
1.1212 - ///
1.1213 - /// Execute the algorithm using \ref Howard for negative
1.1214 - /// cycle detection.
1.1215 - void startMinMean() {
1.1216 - typedef Path<ResDigraph> ResPath;
1.1217 - Howard<ResDigraph, ResidualCostMap> mmc(*_res_graph, _res_cost);
1.1218 - ResPath cycle;
1.1219 + // Execute the "Minimum Mean Cycle Canceling" method
1.1220 + void startMinMeanCycleCanceling() {
1.1221 + typedef SimplePath<StaticDigraph> SPath;
1.1222 + typedef typename SPath::ArcIt SPathArcIt;
1.1223 + typedef typename Howard<StaticDigraph, CostArcMap>
1.1224 + ::template SetPath<SPath>::Create MMC;
1.1225 +
1.1226 + SPath cycle;
1.1227 + MMC mmc(_sgr, _cost_map);
1.1228 + mmc.cycle(cycle);
1.1229 + buildResidualNetwork();
1.1230 + while (mmc.findMinMean() && mmc.cycleLength() < 0) {
1.1231 + // Find the cycle
1.1232 + mmc.findCycle();
1.1233
1.1234 - mmc.cycle(cycle);
1.1235 - if (mmc.findMinMean()) {
1.1236 - while (mmc.cycleLength() < 0) {
1.1237 - // Finding the cycle
1.1238 - mmc.findCycle();
1.1239 + // Compute delta value
1.1240 + Value delta = INF;
1.1241 + for (SPathArcIt a(cycle); a != INVALID; ++a) {
1.1242 + Value d = _res_cap[_id_vec[_sgr.id(a)]];
1.1243 + if (d < delta) delta = d;
1.1244 + }
1.1245
1.1246 - // Finding the largest flow amount that can be augmented
1.1247 - // along the cycle
1.1248 - Capacity delta = 0;
1.1249 - for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e) {
1.1250 - if (delta == 0 || _res_graph->residualCapacity(e) < delta)
1.1251 - delta = _res_graph->residualCapacity(e);
1.1252 + // Augment along the cycle
1.1253 + for (SPathArcIt a(cycle); a != INVALID; ++a) {
1.1254 + int j = _id_vec[_sgr.id(a)];
1.1255 + _res_cap[j] -= delta;
1.1256 + _res_cap[_reverse[j]] += delta;
1.1257 + }
1.1258 +
1.1259 + // Rebuild the residual network
1.1260 + buildResidualNetwork();
1.1261 + }
1.1262 + }
1.1263 +
1.1264 + // Execute the "Cancel And Tighten" method
1.1265 + void startCancelAndTighten() {
1.1266 + // Constants for the min mean cycle computations
1.1267 + const double LIMIT_FACTOR = 1.0;
1.1268 + const int MIN_LIMIT = 5;
1.1269 +
1.1270 + // Contruct auxiliary data vectors
1.1271 + DoubleVector pi(_res_node_num, 0.0);
1.1272 + IntVector level(_res_node_num);
1.1273 + CharVector reached(_res_node_num);
1.1274 + CharVector processed(_res_node_num);
1.1275 + IntVector pred_node(_res_node_num);
1.1276 + IntVector pred_arc(_res_node_num);
1.1277 + std::vector<int> stack(_res_node_num);
1.1278 + std::vector<int> proc_vector(_res_node_num);
1.1279 +
1.1280 + // Initialize epsilon
1.1281 + double epsilon = 0;
1.1282 + for (int a = 0; a != _res_arc_num; ++a) {
1.1283 + if (_res_cap[a] > 0 && -_cost[a] > epsilon)
1.1284 + epsilon = -_cost[a];
1.1285 + }
1.1286 +
1.1287 + // Start phases
1.1288 + Tolerance<double> tol;
1.1289 + tol.epsilon(1e-6);
1.1290 + int limit = int(LIMIT_FACTOR * std::sqrt(double(_res_node_num)));
1.1291 + if (limit < MIN_LIMIT) limit = MIN_LIMIT;
1.1292 + int iter = limit;
1.1293 + while (epsilon * _res_node_num >= 1) {
1.1294 + // Find and cancel cycles in the admissible network using DFS
1.1295 + for (int u = 0; u != _res_node_num; ++u) {
1.1296 + reached[u] = false;
1.1297 + processed[u] = false;
1.1298 + }
1.1299 + int stack_head = -1;
1.1300 + int proc_head = -1;
1.1301 + for (int start = 0; start != _res_node_num; ++start) {
1.1302 + if (reached[start]) continue;
1.1303 +
1.1304 + // New start node
1.1305 + reached[start] = true;
1.1306 + pred_arc[start] = -1;
1.1307 + pred_node[start] = -1;
1.1308 +
1.1309 + // Find the first admissible outgoing arc
1.1310 + double p = pi[start];
1.1311 + int a = _first_out[start];
1.1312 + int last_out = _first_out[start+1];
1.1313 + for (; a != last_out && (_res_cap[a] == 0 ||
1.1314 + !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
1.1315 + if (a == last_out) {
1.1316 + processed[start] = true;
1.1317 + proc_vector[++proc_head] = start;
1.1318 + continue;
1.1319 + }
1.1320 + stack[++stack_head] = a;
1.1321 +
1.1322 + while (stack_head >= 0) {
1.1323 + int sa = stack[stack_head];
1.1324 + int u = _source[sa];
1.1325 + int v = _target[sa];
1.1326 +
1.1327 + if (!reached[v]) {
1.1328 + // A new node is reached
1.1329 + reached[v] = true;
1.1330 + pred_node[v] = u;
1.1331 + pred_arc[v] = sa;
1.1332 + p = pi[v];
1.1333 + a = _first_out[v];
1.1334 + last_out = _first_out[v+1];
1.1335 + for (; a != last_out && (_res_cap[a] == 0 ||
1.1336 + !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
1.1337 + stack[++stack_head] = a == last_out ? -1 : a;
1.1338 + } else {
1.1339 + if (!processed[v]) {
1.1340 + // A cycle is found
1.1341 + int n, w = u;
1.1342 + Value d, delta = _res_cap[sa];
1.1343 + for (n = u; n != v; n = pred_node[n]) {
1.1344 + d = _res_cap[pred_arc[n]];
1.1345 + if (d <= delta) {
1.1346 + delta = d;
1.1347 + w = pred_node[n];
1.1348 + }
1.1349 + }
1.1350 +
1.1351 + // Augment along the cycle
1.1352 + _res_cap[sa] -= delta;
1.1353 + _res_cap[_reverse[sa]] += delta;
1.1354 + for (n = u; n != v; n = pred_node[n]) {
1.1355 + int pa = pred_arc[n];
1.1356 + _res_cap[pa] -= delta;
1.1357 + _res_cap[_reverse[pa]] += delta;
1.1358 + }
1.1359 + for (n = u; stack_head > 0 && n != w; n = pred_node[n]) {
1.1360 + --stack_head;
1.1361 + reached[n] = false;
1.1362 + }
1.1363 + u = w;
1.1364 + }
1.1365 + v = u;
1.1366 +
1.1367 + // Find the next admissible outgoing arc
1.1368 + p = pi[v];
1.1369 + a = stack[stack_head] + 1;
1.1370 + last_out = _first_out[v+1];
1.1371 + for (; a != last_out && (_res_cap[a] == 0 ||
1.1372 + !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
1.1373 + stack[stack_head] = a == last_out ? -1 : a;
1.1374 + }
1.1375 +
1.1376 + while (stack_head >= 0 && stack[stack_head] == -1) {
1.1377 + processed[v] = true;
1.1378 + proc_vector[++proc_head] = v;
1.1379 + if (--stack_head >= 0) {
1.1380 + // Find the next admissible outgoing arc
1.1381 + v = _source[stack[stack_head]];
1.1382 + p = pi[v];
1.1383 + a = stack[stack_head] + 1;
1.1384 + last_out = _first_out[v+1];
1.1385 + for (; a != last_out && (_res_cap[a] == 0 ||
1.1386 + !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
1.1387 + stack[stack_head] = a == last_out ? -1 : a;
1.1388 + }
1.1389 + }
1.1390 + }
1.1391 + }
1.1392 +
1.1393 + // Tighten potentials and epsilon
1.1394 + if (--iter > 0) {
1.1395 + for (int u = 0; u != _res_node_num; ++u) {
1.1396 + level[u] = 0;
1.1397 + }
1.1398 + for (int i = proc_head; i > 0; --i) {
1.1399 + int u = proc_vector[i];
1.1400 + double p = pi[u];
1.1401 + int l = level[u] + 1;
1.1402 + int last_out = _first_out[u+1];
1.1403 + for (int a = _first_out[u]; a != last_out; ++a) {
1.1404 + int v = _target[a];
1.1405 + if (_res_cap[a] > 0 && tol.negative(_cost[a] + p - pi[v]) &&
1.1406 + l > level[v]) level[v] = l;
1.1407 + }
1.1408 }
1.1409
1.1410 - // Augmenting along the cycle
1.1411 - for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e)
1.1412 - _res_graph->augment(e, delta);
1.1413 + // Modify potentials
1.1414 + double q = std::numeric_limits<double>::max();
1.1415 + for (int u = 0; u != _res_node_num; ++u) {
1.1416 + int lu = level[u];
1.1417 + double p, pu = pi[u];
1.1418 + int last_out = _first_out[u+1];
1.1419 + for (int a = _first_out[u]; a != last_out; ++a) {
1.1420 + if (_res_cap[a] == 0) continue;
1.1421 + int v = _target[a];
1.1422 + int ld = lu - level[v];
1.1423 + if (ld > 0) {
1.1424 + p = (_cost[a] + pu - pi[v] + epsilon) / (ld + 1);
1.1425 + if (p < q) q = p;
1.1426 + }
1.1427 + }
1.1428 + }
1.1429 + for (int u = 0; u != _res_node_num; ++u) {
1.1430 + pi[u] -= q * level[u];
1.1431 + }
1.1432
1.1433 - // Finding the minimum cycle mean for the modified residual
1.1434 - // digraph
1.1435 - if (!mmc.findMinMean()) break;
1.1436 + // Modify epsilon
1.1437 + epsilon = 0;
1.1438 + for (int u = 0; u != _res_node_num; ++u) {
1.1439 + double curr, pu = pi[u];
1.1440 + int last_out = _first_out[u+1];
1.1441 + for (int a = _first_out[u]; a != last_out; ++a) {
1.1442 + if (_res_cap[a] == 0) continue;
1.1443 + curr = _cost[a] + pu - pi[_target[a]];
1.1444 + if (-curr > epsilon) epsilon = -curr;
1.1445 + }
1.1446 + }
1.1447 + } else {
1.1448 + typedef Howard<StaticDigraph, CostArcMap> MMC;
1.1449 + typedef typename BellmanFord<StaticDigraph, CostArcMap>
1.1450 + ::template SetDistMap<CostNodeMap>::Create BF;
1.1451 +
1.1452 + // Set epsilon to the minimum cycle mean
1.1453 + buildResidualNetwork();
1.1454 + MMC mmc(_sgr, _cost_map);
1.1455 + mmc.findMinMean();
1.1456 + epsilon = -mmc.cycleMean();
1.1457 + Cost cycle_cost = mmc.cycleLength();
1.1458 + int cycle_size = mmc.cycleArcNum();
1.1459 +
1.1460 + // Compute feasible potentials for the current epsilon
1.1461 + for (int i = 0; i != int(_cost_vec.size()); ++i) {
1.1462 + _cost_vec[i] = cycle_size * _cost_vec[i] - cycle_cost;
1.1463 + }
1.1464 + BF bf(_sgr, _cost_map);
1.1465 + bf.distMap(_pi_map);
1.1466 + bf.init(0);
1.1467 + bf.start();
1.1468 + for (int u = 0; u != _res_node_num; ++u) {
1.1469 + pi[u] = static_cast<double>(_pi[u]) / cycle_size;
1.1470 + }
1.1471 +
1.1472 + iter = limit;
1.1473 }
1.1474 }
1.1475 -
1.1476 - // Computing node potentials
1.1477 - BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost);
1.1478 - bf.init(0); bf.start();
1.1479 - for (NodeIt n(_graph); n != INVALID; ++n)
1.1480 - (*_potential)[n] = bf.dist(n);
1.1481 }
1.1482
1.1483 }; //class CycleCanceling