1.1 --- a/lemon/cancel_and_tighten.h Fri Nov 13 00:09:35 2009 +0100
1.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
1.3 @@ -1,797 +0,0 @@
1.4 -/* -*- C++ -*-
1.5 - *
1.6 - * This file is a part of LEMON, a generic C++ optimization library
1.7 - *
1.8 - * Copyright (C) 2003-2008
1.9 - * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 - * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 - *
1.12 - * Permission to use, modify and distribute this software is granted
1.13 - * provided that this copyright notice appears in all copies. For
1.14 - * precise terms see the accompanying LICENSE file.
1.15 - *
1.16 - * This software is provided "AS IS" with no warranty of any kind,
1.17 - * express or implied, and with no claim as to its suitability for any
1.18 - * purpose.
1.19 - *
1.20 - */
1.21 -
1.22 -#ifndef LEMON_CANCEL_AND_TIGHTEN_H
1.23 -#define LEMON_CANCEL_AND_TIGHTEN_H
1.24 -
1.25 -/// \ingroup min_cost_flow
1.26 -///
1.27 -/// \file
1.28 -/// \brief Cancel and Tighten algorithm for finding a minimum cost flow.
1.29 -
1.30 -#include <vector>
1.31 -
1.32 -#include <lemon/circulation.h>
1.33 -#include <lemon/bellman_ford.h>
1.34 -#include <lemon/howard.h>
1.35 -#include <lemon/adaptors.h>
1.36 -#include <lemon/tolerance.h>
1.37 -#include <lemon/math.h>
1.38 -
1.39 -#include <lemon/static_graph.h>
1.40 -
1.41 -namespace lemon {
1.42 -
1.43 - /// \addtogroup min_cost_flow
1.44 - /// @{
1.45 -
1.46 - /// \brief Implementation of the Cancel and Tighten algorithm for
1.47 - /// finding a minimum cost flow.
1.48 - ///
1.49 - /// \ref CancelAndTighten implements the Cancel and Tighten algorithm for
1.50 - /// finding a minimum cost flow.
1.51 - ///
1.52 - /// \tparam Digraph The digraph type the algorithm runs on.
1.53 - /// \tparam LowerMap The type of the lower bound map.
1.54 - /// \tparam CapacityMap The type of the capacity (upper bound) map.
1.55 - /// \tparam CostMap The type of the cost (length) map.
1.56 - /// \tparam SupplyMap The type of the supply map.
1.57 - ///
1.58 - /// \warning
1.59 - /// - Arc capacities and costs should be \e non-negative \e integers.
1.60 - /// - Supply values should be \e signed \e integers.
1.61 - /// - The value types of the maps should be convertible to each other.
1.62 - /// - \c CostMap::Value must be signed type.
1.63 - ///
1.64 - /// \author Peter Kovacs
1.65 - template < typename Digraph,
1.66 - typename LowerMap = typename Digraph::template ArcMap<int>,
1.67 - typename CapacityMap = typename Digraph::template ArcMap<int>,
1.68 - typename CostMap = typename Digraph::template ArcMap<int>,
1.69 - typename SupplyMap = typename Digraph::template NodeMap<int> >
1.70 - class CancelAndTighten
1.71 - {
1.72 - TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
1.73 -
1.74 - typedef typename CapacityMap::Value Capacity;
1.75 - typedef typename CostMap::Value Cost;
1.76 - typedef typename SupplyMap::Value Supply;
1.77 - typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
1.78 - typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
1.79 -
1.80 - typedef ResidualDigraph< const Digraph,
1.81 - CapacityArcMap, CapacityArcMap > ResDigraph;
1.82 -
1.83 - public:
1.84 -
1.85 - /// The type of the flow map.
1.86 - typedef typename Digraph::template ArcMap<Capacity> FlowMap;
1.87 - /// The type of the potential map.
1.88 - typedef typename Digraph::template NodeMap<Cost> PotentialMap;
1.89 -
1.90 - private:
1.91 -
1.92 - /// \brief Map adaptor class for handling residual arc costs.
1.93 - ///
1.94 - /// Map adaptor class for handling residual arc costs.
1.95 - class ResidualCostMap : public MapBase<typename ResDigraph::Arc, Cost>
1.96 - {
1.97 - typedef typename ResDigraph::Arc Arc;
1.98 -
1.99 - private:
1.100 -
1.101 - const CostMap &_cost_map;
1.102 -
1.103 - public:
1.104 -
1.105 - ///\e
1.106 - ResidualCostMap(const CostMap &cost_map) : _cost_map(cost_map) {}
1.107 -
1.108 - ///\e
1.109 - Cost operator[](const Arc &e) const {
1.110 - return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
1.111 - }
1.112 -
1.113 - }; //class ResidualCostMap
1.114 -
1.115 - /// \brief Map adaptor class for handling reduced arc costs.
1.116 - ///
1.117 - /// Map adaptor class for handling reduced arc costs.
1.118 - class ReducedCostMap : public MapBase<Arc, Cost>
1.119 - {
1.120 - private:
1.121 -
1.122 - const Digraph &_gr;
1.123 - const CostMap &_cost_map;
1.124 - const PotentialMap &_pot_map;
1.125 -
1.126 - public:
1.127 -
1.128 - ///\e
1.129 - ReducedCostMap( const Digraph &gr,
1.130 - const CostMap &cost_map,
1.131 - const PotentialMap &pot_map ) :
1.132 - _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
1.133 -
1.134 - ///\e
1.135 - inline Cost operator[](const Arc &e) const {
1.136 - return _cost_map[e] + _pot_map[_gr.source(e)]
1.137 - - _pot_map[_gr.target(e)];
1.138 - }
1.139 -
1.140 - }; //class ReducedCostMap
1.141 -
1.142 - struct BFOperationTraits {
1.143 - static double zero() { return 0; }
1.144 -
1.145 - static double infinity() {
1.146 - return std::numeric_limits<double>::infinity();
1.147 - }
1.148 -
1.149 - static double plus(const double& left, const double& right) {
1.150 - return left + right;
1.151 - }
1.152 -
1.153 - static bool less(const double& left, const double& right) {
1.154 - return left + 1e-6 < right;
1.155 - }
1.156 - }; // class BFOperationTraits
1.157 -
1.158 - private:
1.159 -
1.160 - // The digraph the algorithm runs on
1.161 - const Digraph &_graph;
1.162 - // The original lower bound map
1.163 - const LowerMap *_lower;
1.164 - // The modified capacity map
1.165 - CapacityArcMap _capacity;
1.166 - // The original cost map
1.167 - const CostMap &_cost;
1.168 - // The modified supply map
1.169 - SupplyNodeMap _supply;
1.170 - bool _valid_supply;
1.171 -
1.172 - // Arc map of the current flow
1.173 - FlowMap *_flow;
1.174 - bool _local_flow;
1.175 - // Node map of the current potentials
1.176 - PotentialMap *_potential;
1.177 - bool _local_potential;
1.178 -
1.179 - // The residual digraph
1.180 - ResDigraph *_res_graph;
1.181 - // The residual cost map
1.182 - ResidualCostMap _res_cost;
1.183 -
1.184 - public:
1.185 -
1.186 - /// \brief General constructor (with lower bounds).
1.187 - ///
1.188 - /// General constructor (with lower bounds).
1.189 - ///
1.190 - /// \param digraph The digraph the algorithm runs on.
1.191 - /// \param lower The lower bounds of the arcs.
1.192 - /// \param capacity The capacities (upper bounds) of the arcs.
1.193 - /// \param cost The cost (length) values of the arcs.
1.194 - /// \param supply The supply values of the nodes (signed).
1.195 - CancelAndTighten( const Digraph &digraph,
1.196 - const LowerMap &lower,
1.197 - const CapacityMap &capacity,
1.198 - const CostMap &cost,
1.199 - const SupplyMap &supply ) :
1.200 - _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost),
1.201 - _supply(digraph), _flow(NULL), _local_flow(false),
1.202 - _potential(NULL), _local_potential(false),
1.203 - _res_graph(NULL), _res_cost(_cost)
1.204 - {
1.205 - // Check the sum of supply values
1.206 - Supply sum = 0;
1.207 - for (NodeIt n(_graph); n != INVALID; ++n) {
1.208 - _supply[n] = supply[n];
1.209 - sum += _supply[n];
1.210 - }
1.211 - _valid_supply = sum == 0;
1.212 -
1.213 - // Remove non-zero lower bounds
1.214 - for (ArcIt e(_graph); e != INVALID; ++e) {
1.215 - _capacity[e] = capacity[e];
1.216 - if (lower[e] != 0) {
1.217 - _capacity[e] -= lower[e];
1.218 - _supply[_graph.source(e)] -= lower[e];
1.219 - _supply[_graph.target(e)] += lower[e];
1.220 - }
1.221 - }
1.222 - }
1.223 -/*
1.224 - /// \brief General constructor (without lower bounds).
1.225 - ///
1.226 - /// General constructor (without lower bounds).
1.227 - ///
1.228 - /// \param digraph The digraph the algorithm runs on.
1.229 - /// \param capacity The capacities (upper bounds) of the arcs.
1.230 - /// \param cost The cost (length) values of the arcs.
1.231 - /// \param supply The supply values of the nodes (signed).
1.232 - CancelAndTighten( const Digraph &digraph,
1.233 - const CapacityMap &capacity,
1.234 - const CostMap &cost,
1.235 - const SupplyMap &supply ) :
1.236 - _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
1.237 - _supply(supply), _flow(NULL), _local_flow(false),
1.238 - _potential(NULL), _local_potential(false),
1.239 - _res_graph(NULL), _res_cost(_cost)
1.240 - {
1.241 - // Check the sum of supply values
1.242 - Supply sum = 0;
1.243 - for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
1.244 - _valid_supply = sum == 0;
1.245 - }
1.246 -
1.247 - /// \brief Simple constructor (with lower bounds).
1.248 - ///
1.249 - /// Simple constructor (with lower bounds).
1.250 - ///
1.251 - /// \param digraph The digraph the algorithm runs on.
1.252 - /// \param lower The lower bounds of the arcs.
1.253 - /// \param capacity The capacities (upper bounds) of the arcs.
1.254 - /// \param cost The cost (length) values of the arcs.
1.255 - /// \param s The source node.
1.256 - /// \param t The target node.
1.257 - /// \param flow_value The required amount of flow from node \c s
1.258 - /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.259 - CancelAndTighten( const Digraph &digraph,
1.260 - const LowerMap &lower,
1.261 - const CapacityMap &capacity,
1.262 - const CostMap &cost,
1.263 - Node s, Node t,
1.264 - Supply flow_value ) :
1.265 - _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost),
1.266 - _supply(digraph, 0), _flow(NULL), _local_flow(false),
1.267 - _potential(NULL), _local_potential(false),
1.268 - _res_graph(NULL), _res_cost(_cost)
1.269 - {
1.270 - // Remove non-zero lower bounds
1.271 - _supply[s] = flow_value;
1.272 - _supply[t] = -flow_value;
1.273 - for (ArcIt e(_graph); e != INVALID; ++e) {
1.274 - if (lower[e] != 0) {
1.275 - _capacity[e] -= lower[e];
1.276 - _supply[_graph.source(e)] -= lower[e];
1.277 - _supply[_graph.target(e)] += lower[e];
1.278 - }
1.279 - }
1.280 - _valid_supply = true;
1.281 - }
1.282 -
1.283 - /// \brief Simple constructor (without lower bounds).
1.284 - ///
1.285 - /// Simple constructor (without lower bounds).
1.286 - ///
1.287 - /// \param digraph The digraph the algorithm runs on.
1.288 - /// \param capacity The capacities (upper bounds) of the arcs.
1.289 - /// \param cost The cost (length) values of the arcs.
1.290 - /// \param s The source node.
1.291 - /// \param t The target node.
1.292 - /// \param flow_value The required amount of flow from node \c s
1.293 - /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.294 - CancelAndTighten( const Digraph &digraph,
1.295 - const CapacityMap &capacity,
1.296 - const CostMap &cost,
1.297 - Node s, Node t,
1.298 - Supply flow_value ) :
1.299 - _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
1.300 - _supply(digraph, 0), _flow(NULL), _local_flow(false),
1.301 - _potential(NULL), _local_potential(false),
1.302 - _res_graph(NULL), _res_cost(_cost)
1.303 - {
1.304 - _supply[s] = flow_value;
1.305 - _supply[t] = -flow_value;
1.306 - _valid_supply = true;
1.307 - }
1.308 -*/
1.309 - /// Destructor.
1.310 - ~CancelAndTighten() {
1.311 - if (_local_flow) delete _flow;
1.312 - if (_local_potential) delete _potential;
1.313 - delete _res_graph;
1.314 - }
1.315 -
1.316 - /// \brief Set the flow map.
1.317 - ///
1.318 - /// Set the flow map.
1.319 - ///
1.320 - /// \return \c (*this)
1.321 - CancelAndTighten& flowMap(FlowMap &map) {
1.322 - if (_local_flow) {
1.323 - delete _flow;
1.324 - _local_flow = false;
1.325 - }
1.326 - _flow = ↦
1.327 - return *this;
1.328 - }
1.329 -
1.330 - /// \brief Set the potential map.
1.331 - ///
1.332 - /// Set the potential map.
1.333 - ///
1.334 - /// \return \c (*this)
1.335 - CancelAndTighten& potentialMap(PotentialMap &map) {
1.336 - if (_local_potential) {
1.337 - delete _potential;
1.338 - _local_potential = false;
1.339 - }
1.340 - _potential = ↦
1.341 - return *this;
1.342 - }
1.343 -
1.344 - /// \name Execution control
1.345 -
1.346 - /// @{
1.347 -
1.348 - /// \brief Run the algorithm.
1.349 - ///
1.350 - /// Run the algorithm.
1.351 - ///
1.352 - /// \return \c true if a feasible flow can be found.
1.353 - bool run() {
1.354 - return init() && start();
1.355 - }
1.356 -
1.357 - /// @}
1.358 -
1.359 - /// \name Query Functions
1.360 - /// The result of the algorithm can be obtained using these
1.361 - /// functions.\n
1.362 - /// \ref lemon::CancelAndTighten::run() "run()" must be called before
1.363 - /// using them.
1.364 -
1.365 - /// @{
1.366 -
1.367 - /// \brief Return a const reference to the arc map storing the
1.368 - /// found flow.
1.369 - ///
1.370 - /// Return a const reference to the arc map storing the found flow.
1.371 - ///
1.372 - /// \pre \ref run() must be called before using this function.
1.373 - const FlowMap& flowMap() const {
1.374 - return *_flow;
1.375 - }
1.376 -
1.377 - /// \brief Return a const reference to the node map storing the
1.378 - /// found potentials (the dual solution).
1.379 - ///
1.380 - /// Return a const reference to the node map storing the found
1.381 - /// potentials (the dual solution).
1.382 - ///
1.383 - /// \pre \ref run() must be called before using this function.
1.384 - const PotentialMap& potentialMap() const {
1.385 - return *_potential;
1.386 - }
1.387 -
1.388 - /// \brief Return the flow on the given arc.
1.389 - ///
1.390 - /// Return the flow on the given arc.
1.391 - ///
1.392 - /// \pre \ref run() must be called before using this function.
1.393 - Capacity flow(const Arc& arc) const {
1.394 - return (*_flow)[arc];
1.395 - }
1.396 -
1.397 - /// \brief Return the potential of the given node.
1.398 - ///
1.399 - /// Return the potential of the given node.
1.400 - ///
1.401 - /// \pre \ref run() must be called before using this function.
1.402 - Cost potential(const Node& node) const {
1.403 - return (*_potential)[node];
1.404 - }
1.405 -
1.406 - /// \brief Return the total cost of the found flow.
1.407 - ///
1.408 - /// Return the total cost of the found flow. The complexity of the
1.409 - /// function is \f$ O(e) \f$.
1.410 - ///
1.411 - /// \pre \ref run() must be called before using this function.
1.412 - Cost totalCost() const {
1.413 - Cost c = 0;
1.414 - for (ArcIt e(_graph); e != INVALID; ++e)
1.415 - c += (*_flow)[e] * _cost[e];
1.416 - return c;
1.417 - }
1.418 -
1.419 - /// @}
1.420 -
1.421 - private:
1.422 -
1.423 - /// Initialize the algorithm.
1.424 - bool init() {
1.425 - if (!_valid_supply) return false;
1.426 -
1.427 - // Initialize flow and potential maps
1.428 - if (!_flow) {
1.429 - _flow = new FlowMap(_graph);
1.430 - _local_flow = true;
1.431 - }
1.432 - if (!_potential) {
1.433 - _potential = new PotentialMap(_graph);
1.434 - _local_potential = true;
1.435 - }
1.436 -
1.437 - _res_graph = new ResDigraph(_graph, _capacity, *_flow);
1.438 -
1.439 - // Find a feasible flow using Circulation
1.440 - Circulation< Digraph, ConstMap<Arc, Capacity>,
1.441 - CapacityArcMap, SupplyMap >
1.442 - circulation( _graph, constMap<Arc>(Capacity(0)),
1.443 - _capacity, _supply );
1.444 - return circulation.flowMap(*_flow).run();
1.445 - }
1.446 -
1.447 - bool start() {
1.448 - const double LIMIT_FACTOR = 0.01;
1.449 - const int MIN_LIMIT = 3;
1.450 -
1.451 - typedef typename Digraph::template NodeMap<double> FloatPotentialMap;
1.452 - typedef typename Digraph::template NodeMap<int> LevelMap;
1.453 - typedef typename Digraph::template NodeMap<bool> BoolNodeMap;
1.454 - typedef typename Digraph::template NodeMap<Node> PredNodeMap;
1.455 - typedef typename Digraph::template NodeMap<Arc> PredArcMap;
1.456 - typedef typename ResDigraph::template ArcMap<double> ResShiftCostMap;
1.457 - FloatPotentialMap pi(_graph);
1.458 - LevelMap level(_graph);
1.459 - BoolNodeMap reached(_graph);
1.460 - BoolNodeMap processed(_graph);
1.461 - PredNodeMap pred_node(_graph);
1.462 - PredArcMap pred_arc(_graph);
1.463 - int node_num = countNodes(_graph);
1.464 - typedef std::pair<Arc, bool> pair;
1.465 - std::vector<pair> stack(node_num);
1.466 - std::vector<Node> proc_vector(node_num);
1.467 - ResShiftCostMap shift_cost(*_res_graph);
1.468 -
1.469 - Tolerance<double> tol;
1.470 - tol.epsilon(1e-6);
1.471 -
1.472 - Timer t1, t2, t3;
1.473 - t1.reset();
1.474 - t2.reset();
1.475 - t3.reset();
1.476 -
1.477 - // Initialize epsilon and the node potentials
1.478 - double epsilon = 0;
1.479 - for (ArcIt e(_graph); e != INVALID; ++e) {
1.480 - if (_capacity[e] - (*_flow)[e] > 0 && _cost[e] < -epsilon)
1.481 - epsilon = -_cost[e];
1.482 - else if ((*_flow)[e] > 0 && _cost[e] > epsilon)
1.483 - epsilon = _cost[e];
1.484 - }
1.485 - for (NodeIt v(_graph); v != INVALID; ++v) {
1.486 - pi[v] = 0;
1.487 - }
1.488 -
1.489 - // Start phases
1.490 - int limit = int(LIMIT_FACTOR * node_num);
1.491 - if (limit < MIN_LIMIT) limit = MIN_LIMIT;
1.492 - int iter = limit;
1.493 - while (epsilon * node_num >= 1) {
1.494 - t1.start();
1.495 - // Find and cancel cycles in the admissible digraph using DFS
1.496 - for (NodeIt n(_graph); n != INVALID; ++n) {
1.497 - reached[n] = false;
1.498 - processed[n] = false;
1.499 - }
1.500 - int stack_head = -1;
1.501 - int proc_head = -1;
1.502 -
1.503 - for (NodeIt start(_graph); start != INVALID; ++start) {
1.504 - if (reached[start]) continue;
1.505 -
1.506 - // New start node
1.507 - reached[start] = true;
1.508 - pred_arc[start] = INVALID;
1.509 - pred_node[start] = INVALID;
1.510 -
1.511 - // Find the first admissible residual outgoing arc
1.512 - double p = pi[start];
1.513 - Arc e;
1.514 - _graph.firstOut(e, start);
1.515 - while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
1.516 - !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
1.517 - _graph.nextOut(e);
1.518 - if (e != INVALID) {
1.519 - stack[++stack_head] = pair(e, true);
1.520 - goto next_step_1;
1.521 - }
1.522 - _graph.firstIn(e, start);
1.523 - while ( e != INVALID && ((*_flow)[e] == 0 ||
1.524 - !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
1.525 - _graph.nextIn(e);
1.526 - if (e != INVALID) {
1.527 - stack[++stack_head] = pair(e, false);
1.528 - goto next_step_1;
1.529 - }
1.530 - processed[start] = true;
1.531 - proc_vector[++proc_head] = start;
1.532 - continue;
1.533 - next_step_1:
1.534 -
1.535 - while (stack_head >= 0) {
1.536 - Arc se = stack[stack_head].first;
1.537 - bool sf = stack[stack_head].second;
1.538 - Node u, v;
1.539 - if (sf) {
1.540 - u = _graph.source(se);
1.541 - v = _graph.target(se);
1.542 - } else {
1.543 - u = _graph.target(se);
1.544 - v = _graph.source(se);
1.545 - }
1.546 -
1.547 - if (!reached[v]) {
1.548 - // A new node is reached
1.549 - reached[v] = true;
1.550 - pred_node[v] = u;
1.551 - pred_arc[v] = se;
1.552 - // Find the first admissible residual outgoing arc
1.553 - double p = pi[v];
1.554 - Arc e;
1.555 - _graph.firstOut(e, v);
1.556 - while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
1.557 - !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
1.558 - _graph.nextOut(e);
1.559 - if (e != INVALID) {
1.560 - stack[++stack_head] = pair(e, true);
1.561 - goto next_step_2;
1.562 - }
1.563 - _graph.firstIn(e, v);
1.564 - while ( e != INVALID && ((*_flow)[e] == 0 ||
1.565 - !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
1.566 - _graph.nextIn(e);
1.567 - stack[++stack_head] = pair(e, false);
1.568 - next_step_2: ;
1.569 - } else {
1.570 - if (!processed[v]) {
1.571 - // A cycle is found
1.572 - Node n, w = u;
1.573 - Capacity d, delta = sf ? _capacity[se] - (*_flow)[se] :
1.574 - (*_flow)[se];
1.575 - for (n = u; n != v; n = pred_node[n]) {
1.576 - d = _graph.target(pred_arc[n]) == n ?
1.577 - _capacity[pred_arc[n]] - (*_flow)[pred_arc[n]] :
1.578 - (*_flow)[pred_arc[n]];
1.579 - if (d <= delta) {
1.580 - delta = d;
1.581 - w = pred_node[n];
1.582 - }
1.583 - }
1.584 -
1.585 -/*
1.586 - std::cout << "CYCLE FOUND: ";
1.587 - if (sf)
1.588 - std::cout << _cost[se] + pi[_graph.source(se)] - pi[_graph.target(se)];
1.589 - else
1.590 - std::cout << _graph.id(se) << ":" << -(_cost[se] + pi[_graph.source(se)] - pi[_graph.target(se)]);
1.591 - for (n = u; n != v; n = pred_node[n]) {
1.592 - if (_graph.target(pred_arc[n]) == n)
1.593 - std::cout << " " << _cost[pred_arc[n]] + pi[_graph.source(pred_arc[n])] - pi[_graph.target(pred_arc[n])];
1.594 - else
1.595 - std::cout << " " << -(_cost[pred_arc[n]] + pi[_graph.source(pred_arc[n])] - pi[_graph.target(pred_arc[n])]);
1.596 - }
1.597 - std::cout << "\n";
1.598 -*/
1.599 - // Augment along the cycle
1.600 - (*_flow)[se] = sf ? (*_flow)[se] + delta :
1.601 - (*_flow)[se] - delta;
1.602 - for (n = u; n != v; n = pred_node[n]) {
1.603 - if (_graph.target(pred_arc[n]) == n)
1.604 - (*_flow)[pred_arc[n]] += delta;
1.605 - else
1.606 - (*_flow)[pred_arc[n]] -= delta;
1.607 - }
1.608 - for (n = u; stack_head > 0 && n != w; n = pred_node[n]) {
1.609 - --stack_head;
1.610 - reached[n] = false;
1.611 - }
1.612 - u = w;
1.613 - }
1.614 - v = u;
1.615 -
1.616 - // Find the next admissible residual outgoing arc
1.617 - double p = pi[v];
1.618 - Arc e = stack[stack_head].first;
1.619 - if (!stack[stack_head].second) {
1.620 - _graph.nextIn(e);
1.621 - goto in_arc_3;
1.622 - }
1.623 - _graph.nextOut(e);
1.624 - while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
1.625 - !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
1.626 - _graph.nextOut(e);
1.627 - if (e != INVALID) {
1.628 - stack[stack_head] = pair(e, true);
1.629 - goto next_step_3;
1.630 - }
1.631 - _graph.firstIn(e, v);
1.632 - in_arc_3:
1.633 - while ( e != INVALID && ((*_flow)[e] == 0 ||
1.634 - !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
1.635 - _graph.nextIn(e);
1.636 - stack[stack_head] = pair(e, false);
1.637 - next_step_3: ;
1.638 - }
1.639 -
1.640 - while (stack_head >= 0 && stack[stack_head].first == INVALID) {
1.641 - processed[v] = true;
1.642 - proc_vector[++proc_head] = v;
1.643 - if (--stack_head >= 0) {
1.644 - v = stack[stack_head].second ?
1.645 - _graph.source(stack[stack_head].first) :
1.646 - _graph.target(stack[stack_head].first);
1.647 - // Find the next admissible residual outgoing arc
1.648 - double p = pi[v];
1.649 - Arc e = stack[stack_head].first;
1.650 - if (!stack[stack_head].second) {
1.651 - _graph.nextIn(e);
1.652 - goto in_arc_4;
1.653 - }
1.654 - _graph.nextOut(e);
1.655 - while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
1.656 - !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
1.657 - _graph.nextOut(e);
1.658 - if (e != INVALID) {
1.659 - stack[stack_head] = pair(e, true);
1.660 - goto next_step_4;
1.661 - }
1.662 - _graph.firstIn(e, v);
1.663 - in_arc_4:
1.664 - while ( e != INVALID && ((*_flow)[e] == 0 ||
1.665 - !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
1.666 - _graph.nextIn(e);
1.667 - stack[stack_head] = pair(e, false);
1.668 - next_step_4: ;
1.669 - }
1.670 - }
1.671 - }
1.672 - }
1.673 - t1.stop();
1.674 -
1.675 - // Tighten potentials and epsilon
1.676 - if (--iter > 0) {
1.677 - // Compute levels
1.678 - t2.start();
1.679 - for (int i = proc_head; i >= 0; --i) {
1.680 - Node v = proc_vector[i];
1.681 - double p = pi[v];
1.682 - int l = 0;
1.683 - for (InArcIt e(_graph, v); e != INVALID; ++e) {
1.684 - Node u = _graph.source(e);
1.685 - if ( _capacity[e] - (*_flow)[e] > 0 &&
1.686 - tol.negative(_cost[e] + pi[u] - p) &&
1.687 - level[u] + 1 > l ) l = level[u] + 1;
1.688 - }
1.689 - for (OutArcIt e(_graph, v); e != INVALID; ++e) {
1.690 - Node u = _graph.target(e);
1.691 - if ( (*_flow)[e] > 0 &&
1.692 - tol.negative(-_cost[e] + pi[u] - p) &&
1.693 - level[u] + 1 > l ) l = level[u] + 1;
1.694 - }
1.695 - level[v] = l;
1.696 - }
1.697 -
1.698 - // Modify potentials
1.699 - double p, q = -1;
1.700 - for (ArcIt e(_graph); e != INVALID; ++e) {
1.701 - Node u = _graph.source(e);
1.702 - Node v = _graph.target(e);
1.703 - if (_capacity[e] - (*_flow)[e] > 0 && level[u] - level[v] > 0) {
1.704 - p = (_cost[e] + pi[u] - pi[v] + epsilon) /
1.705 - (level[u] - level[v] + 1);
1.706 - if (q < 0 || p < q) q = p;
1.707 - }
1.708 - else if ((*_flow)[e] > 0 && level[v] - level[u] > 0) {
1.709 - p = (-_cost[e] - pi[u] + pi[v] + epsilon) /
1.710 - (level[v] - level[u] + 1);
1.711 - if (q < 0 || p < q) q = p;
1.712 - }
1.713 - }
1.714 - for (NodeIt v(_graph); v != INVALID; ++v) {
1.715 - pi[v] -= q * level[v];
1.716 - }
1.717 -
1.718 - // Modify epsilon
1.719 - epsilon = 0;
1.720 - for (ArcIt e(_graph); e != INVALID; ++e) {
1.721 - double curr = _cost[e] + pi[_graph.source(e)]
1.722 - - pi[_graph.target(e)];
1.723 - if (_capacity[e] - (*_flow)[e] > 0 && curr < -epsilon)
1.724 - epsilon = -curr;
1.725 - else if ((*_flow)[e] > 0 && curr > epsilon)
1.726 - epsilon = curr;
1.727 - }
1.728 - t2.stop();
1.729 - } else {
1.730 - // Set epsilon to the minimum cycle mean
1.731 - t3.start();
1.732 -
1.733 -/**/
1.734 - StaticDigraph static_graph;
1.735 - typename ResDigraph::template NodeMap<typename StaticDigraph::Node> node_ref(*_res_graph);
1.736 - typename ResDigraph::template ArcMap<typename StaticDigraph::Arc> arc_ref(*_res_graph);
1.737 - static_graph.build(*_res_graph, node_ref, arc_ref);
1.738 - typename StaticDigraph::template NodeMap<double> static_pi(static_graph);
1.739 - typename StaticDigraph::template ArcMap<double> static_cost(static_graph);
1.740 -
1.741 - for (typename ResDigraph::ArcIt e(*_res_graph); e != INVALID; ++e)
1.742 - static_cost[arc_ref[e]] = _res_cost[e];
1.743 -
1.744 - Howard<StaticDigraph, typename StaticDigraph::template ArcMap<double> >
1.745 - mmc(static_graph, static_cost);
1.746 - mmc.findMinMean();
1.747 - epsilon = -mmc.cycleMean();
1.748 -/**/
1.749 -
1.750 -/*
1.751 - Howard<ResDigraph, ResidualCostMap> mmc(*_res_graph, _res_cost);
1.752 - mmc.findMinMean();
1.753 - epsilon = -mmc.cycleMean();
1.754 -*/
1.755 -
1.756 - // Compute feasible potentials for the current epsilon
1.757 - for (typename StaticDigraph::ArcIt e(static_graph); e != INVALID; ++e)
1.758 - static_cost[e] += epsilon;
1.759 - typename BellmanFord<StaticDigraph, typename StaticDigraph::template ArcMap<double> >::
1.760 - template SetDistMap<typename StaticDigraph::template NodeMap<double> >::
1.761 - template SetOperationTraits<BFOperationTraits>::Create
1.762 - bf(static_graph, static_cost);
1.763 - bf.distMap(static_pi).init(0);
1.764 - bf.start();
1.765 - for (NodeIt n(_graph); n != INVALID; ++n)
1.766 - pi[n] = static_pi[node_ref[n]];
1.767 -
1.768 -/*
1.769 - for (typename ResDigraph::ArcIt e(*_res_graph); e != INVALID; ++e)
1.770 - shift_cost[e] = _res_cost[e] + epsilon;
1.771 - typename BellmanFord<ResDigraph, ResShiftCostMap>::
1.772 - template SetDistMap<FloatPotentialMap>::
1.773 - template SetOperationTraits<BFOperationTraits>::Create
1.774 - bf(*_res_graph, shift_cost);
1.775 - bf.distMap(pi).init(0);
1.776 - bf.start();
1.777 -*/
1.778 -
1.779 - iter = limit;
1.780 - t3.stop();
1.781 - }
1.782 - }
1.783 -
1.784 -// std::cout << t1.realTime() << " " << t2.realTime() << " " << t3.realTime() << "\n";
1.785 -
1.786 - // Handle non-zero lower bounds
1.787 - if (_lower) {
1.788 - for (ArcIt e(_graph); e != INVALID; ++e)
1.789 - (*_flow)[e] += (*_lower)[e];
1.790 - }
1.791 - return true;
1.792 - }
1.793 -
1.794 - }; //class CancelAndTighten
1.795 -
1.796 - ///@}
1.797 -
1.798 -} //namespace lemon
1.799 -
1.800 -#endif //LEMON_CANCEL_AND_TIGHTEN_H