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