1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/cancel_and_tighten.h Mon Jun 01 15:37:51 2009 +0000
1.3 @@ -0,0 +1,795 @@
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/min_mean_cycle.h>
1.35 +#include <lemon/graph_adaptor.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 Graph The directed graph 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 + /// - Edge 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 Graph,
1.66 + typename LowerMap = typename Graph::template EdgeMap<int>,
1.67 + typename CapacityMap = typename Graph::template EdgeMap<int>,
1.68 + typename CostMap = typename Graph::template EdgeMap<int>,
1.69 + typename SupplyMap = typename Graph::template NodeMap<int> >
1.70 + class CancelAndTighten
1.71 + {
1.72 + GRAPH_TYPEDEFS(typename Graph);
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 Graph::template EdgeMap<Capacity> CapacityEdgeMap;
1.78 + typedef typename Graph::template NodeMap<Supply> SupplyNodeMap;
1.79 +
1.80 + typedef ResGraphAdaptor< const Graph, Capacity,
1.81 + CapacityEdgeMap, CapacityEdgeMap > ResGraph;
1.82 +
1.83 + public:
1.84 +
1.85 + /// The type of the flow map.
1.86 + typedef typename Graph::template EdgeMap<Capacity> FlowMap;
1.87 + /// The type of the potential map.
1.88 + typedef typename Graph::template NodeMap<Cost> PotentialMap;
1.89 +
1.90 + private:
1.91 +
1.92 + /// \brief Map adaptor class for handling residual edge costs.
1.93 + ///
1.94 + /// Map adaptor class for handling residual edge costs.
1.95 + class ResidualCostMap : public MapBase<typename ResGraph::Edge, Cost>
1.96 + {
1.97 + typedef typename ResGraph::Edge Edge;
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 Edge &e) const {
1.110 + return ResGraph::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 edge costs.
1.116 + ///
1.117 + /// Map adaptor class for handling reduced edge costs.
1.118 + class ReducedCostMap : public MapBase<Edge, Cost>
1.119 + {
1.120 + private:
1.121 +
1.122 + const Graph &_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 Graph &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 Edge &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 directed graph the algorithm runs on
1.161 + const Graph &_graph;
1.162 + // The original lower bound map
1.163 + const LowerMap *_lower;
1.164 + // The modified capacity map
1.165 + CapacityEdgeMap _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 + // Edge 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 graph
1.180 + ResGraph *_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 graph The directed graph the algorithm runs on.
1.191 + /// \param lower The lower bounds of the edges.
1.192 + /// \param capacity The capacities (upper bounds) of the edges.
1.193 + /// \param cost The cost (length) values of the edges.
1.194 + /// \param supply The supply values of the nodes (signed).
1.195 + CancelAndTighten( const Graph &graph,
1.196 + const LowerMap &lower,
1.197 + const CapacityMap &capacity,
1.198 + const CostMap &cost,
1.199 + const SupplyMap &supply ) :
1.200 + _graph(graph), _lower(&lower), _capacity(capacity), _cost(cost),
1.201 + _supply(supply), _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) sum += _supply[n];
1.208 + _valid_supply = sum == 0;
1.209 +
1.210 + // Remove non-zero lower bounds
1.211 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.212 + if (lower[e] != 0) {
1.213 + _capacity[e] -= lower[e];
1.214 + _supply[_graph.source(e)] -= lower[e];
1.215 + _supply[_graph.target(e)] += lower[e];
1.216 + }
1.217 + }
1.218 + }
1.219 +
1.220 + /// \brief General constructor (without lower bounds).
1.221 + ///
1.222 + /// General constructor (without lower bounds).
1.223 + ///
1.224 + /// \param graph The directed graph the algorithm runs on.
1.225 + /// \param capacity The capacities (upper bounds) of the edges.
1.226 + /// \param cost The cost (length) values of the edges.
1.227 + /// \param supply The supply values of the nodes (signed).
1.228 + CancelAndTighten( const Graph &graph,
1.229 + const CapacityMap &capacity,
1.230 + const CostMap &cost,
1.231 + const SupplyMap &supply ) :
1.232 + _graph(graph), _lower(NULL), _capacity(capacity), _cost(cost),
1.233 + _supply(supply), _flow(NULL), _local_flow(false),
1.234 + _potential(NULL), _local_potential(false),
1.235 + _res_graph(NULL), _res_cost(_cost)
1.236 + {
1.237 + // Check the sum of supply values
1.238 + Supply sum = 0;
1.239 + for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
1.240 + _valid_supply = sum == 0;
1.241 + }
1.242 +
1.243 + /// \brief Simple constructor (with lower bounds).
1.244 + ///
1.245 + /// Simple constructor (with lower bounds).
1.246 + ///
1.247 + /// \param graph The directed graph the algorithm runs on.
1.248 + /// \param lower The lower bounds of the edges.
1.249 + /// \param capacity The capacities (upper bounds) of the edges.
1.250 + /// \param cost The cost (length) values of the edges.
1.251 + /// \param s The source node.
1.252 + /// \param t The target node.
1.253 + /// \param flow_value The required amount of flow from node \c s
1.254 + /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.255 + CancelAndTighten( const Graph &graph,
1.256 + const LowerMap &lower,
1.257 + const CapacityMap &capacity,
1.258 + const CostMap &cost,
1.259 + Node s, Node t,
1.260 + Supply flow_value ) :
1.261 + _graph(graph), _lower(&lower), _capacity(capacity), _cost(cost),
1.262 + _supply(graph, 0), _flow(NULL), _local_flow(false),
1.263 + _potential(NULL), _local_potential(false),
1.264 + _res_graph(NULL), _res_cost(_cost)
1.265 + {
1.266 + // Remove non-zero lower bounds
1.267 + _supply[s] = flow_value;
1.268 + _supply[t] = -flow_value;
1.269 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.270 + if (lower[e] != 0) {
1.271 + _capacity[e] -= lower[e];
1.272 + _supply[_graph.source(e)] -= lower[e];
1.273 + _supply[_graph.target(e)] += lower[e];
1.274 + }
1.275 + }
1.276 + _valid_supply = true;
1.277 + }
1.278 +
1.279 + /// \brief Simple constructor (without lower bounds).
1.280 + ///
1.281 + /// Simple constructor (without lower bounds).
1.282 + ///
1.283 + /// \param graph The directed graph the algorithm runs on.
1.284 + /// \param capacity The capacities (upper bounds) of the edges.
1.285 + /// \param cost The cost (length) values of the edges.
1.286 + /// \param s The source node.
1.287 + /// \param t The target node.
1.288 + /// \param flow_value The required amount of flow from node \c s
1.289 + /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.290 + CancelAndTighten( const Graph &graph,
1.291 + const CapacityMap &capacity,
1.292 + const CostMap &cost,
1.293 + Node s, Node t,
1.294 + Supply flow_value ) :
1.295 + _graph(graph), _lower(NULL), _capacity(capacity), _cost(cost),
1.296 + _supply(graph, 0), _flow(NULL), _local_flow(false),
1.297 + _potential(NULL), _local_potential(false),
1.298 + _res_graph(NULL), _res_cost(_cost)
1.299 + {
1.300 + _supply[s] = flow_value;
1.301 + _supply[t] = -flow_value;
1.302 + _valid_supply = true;
1.303 + }
1.304 +
1.305 + /// Destructor.
1.306 + ~CancelAndTighten() {
1.307 + if (_local_flow) delete _flow;
1.308 + if (_local_potential) delete _potential;
1.309 + delete _res_graph;
1.310 + }
1.311 +
1.312 + /// \brief Set the flow map.
1.313 + ///
1.314 + /// Set the flow map.
1.315 + ///
1.316 + /// \return \c (*this)
1.317 + CancelAndTighten& flowMap(FlowMap &map) {
1.318 + if (_local_flow) {
1.319 + delete _flow;
1.320 + _local_flow = false;
1.321 + }
1.322 + _flow = ↦
1.323 + return *this;
1.324 + }
1.325 +
1.326 + /// \brief Set the potential map.
1.327 + ///
1.328 + /// Set the potential map.
1.329 + ///
1.330 + /// \return \c (*this)
1.331 + CancelAndTighten& potentialMap(PotentialMap &map) {
1.332 + if (_local_potential) {
1.333 + delete _potential;
1.334 + _local_potential = false;
1.335 + }
1.336 + _potential = ↦
1.337 + return *this;
1.338 + }
1.339 +
1.340 + /// \name Execution control
1.341 +
1.342 + /// @{
1.343 +
1.344 + /// \brief Run the algorithm.
1.345 + ///
1.346 + /// Run the algorithm.
1.347 + ///
1.348 + /// \return \c true if a feasible flow can be found.
1.349 + bool run() {
1.350 + return init() && start();
1.351 + }
1.352 +
1.353 + /// @}
1.354 +
1.355 + /// \name Query Functions
1.356 + /// The result of the algorithm can be obtained using these
1.357 + /// functions.\n
1.358 + /// \ref lemon::CancelAndTighten::run() "run()" must be called before
1.359 + /// using them.
1.360 +
1.361 + /// @{
1.362 +
1.363 + /// \brief Return a const reference to the edge map storing the
1.364 + /// found flow.
1.365 + ///
1.366 + /// Return a const reference to the edge map storing the found flow.
1.367 + ///
1.368 + /// \pre \ref run() must be called before using this function.
1.369 + const FlowMap& flowMap() const {
1.370 + return *_flow;
1.371 + }
1.372 +
1.373 + /// \brief Return a const reference to the node map storing the
1.374 + /// found potentials (the dual solution).
1.375 + ///
1.376 + /// Return a const reference to the node map storing the found
1.377 + /// potentials (the dual solution).
1.378 + ///
1.379 + /// \pre \ref run() must be called before using this function.
1.380 + const PotentialMap& potentialMap() const {
1.381 + return *_potential;
1.382 + }
1.383 +
1.384 + /// \brief Return the flow on the given edge.
1.385 + ///
1.386 + /// Return the flow on the given edge.
1.387 + ///
1.388 + /// \pre \ref run() must be called before using this function.
1.389 + Capacity flow(const Edge& edge) const {
1.390 + return (*_flow)[edge];
1.391 + }
1.392 +
1.393 + /// \brief Return the potential of the given node.
1.394 + ///
1.395 + /// Return the potential of the given node.
1.396 + ///
1.397 + /// \pre \ref run() must be called before using this function.
1.398 + Cost potential(const Node& node) const {
1.399 + return (*_potential)[node];
1.400 + }
1.401 +
1.402 + /// \brief Return the total cost of the found flow.
1.403 + ///
1.404 + /// Return the total cost of the found flow. The complexity of the
1.405 + /// function is \f$ O(e) \f$.
1.406 + ///
1.407 + /// \pre \ref run() must be called before using this function.
1.408 + Cost totalCost() const {
1.409 + Cost c = 0;
1.410 + for (EdgeIt e(_graph); e != INVALID; ++e)
1.411 + c += (*_flow)[e] * _cost[e];
1.412 + return c;
1.413 + }
1.414 +
1.415 + /// @}
1.416 +
1.417 + private:
1.418 +
1.419 + /// Initialize the algorithm.
1.420 + bool init() {
1.421 + if (!_valid_supply) return false;
1.422 +
1.423 + // Initialize flow and potential maps
1.424 + if (!_flow) {
1.425 + _flow = new FlowMap(_graph);
1.426 + _local_flow = true;
1.427 + }
1.428 + if (!_potential) {
1.429 + _potential = new PotentialMap(_graph);
1.430 + _local_potential = true;
1.431 + }
1.432 +
1.433 + _res_graph = new ResGraph(_graph, _capacity, *_flow);
1.434 +
1.435 + // Find a feasible flow using Circulation
1.436 + Circulation< Graph, ConstMap<Edge, Capacity>,
1.437 + CapacityEdgeMap, SupplyMap >
1.438 + circulation( _graph, constMap<Edge>(Capacity(0)),
1.439 + _capacity, _supply );
1.440 + return circulation.flowMap(*_flow).run();
1.441 + }
1.442 +
1.443 + bool start() {
1.444 + const double LIMIT_FACTOR = 0.01;
1.445 + const int MIN_LIMIT = 3;
1.446 +
1.447 + typedef typename Graph::template NodeMap<double> FloatPotentialMap;
1.448 + typedef typename Graph::template NodeMap<int> LevelMap;
1.449 + typedef typename Graph::template NodeMap<bool> BoolNodeMap;
1.450 + typedef typename Graph::template NodeMap<Node> PredNodeMap;
1.451 + typedef typename Graph::template NodeMap<Edge> PredEdgeMap;
1.452 + typedef typename ResGraph::template EdgeMap<double> ResShiftCostMap;
1.453 + FloatPotentialMap pi(_graph);
1.454 + LevelMap level(_graph);
1.455 + BoolNodeMap reached(_graph);
1.456 + BoolNodeMap processed(_graph);
1.457 + PredNodeMap pred_node(_graph);
1.458 + PredEdgeMap pred_edge(_graph);
1.459 + int node_num = countNodes(_graph);
1.460 + typedef std::pair<Edge, bool> pair;
1.461 + std::vector<pair> stack(node_num);
1.462 + std::vector<Node> proc_vector(node_num);
1.463 + ResShiftCostMap shift_cost(*_res_graph);
1.464 +
1.465 + Tolerance<double> tol;
1.466 + tol.epsilon(1e-6);
1.467 +
1.468 + Timer t1, t2, t3;
1.469 + t1.reset();
1.470 + t2.reset();
1.471 + t3.reset();
1.472 +
1.473 + // Initialize epsilon and the node potentials
1.474 + double epsilon = 0;
1.475 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.476 + if (_capacity[e] - (*_flow)[e] > 0 && _cost[e] < -epsilon)
1.477 + epsilon = -_cost[e];
1.478 + else if ((*_flow)[e] > 0 && _cost[e] > epsilon)
1.479 + epsilon = _cost[e];
1.480 + }
1.481 + for (NodeIt v(_graph); v != INVALID; ++v) {
1.482 + pi[v] = 0;
1.483 + }
1.484 +
1.485 + // Start phases
1.486 + int limit = int(LIMIT_FACTOR * node_num);
1.487 + if (limit < MIN_LIMIT) limit = MIN_LIMIT;
1.488 + int iter = limit;
1.489 + while (epsilon * node_num >= 1) {
1.490 + t1.start();
1.491 + // Find and cancel cycles in the admissible graph using DFS
1.492 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.493 + reached[n] = false;
1.494 + processed[n] = false;
1.495 + }
1.496 + int stack_head = -1;
1.497 + int proc_head = -1;
1.498 +
1.499 + for (NodeIt start(_graph); start != INVALID; ++start) {
1.500 + if (reached[start]) continue;
1.501 +
1.502 + // New start node
1.503 + reached[start] = true;
1.504 + pred_edge[start] = INVALID;
1.505 + pred_node[start] = INVALID;
1.506 +
1.507 + // Find the first admissible residual outgoing edge
1.508 + double p = pi[start];
1.509 + Edge e;
1.510 + _graph.firstOut(e, start);
1.511 + while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
1.512 + !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
1.513 + _graph.nextOut(e);
1.514 + if (e != INVALID) {
1.515 + stack[++stack_head] = pair(e, true);
1.516 + goto next_step_1;
1.517 + }
1.518 + _graph.firstIn(e, start);
1.519 + while ( e != INVALID && ((*_flow)[e] == 0 ||
1.520 + !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
1.521 + _graph.nextIn(e);
1.522 + if (e != INVALID) {
1.523 + stack[++stack_head] = pair(e, false);
1.524 + goto next_step_1;
1.525 + }
1.526 + processed[start] = true;
1.527 + proc_vector[++proc_head] = start;
1.528 + continue;
1.529 + next_step_1:
1.530 +
1.531 + while (stack_head >= 0) {
1.532 + Edge se = stack[stack_head].first;
1.533 + bool sf = stack[stack_head].second;
1.534 + Node u, v;
1.535 + if (sf) {
1.536 + u = _graph.source(se);
1.537 + v = _graph.target(se);
1.538 + } else {
1.539 + u = _graph.target(se);
1.540 + v = _graph.source(se);
1.541 + }
1.542 +
1.543 + if (!reached[v]) {
1.544 + // A new node is reached
1.545 + reached[v] = true;
1.546 + pred_node[v] = u;
1.547 + pred_edge[v] = se;
1.548 + // Find the first admissible residual outgoing edge
1.549 + double p = pi[v];
1.550 + Edge e;
1.551 + _graph.firstOut(e, v);
1.552 + while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
1.553 + !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
1.554 + _graph.nextOut(e);
1.555 + if (e != INVALID) {
1.556 + stack[++stack_head] = pair(e, true);
1.557 + goto next_step_2;
1.558 + }
1.559 + _graph.firstIn(e, v);
1.560 + while ( e != INVALID && ((*_flow)[e] == 0 ||
1.561 + !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
1.562 + _graph.nextIn(e);
1.563 + stack[++stack_head] = pair(e, false);
1.564 + next_step_2: ;
1.565 + } else {
1.566 + if (!processed[v]) {
1.567 + // A cycle is found
1.568 + Node n, w = u;
1.569 + Capacity d, delta = sf ? _capacity[se] - (*_flow)[se] :
1.570 + (*_flow)[se];
1.571 + for (n = u; n != v; n = pred_node[n]) {
1.572 + d = _graph.target(pred_edge[n]) == n ?
1.573 + _capacity[pred_edge[n]] - (*_flow)[pred_edge[n]] :
1.574 + (*_flow)[pred_edge[n]];
1.575 + if (d <= delta) {
1.576 + delta = d;
1.577 + w = pred_node[n];
1.578 + }
1.579 + }
1.580 +
1.581 +/*
1.582 + std::cout << "CYCLE FOUND: ";
1.583 + if (sf)
1.584 + std::cout << _cost[se] + pi[_graph.source(se)] - pi[_graph.target(se)];
1.585 + else
1.586 + std::cout << _graph.id(se) << ":" << -(_cost[se] + pi[_graph.source(se)] - pi[_graph.target(se)]);
1.587 + for (n = u; n != v; n = pred_node[n]) {
1.588 + if (_graph.target(pred_edge[n]) == n)
1.589 + std::cout << " " << _cost[pred_edge[n]] + pi[_graph.source(pred_edge[n])] - pi[_graph.target(pred_edge[n])];
1.590 + else
1.591 + std::cout << " " << -(_cost[pred_edge[n]] + pi[_graph.source(pred_edge[n])] - pi[_graph.target(pred_edge[n])]);
1.592 + }
1.593 + std::cout << "\n";
1.594 +*/
1.595 + // Augment along the cycle
1.596 + (*_flow)[se] = sf ? (*_flow)[se] + delta :
1.597 + (*_flow)[se] - delta;
1.598 + for (n = u; n != v; n = pred_node[n]) {
1.599 + if (_graph.target(pred_edge[n]) == n)
1.600 + (*_flow)[pred_edge[n]] += delta;
1.601 + else
1.602 + (*_flow)[pred_edge[n]] -= delta;
1.603 + }
1.604 + for (n = u; stack_head > 0 && n != w; n = pred_node[n]) {
1.605 + --stack_head;
1.606 + reached[n] = false;
1.607 + }
1.608 + u = w;
1.609 + }
1.610 + v = u;
1.611 +
1.612 + // Find the next admissible residual outgoing edge
1.613 + double p = pi[v];
1.614 + Edge e = stack[stack_head].first;
1.615 + if (!stack[stack_head].second) {
1.616 + _graph.nextIn(e);
1.617 + goto in_edge_3;
1.618 + }
1.619 + _graph.nextOut(e);
1.620 + while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
1.621 + !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
1.622 + _graph.nextOut(e);
1.623 + if (e != INVALID) {
1.624 + stack[stack_head] = pair(e, true);
1.625 + goto next_step_3;
1.626 + }
1.627 + _graph.firstIn(e, v);
1.628 + in_edge_3:
1.629 + while ( e != INVALID && ((*_flow)[e] == 0 ||
1.630 + !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
1.631 + _graph.nextIn(e);
1.632 + stack[stack_head] = pair(e, false);
1.633 + next_step_3: ;
1.634 + }
1.635 +
1.636 + while (stack_head >= 0 && stack[stack_head].first == INVALID) {
1.637 + processed[v] = true;
1.638 + proc_vector[++proc_head] = v;
1.639 + if (--stack_head >= 0) {
1.640 + v = stack[stack_head].second ?
1.641 + _graph.source(stack[stack_head].first) :
1.642 + _graph.target(stack[stack_head].first);
1.643 + // Find the next admissible residual outgoing edge
1.644 + double p = pi[v];
1.645 + Edge e = stack[stack_head].first;
1.646 + if (!stack[stack_head].second) {
1.647 + _graph.nextIn(e);
1.648 + goto in_edge_4;
1.649 + }
1.650 + _graph.nextOut(e);
1.651 + while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
1.652 + !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
1.653 + _graph.nextOut(e);
1.654 + if (e != INVALID) {
1.655 + stack[stack_head] = pair(e, true);
1.656 + goto next_step_4;
1.657 + }
1.658 + _graph.firstIn(e, v);
1.659 + in_edge_4:
1.660 + while ( e != INVALID && ((*_flow)[e] == 0 ||
1.661 + !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
1.662 + _graph.nextIn(e);
1.663 + stack[stack_head] = pair(e, false);
1.664 + next_step_4: ;
1.665 + }
1.666 + }
1.667 + }
1.668 + }
1.669 + t1.stop();
1.670 +
1.671 + // Tighten potentials and epsilon
1.672 + if (--iter > 0) {
1.673 + // Compute levels
1.674 + t2.start();
1.675 + for (int i = proc_head; i >= 0; --i) {
1.676 + Node v = proc_vector[i];
1.677 + double p = pi[v];
1.678 + int l = 0;
1.679 + for (InEdgeIt e(_graph, v); e != INVALID; ++e) {
1.680 + Node u = _graph.source(e);
1.681 + if ( _capacity[e] - (*_flow)[e] > 0 &&
1.682 + tol.negative(_cost[e] + pi[u] - p) &&
1.683 + level[u] + 1 > l ) l = level[u] + 1;
1.684 + }
1.685 + for (OutEdgeIt e(_graph, v); e != INVALID; ++e) {
1.686 + Node u = _graph.target(e);
1.687 + if ( (*_flow)[e] > 0 &&
1.688 + tol.negative(-_cost[e] + pi[u] - p) &&
1.689 + level[u] + 1 > l ) l = level[u] + 1;
1.690 + }
1.691 + level[v] = l;
1.692 + }
1.693 +
1.694 + // Modify potentials
1.695 + double p, q = -1;
1.696 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.697 + Node u = _graph.source(e);
1.698 + Node v = _graph.target(e);
1.699 + if (_capacity[e] - (*_flow)[e] > 0 && level[u] - level[v] > 0) {
1.700 + p = (_cost[e] + pi[u] - pi[v] + epsilon) /
1.701 + (level[u] - level[v] + 1);
1.702 + if (q < 0 || p < q) q = p;
1.703 + }
1.704 + else if ((*_flow)[e] > 0 && level[v] - level[u] > 0) {
1.705 + p = (-_cost[e] - pi[u] + pi[v] + epsilon) /
1.706 + (level[v] - level[u] + 1);
1.707 + if (q < 0 || p < q) q = p;
1.708 + }
1.709 + }
1.710 + for (NodeIt v(_graph); v != INVALID; ++v) {
1.711 + pi[v] -= q * level[v];
1.712 + }
1.713 +
1.714 + // Modify epsilon
1.715 + epsilon = 0;
1.716 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.717 + double curr = _cost[e] + pi[_graph.source(e)]
1.718 + - pi[_graph.target(e)];
1.719 + if (_capacity[e] - (*_flow)[e] > 0 && curr < -epsilon)
1.720 + epsilon = -curr;
1.721 + else if ((*_flow)[e] > 0 && curr > epsilon)
1.722 + epsilon = curr;
1.723 + }
1.724 + t2.stop();
1.725 + } else {
1.726 + // Set epsilon to the minimum cycle mean
1.727 + t3.start();
1.728 +
1.729 +/**/
1.730 + StaticGraph static_graph;
1.731 + typename ResGraph::template NodeMap<typename StaticGraph::Node> node_ref(*_res_graph);
1.732 + typename ResGraph::template EdgeMap<typename StaticGraph::Edge> edge_ref(*_res_graph);
1.733 + static_graph.build(*_res_graph, node_ref, edge_ref);
1.734 + typename StaticGraph::template NodeMap<double> static_pi(static_graph);
1.735 + typename StaticGraph::template EdgeMap<double> static_cost(static_graph);
1.736 +
1.737 + for (typename ResGraph::EdgeIt e(*_res_graph); e != INVALID; ++e)
1.738 + static_cost[edge_ref[e]] = _res_cost[e];
1.739 +
1.740 + MinMeanCycle<StaticGraph, typename StaticGraph::template EdgeMap<double> >
1.741 + mmc(static_graph, static_cost);
1.742 + mmc.init();
1.743 + mmc.findMinMean();
1.744 + epsilon = -mmc.cycleMean();
1.745 +/**/
1.746 +
1.747 +/*
1.748 + MinMeanCycle<ResGraph, ResidualCostMap> mmc(*_res_graph, _res_cost);
1.749 + mmc.init();
1.750 + mmc.findMinMean();
1.751 + epsilon = -mmc.cycleMean();
1.752 +*/
1.753 +
1.754 + // Compute feasible potentials for the current epsilon
1.755 + for (typename StaticGraph::EdgeIt e(static_graph); e != INVALID; ++e)
1.756 + static_cost[e] += epsilon;
1.757 + typename BellmanFord<StaticGraph, typename StaticGraph::template EdgeMap<double> >::
1.758 + template DefDistMap<typename StaticGraph::template NodeMap<double> >::
1.759 + template DefOperationTraits<BFOperationTraits>::Create
1.760 + bf(static_graph, static_cost);
1.761 + bf.distMap(static_pi).init(0);
1.762 + bf.start();
1.763 + for (NodeIt n(_graph); n != INVALID; ++n)
1.764 + pi[n] = static_pi[node_ref[n]];
1.765 +
1.766 +/*
1.767 + for (typename ResGraph::EdgeIt e(*_res_graph); e != INVALID; ++e)
1.768 + shift_cost[e] = _res_cost[e] + epsilon;
1.769 + typename BellmanFord<ResGraph, ResShiftCostMap>::
1.770 + template DefDistMap<FloatPotentialMap>::
1.771 + template DefOperationTraits<BFOperationTraits>::Create
1.772 + bf(*_res_graph, shift_cost);
1.773 + bf.distMap(pi).init(0);
1.774 + bf.start();
1.775 +*/
1.776 +
1.777 + iter = limit;
1.778 + t3.stop();
1.779 + }
1.780 + }
1.781 +
1.782 +// std::cout << t1.realTime() << " " << t2.realTime() << " " << t3.realTime() << "\n";
1.783 +
1.784 + // Handle non-zero lower bounds
1.785 + if (_lower) {
1.786 + for (EdgeIt e(_graph); e != INVALID; ++e)
1.787 + (*_flow)[e] += (*_lower)[e];
1.788 + }
1.789 + return true;
1.790 + }
1.791 +
1.792 + }; //class CancelAndTighten
1.793 +
1.794 + ///@}
1.795 +
1.796 +} //namespace lemon
1.797 +
1.798 +#endif //LEMON_CANCEL_AND_TIGHTEN_H