1.1 --- a/lemon/Makefile.am Fri Nov 13 00:09:35 2009 +0100
1.2 +++ b/lemon/Makefile.am Fri Nov 13 00:10:33 2009 +0100
1.3 @@ -62,7 +62,6 @@
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 \
2.1 --- a/lemon/cancel_and_tighten.h Fri Nov 13 00:09:35 2009 +0100
2.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
2.3 @@ -1,797 +0,0 @@
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 --- a/lemon/cycle_canceling.h Fri Nov 13 00:09:35 2009 +0100
3.2 +++ b/lemon/cycle_canceling.h Fri Nov 13 00:10:33 2009 +0100
3.3 @@ -19,441 +19,817 @@
3.4 #ifndef LEMON_CYCLE_CANCELING_H
3.5 #define LEMON_CYCLE_CANCELING_H
3.6
3.7 -/// \ingroup min_cost_flow
3.8 -///
3.9 +/// \ingroup min_cost_flow_algs
3.10 /// \file
3.11 -/// \brief Cycle-canceling algorithm for finding a minimum cost flow.
3.12 +/// \brief Cycle-canceling algorithms for finding a minimum cost flow.
3.13
3.14 #include <vector>
3.15 +#include <limits>
3.16 +
3.17 +#include <lemon/core.h>
3.18 +#include <lemon/maps.h>
3.19 +#include <lemon/path.h>
3.20 +#include <lemon/math.h>
3.21 +#include <lemon/static_graph.h>
3.22 #include <lemon/adaptors.h>
3.23 -#include <lemon/path.h>
3.24 -
3.25 #include <lemon/circulation.h>
3.26 #include <lemon/bellman_ford.h>
3.27 #include <lemon/howard.h>
3.28
3.29 namespace lemon {
3.30
3.31 - /// \addtogroup min_cost_flow
3.32 + /// \addtogroup min_cost_flow_algs
3.33 /// @{
3.34
3.35 - /// \brief Implementation of a cycle-canceling algorithm for
3.36 - /// finding a minimum cost flow.
3.37 + /// \brief Implementation of cycle-canceling algorithms for
3.38 + /// finding a \ref min_cost_flow "minimum cost flow".
3.39 ///
3.40 - /// \ref CycleCanceling implements a cycle-canceling algorithm for
3.41 - /// finding a minimum cost flow.
3.42 + /// \ref CycleCanceling implements three different cycle-canceling
3.43 + /// algorithms for finding a \ref min_cost_flow "minimum cost flow".
3.44 + /// The most efficent one (both theoretically and practically)
3.45 + /// is the \ref CANCEL_AND_TIGHTEN "Cancel and Tighten" algorithm,
3.46 + /// thus it is the default method.
3.47 + /// It is strongly polynomial, but in practice, it is typically much
3.48 + /// slower than the scaling algorithms and NetworkSimplex.
3.49 ///
3.50 - /// \tparam Digraph The digraph type the algorithm runs on.
3.51 - /// \tparam LowerMap The type of the lower bound map.
3.52 - /// \tparam CapacityMap The type of the capacity (upper bound) map.
3.53 - /// \tparam CostMap The type of the cost (length) map.
3.54 - /// \tparam SupplyMap The type of the supply map.
3.55 + /// Most of the parameters of the problem (except for the digraph)
3.56 + /// can be given using separate functions, and the algorithm can be
3.57 + /// executed using the \ref run() function. If some parameters are not
3.58 + /// specified, then default values will be used.
3.59 ///
3.60 - /// \warning
3.61 - /// - Arc capacities and costs should be \e non-negative \e integers.
3.62 - /// - Supply values should be \e signed \e integers.
3.63 - /// - The value types of the maps should be convertible to each other.
3.64 - /// - \c CostMap::Value must be signed type.
3.65 + /// \tparam GR The digraph type the algorithm runs on.
3.66 + /// \tparam V The number type used for flow amounts, capacity bounds
3.67 + /// and supply values in the algorithm. By default, it is \c int.
3.68 + /// \tparam C The number type used for costs and potentials in the
3.69 + /// algorithm. By default, it is the same as \c V.
3.70 ///
3.71 - /// \note By default the \ref BellmanFord "Bellman-Ford" algorithm is
3.72 - /// used for negative cycle detection with limited iteration number.
3.73 - /// However \ref CycleCanceling also provides the "Minimum Mean
3.74 - /// Cycle-Canceling" algorithm, which is \e strongly \e polynomial,
3.75 - /// but rather slower in practice.
3.76 - /// To use this version of the algorithm, call \ref run() with \c true
3.77 - /// parameter.
3.78 + /// \warning Both number types must be signed and all input data must
3.79 + /// be integer.
3.80 + /// \warning This algorithm does not support negative costs for such
3.81 + /// arcs that have infinite upper bound.
3.82 ///
3.83 - /// \author Peter Kovacs
3.84 - template < typename Digraph,
3.85 - typename LowerMap = typename Digraph::template ArcMap<int>,
3.86 - typename CapacityMap = typename Digraph::template ArcMap<int>,
3.87 - typename CostMap = typename Digraph::template ArcMap<int>,
3.88 - typename SupplyMap = typename Digraph::template NodeMap<int> >
3.89 + /// \note For more information about the three available methods,
3.90 + /// see \ref Method.
3.91 +#ifdef DOXYGEN
3.92 + template <typename GR, typename V, typename C>
3.93 +#else
3.94 + template <typename GR, typename V = int, typename C = V>
3.95 +#endif
3.96 class CycleCanceling
3.97 {
3.98 - TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
3.99 + public:
3.100
3.101 - typedef typename CapacityMap::Value Capacity;
3.102 - typedef typename CostMap::Value Cost;
3.103 - typedef typename SupplyMap::Value Supply;
3.104 - typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
3.105 - typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
3.106 -
3.107 - typedef ResidualDigraph< const Digraph,
3.108 - CapacityArcMap, CapacityArcMap > ResDigraph;
3.109 - typedef typename ResDigraph::Node ResNode;
3.110 - typedef typename ResDigraph::NodeIt ResNodeIt;
3.111 - typedef typename ResDigraph::Arc ResArc;
3.112 - typedef typename ResDigraph::ArcIt ResArcIt;
3.113 + /// The type of the digraph
3.114 + typedef GR Digraph;
3.115 + /// The type of the flow amounts, capacity bounds and supply values
3.116 + typedef V Value;
3.117 + /// The type of the arc costs
3.118 + typedef C Cost;
3.119
3.120 public:
3.121
3.122 - /// The type of the flow map.
3.123 - typedef typename Digraph::template ArcMap<Capacity> FlowMap;
3.124 - /// The type of the potential map.
3.125 - typedef typename Digraph::template NodeMap<Cost> PotentialMap;
3.126 + /// \brief Problem type constants for the \c run() function.
3.127 + ///
3.128 + /// Enum type containing the problem type constants that can be
3.129 + /// returned by the \ref run() function of the algorithm.
3.130 + enum ProblemType {
3.131 + /// The problem has no feasible solution (flow).
3.132 + INFEASIBLE,
3.133 + /// The problem has optimal solution (i.e. it is feasible and
3.134 + /// bounded), and the algorithm has found optimal flow and node
3.135 + /// potentials (primal and dual solutions).
3.136 + OPTIMAL,
3.137 + /// The digraph contains an arc of negative cost and infinite
3.138 + /// upper bound. It means that the objective function is unbounded
3.139 + /// on that arc, however, note that it could actually be bounded
3.140 + /// over the feasible flows, but this algroithm cannot handle
3.141 + /// these cases.
3.142 + UNBOUNDED
3.143 + };
3.144 +
3.145 + /// \brief Constants for selecting the used method.
3.146 + ///
3.147 + /// Enum type containing constants for selecting the used method
3.148 + /// for the \ref run() function.
3.149 + ///
3.150 + /// \ref CycleCanceling provides three different cycle-canceling
3.151 + /// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel and Tighten"
3.152 + /// is used, which proved to be the most efficient and the most robust
3.153 + /// on various test inputs.
3.154 + /// However, the other methods can be selected using the \ref run()
3.155 + /// function with the proper parameter.
3.156 + enum Method {
3.157 + /// A simple cycle-canceling method, which uses the
3.158 + /// \ref BellmanFord "Bellman-Ford" algorithm with limited iteration
3.159 + /// number for detecting negative cycles in the residual network.
3.160 + SIMPLE_CYCLE_CANCELING,
3.161 + /// The "Minimum Mean Cycle-Canceling" algorithm, which is a
3.162 + /// well-known strongly polynomial method. It improves along a
3.163 + /// \ref min_mean_cycle "minimum mean cycle" in each iteration.
3.164 + /// Its running time complexity is O(n<sup>2</sup>m<sup>3</sup>log(n)).
3.165 + MINIMUM_MEAN_CYCLE_CANCELING,
3.166 + /// The "Cancel And Tighten" algorithm, which can be viewed as an
3.167 + /// improved version of the previous method.
3.168 + /// It is faster both in theory and in practice, its running time
3.169 + /// complexity is O(n<sup>2</sup>m<sup>2</sup>log(n)).
3.170 + CANCEL_AND_TIGHTEN
3.171 + };
3.172
3.173 private:
3.174
3.175 - /// \brief Map adaptor class for handling residual arc costs.
3.176 - ///
3.177 - /// Map adaptor class for handling residual arc costs.
3.178 - class ResidualCostMap : public MapBase<ResArc, Cost>
3.179 - {
3.180 - private:
3.181 + TEMPLATE_DIGRAPH_TYPEDEFS(GR);
3.182 +
3.183 + typedef std::vector<int> IntVector;
3.184 + typedef std::vector<char> CharVector;
3.185 + typedef std::vector<double> DoubleVector;
3.186 + typedef std::vector<Value> ValueVector;
3.187 + typedef std::vector<Cost> CostVector;
3.188
3.189 - const CostMap &_cost_map;
3.190 -
3.191 + private:
3.192 +
3.193 + template <typename KT, typename VT>
3.194 + class VectorMap {
3.195 public:
3.196 -
3.197 - ///\e
3.198 - ResidualCostMap(const CostMap &cost_map) : _cost_map(cost_map) {}
3.199 -
3.200 - ///\e
3.201 - Cost operator[](const ResArc &e) const {
3.202 - return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
3.203 + typedef KT Key;
3.204 + typedef VT Value;
3.205 +
3.206 + VectorMap(std::vector<Value>& v) : _v(v) {}
3.207 +
3.208 + const Value& operator[](const Key& key) const {
3.209 + return _v[StaticDigraph::id(key)];
3.210 }
3.211
3.212 - }; //class ResidualCostMap
3.213 + Value& operator[](const Key& key) {
3.214 + return _v[StaticDigraph::id(key)];
3.215 + }
3.216 +
3.217 + void set(const Key& key, const Value& val) {
3.218 + _v[StaticDigraph::id(key)] = val;
3.219 + }
3.220 +
3.221 + private:
3.222 + std::vector<Value>& _v;
3.223 + };
3.224 +
3.225 + typedef VectorMap<StaticDigraph::Node, Cost> CostNodeMap;
3.226 + typedef VectorMap<StaticDigraph::Arc, Cost> CostArcMap;
3.227
3.228 private:
3.229
3.230 - // The maximum number of iterations for the first execution of the
3.231 - // Bellman-Ford algorithm. It should be at least 2.
3.232 - static const int BF_FIRST_LIMIT = 2;
3.233 - // The iteration limit for the Bellman-Ford algorithm is multiplied
3.234 - // by BF_LIMIT_FACTOR/100 in every round.
3.235 - static const int BF_LIMIT_FACTOR = 150;
3.236
3.237 - private:
3.238 + // Data related to the underlying digraph
3.239 + const GR &_graph;
3.240 + int _node_num;
3.241 + int _arc_num;
3.242 + int _res_node_num;
3.243 + int _res_arc_num;
3.244 + int _root;
3.245
3.246 - // The digraph the algorithm runs on
3.247 - const Digraph &_graph;
3.248 - // The original lower bound map
3.249 - const LowerMap *_lower;
3.250 - // The modified capacity map
3.251 - CapacityArcMap _capacity;
3.252 - // The original cost map
3.253 - const CostMap &_cost;
3.254 - // The modified supply map
3.255 - SupplyNodeMap _supply;
3.256 - bool _valid_supply;
3.257 + // Parameters of the problem
3.258 + bool _have_lower;
3.259 + Value _sum_supply;
3.260
3.261 - // Arc map of the current flow
3.262 - FlowMap *_flow;
3.263 - bool _local_flow;
3.264 - // Node map of the current potentials
3.265 - PotentialMap *_potential;
3.266 - bool _local_potential;
3.267 + // Data structures for storing the digraph
3.268 + IntNodeMap _node_id;
3.269 + IntArcMap _arc_idf;
3.270 + IntArcMap _arc_idb;
3.271 + IntVector _first_out;
3.272 + CharVector _forward;
3.273 + IntVector _source;
3.274 + IntVector _target;
3.275 + IntVector _reverse;
3.276
3.277 - // The residual digraph
3.278 - ResDigraph *_res_graph;
3.279 - // The residual cost map
3.280 - ResidualCostMap _res_cost;
3.281 + // Node and arc data
3.282 + ValueVector _lower;
3.283 + ValueVector _upper;
3.284 + CostVector _cost;
3.285 + ValueVector _supply;
3.286 +
3.287 + ValueVector _res_cap;
3.288 + CostVector _pi;
3.289 +
3.290 + // Data for a StaticDigraph structure
3.291 + typedef std::pair<int, int> IntPair;
3.292 + StaticDigraph _sgr;
3.293 + std::vector<IntPair> _arc_vec;
3.294 + std::vector<Cost> _cost_vec;
3.295 + IntVector _id_vec;
3.296 + CostArcMap _cost_map;
3.297 + CostNodeMap _pi_map;
3.298 +
3.299 + public:
3.300 +
3.301 + /// \brief Constant for infinite upper bounds (capacities).
3.302 + ///
3.303 + /// Constant for infinite upper bounds (capacities).
3.304 + /// It is \c std::numeric_limits<Value>::infinity() if available,
3.305 + /// \c std::numeric_limits<Value>::max() otherwise.
3.306 + const Value INF;
3.307
3.308 public:
3.309
3.310 - /// \brief General constructor (with lower bounds).
3.311 + /// \brief Constructor.
3.312 ///
3.313 - /// General constructor (with lower bounds).
3.314 + /// The constructor of the class.
3.315 ///
3.316 - /// \param digraph The digraph the algorithm runs on.
3.317 - /// \param lower The lower bounds of the arcs.
3.318 - /// \param capacity The capacities (upper bounds) of the arcs.
3.319 - /// \param cost The cost (length) values of the arcs.
3.320 - /// \param supply The supply values of the nodes (signed).
3.321 - CycleCanceling( const Digraph &digraph,
3.322 - const LowerMap &lower,
3.323 - const CapacityMap &capacity,
3.324 - const CostMap &cost,
3.325 - const SupplyMap &supply ) :
3.326 - _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost),
3.327 - _supply(digraph), _flow(NULL), _local_flow(false),
3.328 - _potential(NULL), _local_potential(false),
3.329 - _res_graph(NULL), _res_cost(_cost)
3.330 + /// \param graph The digraph the algorithm runs on.
3.331 + CycleCanceling(const GR& graph) :
3.332 + _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
3.333 + _cost_map(_cost_vec), _pi_map(_pi),
3.334 + INF(std::numeric_limits<Value>::has_infinity ?
3.335 + std::numeric_limits<Value>::infinity() :
3.336 + std::numeric_limits<Value>::max())
3.337 {
3.338 - // Check the sum of supply values
3.339 - Supply sum = 0;
3.340 - for (NodeIt n(_graph); n != INVALID; ++n) {
3.341 - _supply[n] = supply[n];
3.342 - sum += _supply[n];
3.343 + // Check the number types
3.344 + LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
3.345 + "The flow type of CycleCanceling must be signed");
3.346 + LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
3.347 + "The cost type of CycleCanceling must be signed");
3.348 +
3.349 + // Resize vectors
3.350 + _node_num = countNodes(_graph);
3.351 + _arc_num = countArcs(_graph);
3.352 + _res_node_num = _node_num + 1;
3.353 + _res_arc_num = 2 * (_arc_num + _node_num);
3.354 + _root = _node_num;
3.355 +
3.356 + _first_out.resize(_res_node_num + 1);
3.357 + _forward.resize(_res_arc_num);
3.358 + _source.resize(_res_arc_num);
3.359 + _target.resize(_res_arc_num);
3.360 + _reverse.resize(_res_arc_num);
3.361 +
3.362 + _lower.resize(_res_arc_num);
3.363 + _upper.resize(_res_arc_num);
3.364 + _cost.resize(_res_arc_num);
3.365 + _supply.resize(_res_node_num);
3.366 +
3.367 + _res_cap.resize(_res_arc_num);
3.368 + _pi.resize(_res_node_num);
3.369 +
3.370 + _arc_vec.reserve(_res_arc_num);
3.371 + _cost_vec.reserve(_res_arc_num);
3.372 + _id_vec.reserve(_res_arc_num);
3.373 +
3.374 + // Copy the graph
3.375 + int i = 0, j = 0, k = 2 * _arc_num + _node_num;
3.376 + for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
3.377 + _node_id[n] = i;
3.378 }
3.379 - _valid_supply = sum == 0;
3.380 -
3.381 - // Remove non-zero lower bounds
3.382 - for (ArcIt e(_graph); e != INVALID; ++e) {
3.383 - _capacity[e] = capacity[e];
3.384 - if (lower[e] != 0) {
3.385 - _capacity[e] -= lower[e];
3.386 - _supply[_graph.source(e)] -= lower[e];
3.387 - _supply[_graph.target(e)] += lower[e];
3.388 + i = 0;
3.389 + for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
3.390 + _first_out[i] = j;
3.391 + for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
3.392 + _arc_idf[a] = j;
3.393 + _forward[j] = true;
3.394 + _source[j] = i;
3.395 + _target[j] = _node_id[_graph.runningNode(a)];
3.396 }
3.397 + for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
3.398 + _arc_idb[a] = j;
3.399 + _forward[j] = false;
3.400 + _source[j] = i;
3.401 + _target[j] = _node_id[_graph.runningNode(a)];
3.402 + }
3.403 + _forward[j] = false;
3.404 + _source[j] = i;
3.405 + _target[j] = _root;
3.406 + _reverse[j] = k;
3.407 + _forward[k] = true;
3.408 + _source[k] = _root;
3.409 + _target[k] = i;
3.410 + _reverse[k] = j;
3.411 + ++j; ++k;
3.412 }
3.413 - }
3.414 -/*
3.415 - /// \brief General constructor (without lower bounds).
3.416 - ///
3.417 - /// General constructor (without lower bounds).
3.418 - ///
3.419 - /// \param digraph The digraph the algorithm runs on.
3.420 - /// \param capacity The capacities (upper bounds) of the arcs.
3.421 - /// \param cost The cost (length) values of the arcs.
3.422 - /// \param supply The supply values of the nodes (signed).
3.423 - CycleCanceling( const Digraph &digraph,
3.424 - const CapacityMap &capacity,
3.425 - const CostMap &cost,
3.426 - const SupplyMap &supply ) :
3.427 - _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
3.428 - _supply(supply), _flow(NULL), _local_flow(false),
3.429 - _potential(NULL), _local_potential(false), _res_graph(NULL),
3.430 - _res_cost(_cost)
3.431 - {
3.432 - // Check the sum of supply values
3.433 - Supply sum = 0;
3.434 - for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
3.435 - _valid_supply = sum == 0;
3.436 + _first_out[i] = j;
3.437 + _first_out[_res_node_num] = k;
3.438 + for (ArcIt a(_graph); a != INVALID; ++a) {
3.439 + int fi = _arc_idf[a];
3.440 + int bi = _arc_idb[a];
3.441 + _reverse[fi] = bi;
3.442 + _reverse[bi] = fi;
3.443 + }
3.444 +
3.445 + // Reset parameters
3.446 + reset();
3.447 }
3.448
3.449 - /// \brief Simple constructor (with lower bounds).
3.450 + /// \name Parameters
3.451 + /// The parameters of the algorithm can be specified using these
3.452 + /// functions.
3.453 +
3.454 + /// @{
3.455 +
3.456 + /// \brief Set the lower bounds on the arcs.
3.457 ///
3.458 - /// Simple constructor (with lower bounds).
3.459 + /// This function sets the lower bounds on the arcs.
3.460 + /// If it is not used before calling \ref run(), the lower bounds
3.461 + /// will be set to zero on all arcs.
3.462 ///
3.463 - /// \param digraph The digraph the algorithm runs on.
3.464 - /// \param lower The lower bounds of the arcs.
3.465 - /// \param capacity The capacities (upper bounds) of the arcs.
3.466 - /// \param cost The cost (length) values of the arcs.
3.467 - /// \param s The source node.
3.468 - /// \param t The target node.
3.469 - /// \param flow_value The required amount of flow from node \c s
3.470 - /// to node \c t (i.e. the supply of \c s and the demand of \c t).
3.471 - CycleCanceling( const Digraph &digraph,
3.472 - const LowerMap &lower,
3.473 - const CapacityMap &capacity,
3.474 - const CostMap &cost,
3.475 - Node s, Node t,
3.476 - Supply flow_value ) :
3.477 - _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost),
3.478 - _supply(digraph, 0), _flow(NULL), _local_flow(false),
3.479 - _potential(NULL), _local_potential(false), _res_graph(NULL),
3.480 - _res_cost(_cost)
3.481 - {
3.482 - // Remove non-zero lower bounds
3.483 - _supply[s] = flow_value;
3.484 - _supply[t] = -flow_value;
3.485 - for (ArcIt e(_graph); e != INVALID; ++e) {
3.486 - if (lower[e] != 0) {
3.487 - _capacity[e] -= lower[e];
3.488 - _supply[_graph.source(e)] -= lower[e];
3.489 - _supply[_graph.target(e)] += lower[e];
3.490 - }
3.491 + /// \param map An arc map storing the lower bounds.
3.492 + /// Its \c Value type must be convertible to the \c Value type
3.493 + /// of the algorithm.
3.494 + ///
3.495 + /// \return <tt>(*this)</tt>
3.496 + template <typename LowerMap>
3.497 + CycleCanceling& lowerMap(const LowerMap& map) {
3.498 + _have_lower = true;
3.499 + for (ArcIt a(_graph); a != INVALID; ++a) {
3.500 + _lower[_arc_idf[a]] = map[a];
3.501 + _lower[_arc_idb[a]] = map[a];
3.502 }
3.503 - _valid_supply = true;
3.504 - }
3.505 -
3.506 - /// \brief Simple constructor (without lower bounds).
3.507 - ///
3.508 - /// Simple constructor (without lower bounds).
3.509 - ///
3.510 - /// \param digraph The digraph the algorithm runs on.
3.511 - /// \param capacity The capacities (upper bounds) of the arcs.
3.512 - /// \param cost The cost (length) values of the arcs.
3.513 - /// \param s The source node.
3.514 - /// \param t The target node.
3.515 - /// \param flow_value The required amount of flow from node \c s
3.516 - /// to node \c t (i.e. the supply of \c s and the demand of \c t).
3.517 - CycleCanceling( const Digraph &digraph,
3.518 - const CapacityMap &capacity,
3.519 - const CostMap &cost,
3.520 - Node s, Node t,
3.521 - Supply flow_value ) :
3.522 - _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
3.523 - _supply(digraph, 0), _flow(NULL), _local_flow(false),
3.524 - _potential(NULL), _local_potential(false), _res_graph(NULL),
3.525 - _res_cost(_cost)
3.526 - {
3.527 - _supply[s] = flow_value;
3.528 - _supply[t] = -flow_value;
3.529 - _valid_supply = true;
3.530 - }
3.531 -*/
3.532 - /// Destructor.
3.533 - ~CycleCanceling() {
3.534 - if (_local_flow) delete _flow;
3.535 - if (_local_potential) delete _potential;
3.536 - delete _res_graph;
3.537 - }
3.538 -
3.539 - /// \brief Set the flow map.
3.540 - ///
3.541 - /// Set the flow map.
3.542 - ///
3.543 - /// \return \c (*this)
3.544 - CycleCanceling& flowMap(FlowMap &map) {
3.545 - if (_local_flow) {
3.546 - delete _flow;
3.547 - _local_flow = false;
3.548 - }
3.549 - _flow = ↦
3.550 return *this;
3.551 }
3.552
3.553 - /// \brief Set the potential map.
3.554 + /// \brief Set the upper bounds (capacities) on the arcs.
3.555 ///
3.556 - /// Set the potential map.
3.557 + /// This function sets the upper bounds (capacities) on the arcs.
3.558 + /// If it is not used before calling \ref run(), the upper bounds
3.559 + /// will be set to \ref INF on all arcs (i.e. the flow value will be
3.560 + /// unbounded from above).
3.561 ///
3.562 - /// \return \c (*this)
3.563 - CycleCanceling& potentialMap(PotentialMap &map) {
3.564 - if (_local_potential) {
3.565 - delete _potential;
3.566 - _local_potential = false;
3.567 + /// \param map An arc map storing the upper bounds.
3.568 + /// Its \c Value type must be convertible to the \c Value type
3.569 + /// of the algorithm.
3.570 + ///
3.571 + /// \return <tt>(*this)</tt>
3.572 + template<typename UpperMap>
3.573 + CycleCanceling& upperMap(const UpperMap& map) {
3.574 + for (ArcIt a(_graph); a != INVALID; ++a) {
3.575 + _upper[_arc_idf[a]] = map[a];
3.576 }
3.577 - _potential = ↦
3.578 return *this;
3.579 }
3.580
3.581 + /// \brief Set the costs of the arcs.
3.582 + ///
3.583 + /// This function sets the costs of the arcs.
3.584 + /// If it is not used before calling \ref run(), the costs
3.585 + /// will be set to \c 1 on all arcs.
3.586 + ///
3.587 + /// \param map An arc map storing the costs.
3.588 + /// Its \c Value type must be convertible to the \c Cost type
3.589 + /// of the algorithm.
3.590 + ///
3.591 + /// \return <tt>(*this)</tt>
3.592 + template<typename CostMap>
3.593 + CycleCanceling& costMap(const CostMap& map) {
3.594 + for (ArcIt a(_graph); a != INVALID; ++a) {
3.595 + _cost[_arc_idf[a]] = map[a];
3.596 + _cost[_arc_idb[a]] = -map[a];
3.597 + }
3.598 + return *this;
3.599 + }
3.600 +
3.601 + /// \brief Set the supply values of the nodes.
3.602 + ///
3.603 + /// This function sets the supply values of the nodes.
3.604 + /// If neither this function nor \ref stSupply() is used before
3.605 + /// calling \ref run(), the supply of each node will be set to zero.
3.606 + ///
3.607 + /// \param map A node map storing the supply values.
3.608 + /// Its \c Value type must be convertible to the \c Value type
3.609 + /// of the algorithm.
3.610 + ///
3.611 + /// \return <tt>(*this)</tt>
3.612 + template<typename SupplyMap>
3.613 + CycleCanceling& supplyMap(const SupplyMap& map) {
3.614 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.615 + _supply[_node_id[n]] = map[n];
3.616 + }
3.617 + return *this;
3.618 + }
3.619 +
3.620 + /// \brief Set single source and target nodes and a supply value.
3.621 + ///
3.622 + /// This function sets a single source node and a single target node
3.623 + /// and the required flow value.
3.624 + /// If neither this function nor \ref supplyMap() is used before
3.625 + /// calling \ref run(), the supply of each node will be set to zero.
3.626 + ///
3.627 + /// Using this function has the same effect as using \ref supplyMap()
3.628 + /// with such a map in which \c k is assigned to \c s, \c -k is
3.629 + /// assigned to \c t and all other nodes have zero supply value.
3.630 + ///
3.631 + /// \param s The source node.
3.632 + /// \param t The target node.
3.633 + /// \param k The required amount of flow from node \c s to node \c t
3.634 + /// (i.e. the supply of \c s and the demand of \c t).
3.635 + ///
3.636 + /// \return <tt>(*this)</tt>
3.637 + CycleCanceling& stSupply(const Node& s, const Node& t, Value k) {
3.638 + for (int i = 0; i != _res_node_num; ++i) {
3.639 + _supply[i] = 0;
3.640 + }
3.641 + _supply[_node_id[s]] = k;
3.642 + _supply[_node_id[t]] = -k;
3.643 + return *this;
3.644 + }
3.645 +
3.646 + /// @}
3.647 +
3.648 /// \name Execution control
3.649 + /// The algorithm can be executed using \ref run().
3.650
3.651 /// @{
3.652
3.653 /// \brief Run the algorithm.
3.654 ///
3.655 - /// Run the algorithm.
3.656 + /// This function runs the algorithm.
3.657 + /// The paramters can be specified using functions \ref lowerMap(),
3.658 + /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
3.659 + /// For example,
3.660 + /// \code
3.661 + /// CycleCanceling<ListDigraph> cc(graph);
3.662 + /// cc.lowerMap(lower).upperMap(upper).costMap(cost)
3.663 + /// .supplyMap(sup).run();
3.664 + /// \endcode
3.665 ///
3.666 - /// \param min_mean_cc Set this parameter to \c true to run the
3.667 - /// "Minimum Mean Cycle-Canceling" algorithm, which is strongly
3.668 - /// polynomial, but rather slower in practice.
3.669 + /// This function can be called more than once. All the parameters
3.670 + /// that have been given are kept for the next call, unless
3.671 + /// \ref reset() is called, thus only the modified parameters
3.672 + /// have to be set again. See \ref reset() for examples.
3.673 + /// However, the underlying digraph must not be modified after this
3.674 + /// class have been constructed, since it copies and extends the graph.
3.675 ///
3.676 - /// \return \c true if a feasible flow can be found.
3.677 - bool run(bool min_mean_cc = false) {
3.678 - return init() && start(min_mean_cc);
3.679 + /// \param method The cycle-canceling method that will be used.
3.680 + /// For more information, see \ref Method.
3.681 + ///
3.682 + /// \return \c INFEASIBLE if no feasible flow exists,
3.683 + /// \n \c OPTIMAL if the problem has optimal solution
3.684 + /// (i.e. it is feasible and bounded), and the algorithm has found
3.685 + /// optimal flow and node potentials (primal and dual solutions),
3.686 + /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
3.687 + /// and infinite upper bound. It means that the objective function
3.688 + /// is unbounded on that arc, however, note that it could actually be
3.689 + /// bounded over the feasible flows, but this algroithm cannot handle
3.690 + /// these cases.
3.691 + ///
3.692 + /// \see ProblemType, Method
3.693 + ProblemType run(Method method = CANCEL_AND_TIGHTEN) {
3.694 + ProblemType pt = init();
3.695 + if (pt != OPTIMAL) return pt;
3.696 + start(method);
3.697 + return OPTIMAL;
3.698 + }
3.699 +
3.700 + /// \brief Reset all the parameters that have been given before.
3.701 + ///
3.702 + /// This function resets all the paramaters that have been given
3.703 + /// before using functions \ref lowerMap(), \ref upperMap(),
3.704 + /// \ref costMap(), \ref supplyMap(), \ref stSupply().
3.705 + ///
3.706 + /// It is useful for multiple run() calls. If this function is not
3.707 + /// used, all the parameters given before are kept for the next
3.708 + /// \ref run() call.
3.709 + /// However, the underlying digraph must not be modified after this
3.710 + /// class have been constructed, since it copies and extends the graph.
3.711 + ///
3.712 + /// For example,
3.713 + /// \code
3.714 + /// CycleCanceling<ListDigraph> cs(graph);
3.715 + ///
3.716 + /// // First run
3.717 + /// cc.lowerMap(lower).upperMap(upper).costMap(cost)
3.718 + /// .supplyMap(sup).run();
3.719 + ///
3.720 + /// // Run again with modified cost map (reset() is not called,
3.721 + /// // so only the cost map have to be set again)
3.722 + /// cost[e] += 100;
3.723 + /// cc.costMap(cost).run();
3.724 + ///
3.725 + /// // Run again from scratch using reset()
3.726 + /// // (the lower bounds will be set to zero on all arcs)
3.727 + /// cc.reset();
3.728 + /// cc.upperMap(capacity).costMap(cost)
3.729 + /// .supplyMap(sup).run();
3.730 + /// \endcode
3.731 + ///
3.732 + /// \return <tt>(*this)</tt>
3.733 + CycleCanceling& reset() {
3.734 + for (int i = 0; i != _res_node_num; ++i) {
3.735 + _supply[i] = 0;
3.736 + }
3.737 + int limit = _first_out[_root];
3.738 + for (int j = 0; j != limit; ++j) {
3.739 + _lower[j] = 0;
3.740 + _upper[j] = INF;
3.741 + _cost[j] = _forward[j] ? 1 : -1;
3.742 + }
3.743 + for (int j = limit; j != _res_arc_num; ++j) {
3.744 + _lower[j] = 0;
3.745 + _upper[j] = INF;
3.746 + _cost[j] = 0;
3.747 + _cost[_reverse[j]] = 0;
3.748 + }
3.749 + _have_lower = false;
3.750 + return *this;
3.751 }
3.752
3.753 /// @}
3.754
3.755 /// \name Query Functions
3.756 - /// The result of the algorithm can be obtained using these
3.757 + /// The results of the algorithm can be obtained using these
3.758 /// functions.\n
3.759 - /// \ref lemon::CycleCanceling::run() "run()" must be called before
3.760 - /// using them.
3.761 + /// The \ref run() function must be called before using them.
3.762
3.763 /// @{
3.764
3.765 - /// \brief Return a const reference to the arc map storing the
3.766 - /// found flow.
3.767 + /// \brief Return the total cost of the found flow.
3.768 ///
3.769 - /// Return a const reference to the arc map storing the found flow.
3.770 + /// This function returns the total cost of the found flow.
3.771 + /// Its complexity is O(e).
3.772 + ///
3.773 + /// \note The return type of the function can be specified as a
3.774 + /// template parameter. For example,
3.775 + /// \code
3.776 + /// cc.totalCost<double>();
3.777 + /// \endcode
3.778 + /// It is useful if the total cost cannot be stored in the \c Cost
3.779 + /// type of the algorithm, which is the default return type of the
3.780 + /// function.
3.781 ///
3.782 /// \pre \ref run() must be called before using this function.
3.783 - const FlowMap& flowMap() const {
3.784 - return *_flow;
3.785 + template <typename Number>
3.786 + Number totalCost() const {
3.787 + Number c = 0;
3.788 + for (ArcIt a(_graph); a != INVALID; ++a) {
3.789 + int i = _arc_idb[a];
3.790 + c += static_cast<Number>(_res_cap[i]) *
3.791 + (-static_cast<Number>(_cost[i]));
3.792 + }
3.793 + return c;
3.794 }
3.795
3.796 - /// \brief Return a const reference to the node map storing the
3.797 - /// found potentials (the dual solution).
3.798 - ///
3.799 - /// Return a const reference to the node map storing the found
3.800 - /// potentials (the dual solution).
3.801 - ///
3.802 - /// \pre \ref run() must be called before using this function.
3.803 - const PotentialMap& potentialMap() const {
3.804 - return *_potential;
3.805 +#ifndef DOXYGEN
3.806 + Cost totalCost() const {
3.807 + return totalCost<Cost>();
3.808 }
3.809 +#endif
3.810
3.811 /// \brief Return the flow on the given arc.
3.812 ///
3.813 - /// Return the flow on the given arc.
3.814 + /// This function returns the flow on the given arc.
3.815 ///
3.816 /// \pre \ref run() must be called before using this function.
3.817 - Capacity flow(const Arc& arc) const {
3.818 - return (*_flow)[arc];
3.819 + Value flow(const Arc& a) const {
3.820 + return _res_cap[_arc_idb[a]];
3.821 }
3.822
3.823 - /// \brief Return the potential of the given node.
3.824 + /// \brief Return the flow map (the primal solution).
3.825 ///
3.826 - /// Return the potential of the given node.
3.827 + /// This function copies the flow value on each arc into the given
3.828 + /// map. The \c Value type of the algorithm must be convertible to
3.829 + /// the \c Value type of the map.
3.830 ///
3.831 /// \pre \ref run() must be called before using this function.
3.832 - Cost potential(const Node& node) const {
3.833 - return (*_potential)[node];
3.834 + template <typename FlowMap>
3.835 + void flowMap(FlowMap &map) const {
3.836 + for (ArcIt a(_graph); a != INVALID; ++a) {
3.837 + map.set(a, _res_cap[_arc_idb[a]]);
3.838 + }
3.839 }
3.840
3.841 - /// \brief Return the total cost of the found flow.
3.842 + /// \brief Return the potential (dual value) of the given node.
3.843 ///
3.844 - /// Return the total cost of the found flow. The complexity of the
3.845 - /// function is \f$ O(e) \f$.
3.846 + /// This function returns the potential (dual value) of the
3.847 + /// given node.
3.848 ///
3.849 /// \pre \ref run() must be called before using this function.
3.850 - Cost totalCost() const {
3.851 - Cost c = 0;
3.852 - for (ArcIt e(_graph); e != INVALID; ++e)
3.853 - c += (*_flow)[e] * _cost[e];
3.854 - return c;
3.855 + Cost potential(const Node& n) const {
3.856 + return static_cast<Cost>(_pi[_node_id[n]]);
3.857 + }
3.858 +
3.859 + /// \brief Return the potential map (the dual solution).
3.860 + ///
3.861 + /// This function copies the potential (dual value) of each node
3.862 + /// into the given map.
3.863 + /// The \c Cost type of the algorithm must be convertible to the
3.864 + /// \c Value type of the map.
3.865 + ///
3.866 + /// \pre \ref run() must be called before using this function.
3.867 + template <typename PotentialMap>
3.868 + void potentialMap(PotentialMap &map) const {
3.869 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.870 + map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
3.871 + }
3.872 }
3.873
3.874 /// @}
3.875
3.876 private:
3.877
3.878 - /// Initialize the algorithm.
3.879 - bool init() {
3.880 - if (!_valid_supply) return false;
3.881 + // Initialize the algorithm
3.882 + ProblemType init() {
3.883 + if (_res_node_num <= 1) return INFEASIBLE;
3.884
3.885 - // Initializing flow and potential maps
3.886 - if (!_flow) {
3.887 - _flow = new FlowMap(_graph);
3.888 - _local_flow = true;
3.889 + // Check the sum of supply values
3.890 + _sum_supply = 0;
3.891 + for (int i = 0; i != _root; ++i) {
3.892 + _sum_supply += _supply[i];
3.893 }
3.894 - if (!_potential) {
3.895 - _potential = new PotentialMap(_graph);
3.896 - _local_potential = true;
3.897 + if (_sum_supply > 0) return INFEASIBLE;
3.898 +
3.899 +
3.900 + // Initialize vectors
3.901 + for (int i = 0; i != _res_node_num; ++i) {
3.902 + _pi[i] = 0;
3.903 + }
3.904 + ValueVector excess(_supply);
3.905 +
3.906 + // Remove infinite upper bounds and check negative arcs
3.907 + const Value MAX = std::numeric_limits<Value>::max();
3.908 + int last_out;
3.909 + if (_have_lower) {
3.910 + for (int i = 0; i != _root; ++i) {
3.911 + last_out = _first_out[i+1];
3.912 + for (int j = _first_out[i]; j != last_out; ++j) {
3.913 + if (_forward[j]) {
3.914 + Value c = _cost[j] < 0 ? _upper[j] : _lower[j];
3.915 + if (c >= MAX) return UNBOUNDED;
3.916 + excess[i] -= c;
3.917 + excess[_target[j]] += c;
3.918 + }
3.919 + }
3.920 + }
3.921 + } else {
3.922 + for (int i = 0; i != _root; ++i) {
3.923 + last_out = _first_out[i+1];
3.924 + for (int j = _first_out[i]; j != last_out; ++j) {
3.925 + if (_forward[j] && _cost[j] < 0) {
3.926 + Value c = _upper[j];
3.927 + if (c >= MAX) return UNBOUNDED;
3.928 + excess[i] -= c;
3.929 + excess[_target[j]] += c;
3.930 + }
3.931 + }
3.932 + }
3.933 + }
3.934 + Value ex, max_cap = 0;
3.935 + for (int i = 0; i != _res_node_num; ++i) {
3.936 + ex = excess[i];
3.937 + if (ex < 0) max_cap -= ex;
3.938 + }
3.939 + for (int j = 0; j != _res_arc_num; ++j) {
3.940 + if (_upper[j] >= MAX) _upper[j] = max_cap;
3.941 }
3.942
3.943 - _res_graph = new ResDigraph(_graph, _capacity, *_flow);
3.944 + // Initialize maps for Circulation and remove non-zero lower bounds
3.945 + ConstMap<Arc, Value> low(0);
3.946 + typedef typename Digraph::template ArcMap<Value> ValueArcMap;
3.947 + typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
3.948 + ValueArcMap cap(_graph), flow(_graph);
3.949 + ValueNodeMap sup(_graph);
3.950 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.951 + sup[n] = _supply[_node_id[n]];
3.952 + }
3.953 + if (_have_lower) {
3.954 + for (ArcIt a(_graph); a != INVALID; ++a) {
3.955 + int j = _arc_idf[a];
3.956 + Value c = _lower[j];
3.957 + cap[a] = _upper[j] - c;
3.958 + sup[_graph.source(a)] -= c;
3.959 + sup[_graph.target(a)] += c;
3.960 + }
3.961 + } else {
3.962 + for (ArcIt a(_graph); a != INVALID; ++a) {
3.963 + cap[a] = _upper[_arc_idf[a]];
3.964 + }
3.965 + }
3.966
3.967 - // Finding a feasible flow using Circulation
3.968 - Circulation< Digraph, ConstMap<Arc, Capacity>, CapacityArcMap,
3.969 - SupplyMap >
3.970 - circulation( _graph, constMap<Arc>(Capacity(0)), _capacity,
3.971 - _supply );
3.972 - return circulation.flowMap(*_flow).run();
3.973 + // Find a feasible flow using Circulation
3.974 + Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
3.975 + circ(_graph, low, cap, sup);
3.976 + if (!circ.flowMap(flow).run()) return INFEASIBLE;
3.977 +
3.978 + // Set residual capacities and handle GEQ supply type
3.979 + if (_sum_supply < 0) {
3.980 + for (ArcIt a(_graph); a != INVALID; ++a) {
3.981 + Value fa = flow[a];
3.982 + _res_cap[_arc_idf[a]] = cap[a] - fa;
3.983 + _res_cap[_arc_idb[a]] = fa;
3.984 + sup[_graph.source(a)] -= fa;
3.985 + sup[_graph.target(a)] += fa;
3.986 + }
3.987 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.988 + excess[_node_id[n]] = sup[n];
3.989 + }
3.990 + for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
3.991 + int u = _target[a];
3.992 + int ra = _reverse[a];
3.993 + _res_cap[a] = -_sum_supply + 1;
3.994 + _res_cap[ra] = -excess[u];
3.995 + _cost[a] = 0;
3.996 + _cost[ra] = 0;
3.997 + }
3.998 + } else {
3.999 + for (ArcIt a(_graph); a != INVALID; ++a) {
3.1000 + Value fa = flow[a];
3.1001 + _res_cap[_arc_idf[a]] = cap[a] - fa;
3.1002 + _res_cap[_arc_idb[a]] = fa;
3.1003 + }
3.1004 + for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
3.1005 + int ra = _reverse[a];
3.1006 + _res_cap[a] = 1;
3.1007 + _res_cap[ra] = 0;
3.1008 + _cost[a] = 0;
3.1009 + _cost[ra] = 0;
3.1010 + }
3.1011 + }
3.1012 +
3.1013 + return OPTIMAL;
3.1014 + }
3.1015 +
3.1016 + // Build a StaticDigraph structure containing the current
3.1017 + // residual network
3.1018 + void buildResidualNetwork() {
3.1019 + _arc_vec.clear();
3.1020 + _cost_vec.clear();
3.1021 + _id_vec.clear();
3.1022 + for (int j = 0; j != _res_arc_num; ++j) {
3.1023 + if (_res_cap[j] > 0) {
3.1024 + _arc_vec.push_back(IntPair(_source[j], _target[j]));
3.1025 + _cost_vec.push_back(_cost[j]);
3.1026 + _id_vec.push_back(j);
3.1027 + }
3.1028 + }
3.1029 + _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
3.1030 }
3.1031
3.1032 - bool start(bool min_mean_cc) {
3.1033 - if (min_mean_cc)
3.1034 - startMinMean();
3.1035 - else
3.1036 - start();
3.1037 + // Execute the algorithm and transform the results
3.1038 + void start(Method method) {
3.1039 + // Execute the algorithm
3.1040 + switch (method) {
3.1041 + case SIMPLE_CYCLE_CANCELING:
3.1042 + startSimpleCycleCanceling();
3.1043 + break;
3.1044 + case MINIMUM_MEAN_CYCLE_CANCELING:
3.1045 + startMinMeanCycleCanceling();
3.1046 + break;
3.1047 + case CANCEL_AND_TIGHTEN:
3.1048 + startCancelAndTighten();
3.1049 + break;
3.1050 + }
3.1051
3.1052 - // Handling non-zero lower bounds
3.1053 - if (_lower) {
3.1054 - for (ArcIt e(_graph); e != INVALID; ++e)
3.1055 - (*_flow)[e] += (*_lower)[e];
3.1056 + // Compute node potentials
3.1057 + if (method != SIMPLE_CYCLE_CANCELING) {
3.1058 + buildResidualNetwork();
3.1059 + typename BellmanFord<StaticDigraph, CostArcMap>
3.1060 + ::template SetDistMap<CostNodeMap>::Create bf(_sgr, _cost_map);
3.1061 + bf.distMap(_pi_map);
3.1062 + bf.init(0);
3.1063 + bf.start();
3.1064 }
3.1065 - return true;
3.1066 +
3.1067 + // Handle non-zero lower bounds
3.1068 + if (_have_lower) {
3.1069 + int limit = _first_out[_root];
3.1070 + for (int j = 0; j != limit; ++j) {
3.1071 + if (!_forward[j]) _res_cap[j] += _lower[j];
3.1072 + }
3.1073 + }
3.1074 }
3.1075
3.1076 - /// \brief Execute the algorithm using \ref BellmanFord.
3.1077 - ///
3.1078 - /// Execute the algorithm using the \ref BellmanFord
3.1079 - /// "Bellman-Ford" algorithm for negative cycle detection with
3.1080 - /// successively larger limit for the number of iterations.
3.1081 - void start() {
3.1082 - typename BellmanFord<ResDigraph, ResidualCostMap>::PredMap pred(*_res_graph);
3.1083 - typename ResDigraph::template NodeMap<int> visited(*_res_graph);
3.1084 - std::vector<ResArc> cycle;
3.1085 - int node_num = countNodes(_graph);
3.1086 + // Execute the "Simple Cycle Canceling" method
3.1087 + void startSimpleCycleCanceling() {
3.1088 + // Constants for computing the iteration limits
3.1089 + const int BF_FIRST_LIMIT = 2;
3.1090 + const double BF_LIMIT_FACTOR = 1.5;
3.1091 +
3.1092 + typedef VectorMap<StaticDigraph::Arc, Value> FilterMap;
3.1093 + typedef FilterArcs<StaticDigraph, FilterMap> ResDigraph;
3.1094 + typedef VectorMap<StaticDigraph::Node, StaticDigraph::Arc> PredMap;
3.1095 + typedef typename BellmanFord<ResDigraph, CostArcMap>
3.1096 + ::template SetDistMap<CostNodeMap>
3.1097 + ::template SetPredMap<PredMap>::Create BF;
3.1098 +
3.1099 + // Build the residual network
3.1100 + _arc_vec.clear();
3.1101 + _cost_vec.clear();
3.1102 + for (int j = 0; j != _res_arc_num; ++j) {
3.1103 + _arc_vec.push_back(IntPair(_source[j], _target[j]));
3.1104 + _cost_vec.push_back(_cost[j]);
3.1105 + }
3.1106 + _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
3.1107 +
3.1108 + FilterMap filter_map(_res_cap);
3.1109 + ResDigraph rgr(_sgr, filter_map);
3.1110 + std::vector<int> cycle;
3.1111 + std::vector<StaticDigraph::Arc> pred(_res_arc_num);
3.1112 + PredMap pred_map(pred);
3.1113 + BF bf(rgr, _cost_map);
3.1114 + bf.distMap(_pi_map).predMap(pred_map);
3.1115
3.1116 int length_bound = BF_FIRST_LIMIT;
3.1117 bool optimal = false;
3.1118 while (!optimal) {
3.1119 - BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost);
3.1120 - bf.predMap(pred);
3.1121 bf.init(0);
3.1122 int iter_num = 0;
3.1123 bool cycle_found = false;
3.1124 while (!cycle_found) {
3.1125 - int curr_iter_num = iter_num + length_bound <= node_num ?
3.1126 - length_bound : node_num - iter_num;
3.1127 + // Perform some iterations of the Bellman-Ford algorithm
3.1128 + int curr_iter_num = iter_num + length_bound <= _node_num ?
3.1129 + length_bound : _node_num - iter_num;
3.1130 iter_num += curr_iter_num;
3.1131 int real_iter_num = curr_iter_num;
3.1132 for (int i = 0; i < curr_iter_num; ++i) {
3.1133 @@ -465,89 +841,290 @@
3.1134 if (real_iter_num < curr_iter_num) {
3.1135 // Optimal flow is found
3.1136 optimal = true;
3.1137 - // Setting node potentials
3.1138 - for (NodeIt n(_graph); n != INVALID; ++n)
3.1139 - (*_potential)[n] = bf.dist(n);
3.1140 break;
3.1141 } else {
3.1142 - // Searching for node disjoint negative cycles
3.1143 - for (ResNodeIt n(*_res_graph); n != INVALID; ++n)
3.1144 - visited[n] = 0;
3.1145 + // Search for node disjoint negative cycles
3.1146 + std::vector<int> state(_res_node_num, 0);
3.1147 int id = 0;
3.1148 - for (ResNodeIt n(*_res_graph); n != INVALID; ++n) {
3.1149 - if (visited[n] > 0) continue;
3.1150 - visited[n] = ++id;
3.1151 - ResNode u = pred[n] == INVALID ?
3.1152 - INVALID : _res_graph->source(pred[n]);
3.1153 - while (u != INVALID && visited[u] == 0) {
3.1154 - visited[u] = id;
3.1155 - u = pred[u] == INVALID ?
3.1156 - INVALID : _res_graph->source(pred[u]);
3.1157 + for (int u = 0; u != _res_node_num; ++u) {
3.1158 + if (state[u] != 0) continue;
3.1159 + ++id;
3.1160 + int v = u;
3.1161 + for (; v != -1 && state[v] == 0; v = pred[v] == INVALID ?
3.1162 + -1 : rgr.id(rgr.source(pred[v]))) {
3.1163 + state[v] = id;
3.1164 }
3.1165 - if (u != INVALID && visited[u] == id) {
3.1166 - // Finding the negative cycle
3.1167 + if (v != -1 && state[v] == id) {
3.1168 + // A negative cycle is found
3.1169 cycle_found = true;
3.1170 cycle.clear();
3.1171 - ResArc e = pred[u];
3.1172 - cycle.push_back(e);
3.1173 - Capacity d = _res_graph->residualCapacity(e);
3.1174 - while (_res_graph->source(e) != u) {
3.1175 - cycle.push_back(e = pred[_res_graph->source(e)]);
3.1176 - if (_res_graph->residualCapacity(e) < d)
3.1177 - d = _res_graph->residualCapacity(e);
3.1178 + StaticDigraph::Arc a = pred[v];
3.1179 + Value d, delta = _res_cap[rgr.id(a)];
3.1180 + cycle.push_back(rgr.id(a));
3.1181 + while (rgr.id(rgr.source(a)) != v) {
3.1182 + a = pred_map[rgr.source(a)];
3.1183 + d = _res_cap[rgr.id(a)];
3.1184 + if (d < delta) delta = d;
3.1185 + cycle.push_back(rgr.id(a));
3.1186 }
3.1187
3.1188 - // Augmenting along the cycle
3.1189 - for (int i = 0; i < int(cycle.size()); ++i)
3.1190 - _res_graph->augment(cycle[i], d);
3.1191 + // Augment along the cycle
3.1192 + for (int i = 0; i < int(cycle.size()); ++i) {
3.1193 + int j = cycle[i];
3.1194 + _res_cap[j] -= delta;
3.1195 + _res_cap[_reverse[j]] += delta;
3.1196 + }
3.1197 }
3.1198 }
3.1199 }
3.1200
3.1201 - if (!cycle_found)
3.1202 - length_bound = length_bound * BF_LIMIT_FACTOR / 100;
3.1203 + // Increase iteration limit if no cycle is found
3.1204 + if (!cycle_found) {
3.1205 + length_bound = static_cast<int>(length_bound * BF_LIMIT_FACTOR);
3.1206 + }
3.1207 }
3.1208 }
3.1209 }
3.1210
3.1211 - /// \brief Execute the algorithm using \ref Howard.
3.1212 - ///
3.1213 - /// Execute the algorithm using \ref Howard for negative
3.1214 - /// cycle detection.
3.1215 - void startMinMean() {
3.1216 - typedef Path<ResDigraph> ResPath;
3.1217 - Howard<ResDigraph, ResidualCostMap> mmc(*_res_graph, _res_cost);
3.1218 - ResPath cycle;
3.1219 + // Execute the "Minimum Mean Cycle Canceling" method
3.1220 + void startMinMeanCycleCanceling() {
3.1221 + typedef SimplePath<StaticDigraph> SPath;
3.1222 + typedef typename SPath::ArcIt SPathArcIt;
3.1223 + typedef typename Howard<StaticDigraph, CostArcMap>
3.1224 + ::template SetPath<SPath>::Create MMC;
3.1225 +
3.1226 + SPath cycle;
3.1227 + MMC mmc(_sgr, _cost_map);
3.1228 + mmc.cycle(cycle);
3.1229 + buildResidualNetwork();
3.1230 + while (mmc.findMinMean() && mmc.cycleLength() < 0) {
3.1231 + // Find the cycle
3.1232 + mmc.findCycle();
3.1233
3.1234 - mmc.cycle(cycle);
3.1235 - if (mmc.findMinMean()) {
3.1236 - while (mmc.cycleLength() < 0) {
3.1237 - // Finding the cycle
3.1238 - mmc.findCycle();
3.1239 + // Compute delta value
3.1240 + Value delta = INF;
3.1241 + for (SPathArcIt a(cycle); a != INVALID; ++a) {
3.1242 + Value d = _res_cap[_id_vec[_sgr.id(a)]];
3.1243 + if (d < delta) delta = d;
3.1244 + }
3.1245
3.1246 - // Finding the largest flow amount that can be augmented
3.1247 - // along the cycle
3.1248 - Capacity delta = 0;
3.1249 - for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e) {
3.1250 - if (delta == 0 || _res_graph->residualCapacity(e) < delta)
3.1251 - delta = _res_graph->residualCapacity(e);
3.1252 + // Augment along the cycle
3.1253 + for (SPathArcIt a(cycle); a != INVALID; ++a) {
3.1254 + int j = _id_vec[_sgr.id(a)];
3.1255 + _res_cap[j] -= delta;
3.1256 + _res_cap[_reverse[j]] += delta;
3.1257 + }
3.1258 +
3.1259 + // Rebuild the residual network
3.1260 + buildResidualNetwork();
3.1261 + }
3.1262 + }
3.1263 +
3.1264 + // Execute the "Cancel And Tighten" method
3.1265 + void startCancelAndTighten() {
3.1266 + // Constants for the min mean cycle computations
3.1267 + const double LIMIT_FACTOR = 1.0;
3.1268 + const int MIN_LIMIT = 5;
3.1269 +
3.1270 + // Contruct auxiliary data vectors
3.1271 + DoubleVector pi(_res_node_num, 0.0);
3.1272 + IntVector level(_res_node_num);
3.1273 + CharVector reached(_res_node_num);
3.1274 + CharVector processed(_res_node_num);
3.1275 + IntVector pred_node(_res_node_num);
3.1276 + IntVector pred_arc(_res_node_num);
3.1277 + std::vector<int> stack(_res_node_num);
3.1278 + std::vector<int> proc_vector(_res_node_num);
3.1279 +
3.1280 + // Initialize epsilon
3.1281 + double epsilon = 0;
3.1282 + for (int a = 0; a != _res_arc_num; ++a) {
3.1283 + if (_res_cap[a] > 0 && -_cost[a] > epsilon)
3.1284 + epsilon = -_cost[a];
3.1285 + }
3.1286 +
3.1287 + // Start phases
3.1288 + Tolerance<double> tol;
3.1289 + tol.epsilon(1e-6);
3.1290 + int limit = int(LIMIT_FACTOR * std::sqrt(double(_res_node_num)));
3.1291 + if (limit < MIN_LIMIT) limit = MIN_LIMIT;
3.1292 + int iter = limit;
3.1293 + while (epsilon * _res_node_num >= 1) {
3.1294 + // Find and cancel cycles in the admissible network using DFS
3.1295 + for (int u = 0; u != _res_node_num; ++u) {
3.1296 + reached[u] = false;
3.1297 + processed[u] = false;
3.1298 + }
3.1299 + int stack_head = -1;
3.1300 + int proc_head = -1;
3.1301 + for (int start = 0; start != _res_node_num; ++start) {
3.1302 + if (reached[start]) continue;
3.1303 +
3.1304 + // New start node
3.1305 + reached[start] = true;
3.1306 + pred_arc[start] = -1;
3.1307 + pred_node[start] = -1;
3.1308 +
3.1309 + // Find the first admissible outgoing arc
3.1310 + double p = pi[start];
3.1311 + int a = _first_out[start];
3.1312 + int last_out = _first_out[start+1];
3.1313 + for (; a != last_out && (_res_cap[a] == 0 ||
3.1314 + !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
3.1315 + if (a == last_out) {
3.1316 + processed[start] = true;
3.1317 + proc_vector[++proc_head] = start;
3.1318 + continue;
3.1319 + }
3.1320 + stack[++stack_head] = a;
3.1321 +
3.1322 + while (stack_head >= 0) {
3.1323 + int sa = stack[stack_head];
3.1324 + int u = _source[sa];
3.1325 + int v = _target[sa];
3.1326 +
3.1327 + if (!reached[v]) {
3.1328 + // A new node is reached
3.1329 + reached[v] = true;
3.1330 + pred_node[v] = u;
3.1331 + pred_arc[v] = sa;
3.1332 + p = pi[v];
3.1333 + a = _first_out[v];
3.1334 + last_out = _first_out[v+1];
3.1335 + for (; a != last_out && (_res_cap[a] == 0 ||
3.1336 + !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
3.1337 + stack[++stack_head] = a == last_out ? -1 : a;
3.1338 + } else {
3.1339 + if (!processed[v]) {
3.1340 + // A cycle is found
3.1341 + int n, w = u;
3.1342 + Value d, delta = _res_cap[sa];
3.1343 + for (n = u; n != v; n = pred_node[n]) {
3.1344 + d = _res_cap[pred_arc[n]];
3.1345 + if (d <= delta) {
3.1346 + delta = d;
3.1347 + w = pred_node[n];
3.1348 + }
3.1349 + }
3.1350 +
3.1351 + // Augment along the cycle
3.1352 + _res_cap[sa] -= delta;
3.1353 + _res_cap[_reverse[sa]] += delta;
3.1354 + for (n = u; n != v; n = pred_node[n]) {
3.1355 + int pa = pred_arc[n];
3.1356 + _res_cap[pa] -= delta;
3.1357 + _res_cap[_reverse[pa]] += delta;
3.1358 + }
3.1359 + for (n = u; stack_head > 0 && n != w; n = pred_node[n]) {
3.1360 + --stack_head;
3.1361 + reached[n] = false;
3.1362 + }
3.1363 + u = w;
3.1364 + }
3.1365 + v = u;
3.1366 +
3.1367 + // Find the next admissible outgoing arc
3.1368 + p = pi[v];
3.1369 + a = stack[stack_head] + 1;
3.1370 + last_out = _first_out[v+1];
3.1371 + for (; a != last_out && (_res_cap[a] == 0 ||
3.1372 + !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
3.1373 + stack[stack_head] = a == last_out ? -1 : a;
3.1374 + }
3.1375 +
3.1376 + while (stack_head >= 0 && stack[stack_head] == -1) {
3.1377 + processed[v] = true;
3.1378 + proc_vector[++proc_head] = v;
3.1379 + if (--stack_head >= 0) {
3.1380 + // Find the next admissible outgoing arc
3.1381 + v = _source[stack[stack_head]];
3.1382 + p = pi[v];
3.1383 + a = stack[stack_head] + 1;
3.1384 + last_out = _first_out[v+1];
3.1385 + for (; a != last_out && (_res_cap[a] == 0 ||
3.1386 + !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
3.1387 + stack[stack_head] = a == last_out ? -1 : a;
3.1388 + }
3.1389 + }
3.1390 + }
3.1391 + }
3.1392 +
3.1393 + // Tighten potentials and epsilon
3.1394 + if (--iter > 0) {
3.1395 + for (int u = 0; u != _res_node_num; ++u) {
3.1396 + level[u] = 0;
3.1397 + }
3.1398 + for (int i = proc_head; i > 0; --i) {
3.1399 + int u = proc_vector[i];
3.1400 + double p = pi[u];
3.1401 + int l = level[u] + 1;
3.1402 + int last_out = _first_out[u+1];
3.1403 + for (int a = _first_out[u]; a != last_out; ++a) {
3.1404 + int v = _target[a];
3.1405 + if (_res_cap[a] > 0 && tol.negative(_cost[a] + p - pi[v]) &&
3.1406 + l > level[v]) level[v] = l;
3.1407 + }
3.1408 }
3.1409
3.1410 - // Augmenting along the cycle
3.1411 - for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e)
3.1412 - _res_graph->augment(e, delta);
3.1413 + // Modify potentials
3.1414 + double q = std::numeric_limits<double>::max();
3.1415 + for (int u = 0; u != _res_node_num; ++u) {
3.1416 + int lu = level[u];
3.1417 + double p, pu = pi[u];
3.1418 + int last_out = _first_out[u+1];
3.1419 + for (int a = _first_out[u]; a != last_out; ++a) {
3.1420 + if (_res_cap[a] == 0) continue;
3.1421 + int v = _target[a];
3.1422 + int ld = lu - level[v];
3.1423 + if (ld > 0) {
3.1424 + p = (_cost[a] + pu - pi[v] + epsilon) / (ld + 1);
3.1425 + if (p < q) q = p;
3.1426 + }
3.1427 + }
3.1428 + }
3.1429 + for (int u = 0; u != _res_node_num; ++u) {
3.1430 + pi[u] -= q * level[u];
3.1431 + }
3.1432
3.1433 - // Finding the minimum cycle mean for the modified residual
3.1434 - // digraph
3.1435 - if (!mmc.findMinMean()) break;
3.1436 + // Modify epsilon
3.1437 + epsilon = 0;
3.1438 + for (int u = 0; u != _res_node_num; ++u) {
3.1439 + double curr, pu = pi[u];
3.1440 + int last_out = _first_out[u+1];
3.1441 + for (int a = _first_out[u]; a != last_out; ++a) {
3.1442 + if (_res_cap[a] == 0) continue;
3.1443 + curr = _cost[a] + pu - pi[_target[a]];
3.1444 + if (-curr > epsilon) epsilon = -curr;
3.1445 + }
3.1446 + }
3.1447 + } else {
3.1448 + typedef Howard<StaticDigraph, CostArcMap> MMC;
3.1449 + typedef typename BellmanFord<StaticDigraph, CostArcMap>
3.1450 + ::template SetDistMap<CostNodeMap>::Create BF;
3.1451 +
3.1452 + // Set epsilon to the minimum cycle mean
3.1453 + buildResidualNetwork();
3.1454 + MMC mmc(_sgr, _cost_map);
3.1455 + mmc.findMinMean();
3.1456 + epsilon = -mmc.cycleMean();
3.1457 + Cost cycle_cost = mmc.cycleLength();
3.1458 + int cycle_size = mmc.cycleArcNum();
3.1459 +
3.1460 + // Compute feasible potentials for the current epsilon
3.1461 + for (int i = 0; i != int(_cost_vec.size()); ++i) {
3.1462 + _cost_vec[i] = cycle_size * _cost_vec[i] - cycle_cost;
3.1463 + }
3.1464 + BF bf(_sgr, _cost_map);
3.1465 + bf.distMap(_pi_map);
3.1466 + bf.init(0);
3.1467 + bf.start();
3.1468 + for (int u = 0; u != _res_node_num; ++u) {
3.1469 + pi[u] = static_cast<double>(_pi[u]) / cycle_size;
3.1470 + }
3.1471 +
3.1472 + iter = limit;
3.1473 }
3.1474 }
3.1475 -
3.1476 - // Computing node potentials
3.1477 - BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost);
3.1478 - bf.init(0); bf.start();
3.1479 - for (NodeIt n(_graph); n != INVALID; ++n)
3.1480 - (*_potential)[n] = bf.dist(n);
3.1481 }
3.1482
3.1483 }; //class CycleCanceling