1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/cost_scaling.h Thu Nov 12 23:29:42 2009 +0100
1.3 @@ -0,0 +1,850 @@
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_COST_SCALING_H
1.23 +#define LEMON_COST_SCALING_H
1.24 +
1.25 +/// \ingroup min_cost_flow_algs
1.26 +/// \file
1.27 +/// \brief Cost scaling algorithm for finding a minimum cost flow.
1.28 +
1.29 +#include <vector>
1.30 +#include <deque>
1.31 +#include <limits>
1.32 +
1.33 +#include <lemon/core.h>
1.34 +#include <lemon/maps.h>
1.35 +#include <lemon/math.h>
1.36 +#include <lemon/adaptors.h>
1.37 +#include <lemon/circulation.h>
1.38 +#include <lemon/bellman_ford.h>
1.39 +
1.40 +namespace lemon {
1.41 +
1.42 + /// \addtogroup min_cost_flow_algs
1.43 + /// @{
1.44 +
1.45 + /// \brief Implementation of the cost scaling algorithm for finding a
1.46 + /// minimum cost flow.
1.47 + ///
1.48 + /// \ref CostScaling implements the cost scaling algorithm performing
1.49 + /// augment/push and relabel operations for finding a minimum cost
1.50 + /// 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 + /// \note Arc costs are multiplied with the number of nodes during
1.65 + /// the algorithm so overflow problems may arise more easily than with
1.66 + /// other minimum cost flow algorithms.
1.67 + /// If it is available, <tt>long long int</tt> type is used instead of
1.68 + /// <tt>long int</tt> in the inside computations.
1.69 + ///
1.70 + /// \author Peter Kovacs
1.71 + template < typename Digraph,
1.72 + typename LowerMap = typename Digraph::template ArcMap<int>,
1.73 + typename CapacityMap = typename Digraph::template ArcMap<int>,
1.74 + typename CostMap = typename Digraph::template ArcMap<int>,
1.75 + typename SupplyMap = typename Digraph::template NodeMap<int> >
1.76 + class CostScaling
1.77 + {
1.78 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
1.79 +
1.80 + typedef typename CapacityMap::Value Capacity;
1.81 + typedef typename CostMap::Value Cost;
1.82 + typedef typename SupplyMap::Value Supply;
1.83 + typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
1.84 + typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
1.85 +
1.86 + typedef ResidualDigraph< const Digraph,
1.87 + CapacityArcMap, CapacityArcMap > ResDigraph;
1.88 + typedef typename ResDigraph::Arc ResArc;
1.89 +
1.90 +#if defined __GNUC__ && !defined __STRICT_ANSI__
1.91 + typedef long long int LCost;
1.92 +#else
1.93 + typedef long int LCost;
1.94 +#endif
1.95 + typedef typename Digraph::template ArcMap<LCost> LargeCostMap;
1.96 +
1.97 + public:
1.98 +
1.99 + /// The type of the flow map.
1.100 + typedef typename Digraph::template ArcMap<Capacity> FlowMap;
1.101 + /// The type of the potential map.
1.102 + typedef typename Digraph::template NodeMap<LCost> PotentialMap;
1.103 +
1.104 + private:
1.105 +
1.106 + /// \brief Map adaptor class for handling residual arc costs.
1.107 + ///
1.108 + /// Map adaptor class for handling residual arc costs.
1.109 + template <typename Map>
1.110 + class ResidualCostMap : public MapBase<ResArc, typename Map::Value>
1.111 + {
1.112 + private:
1.113 +
1.114 + const Map &_cost_map;
1.115 +
1.116 + public:
1.117 +
1.118 + ///\e
1.119 + ResidualCostMap(const Map &cost_map) :
1.120 + _cost_map(cost_map) {}
1.121 +
1.122 + ///\e
1.123 + inline typename Map::Value operator[](const ResArc &e) const {
1.124 + return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
1.125 + }
1.126 +
1.127 + }; //class ResidualCostMap
1.128 +
1.129 + /// \brief Map adaptor class for handling reduced arc costs.
1.130 + ///
1.131 + /// Map adaptor class for handling reduced arc costs.
1.132 + class ReducedCostMap : public MapBase<Arc, LCost>
1.133 + {
1.134 + private:
1.135 +
1.136 + const Digraph &_gr;
1.137 + const LargeCostMap &_cost_map;
1.138 + const PotentialMap &_pot_map;
1.139 +
1.140 + public:
1.141 +
1.142 + ///\e
1.143 + ReducedCostMap( const Digraph &gr,
1.144 + const LargeCostMap &cost_map,
1.145 + const PotentialMap &pot_map ) :
1.146 + _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
1.147 +
1.148 + ///\e
1.149 + inline LCost operator[](const Arc &e) const {
1.150 + return _cost_map[e] + _pot_map[_gr.source(e)]
1.151 + - _pot_map[_gr.target(e)];
1.152 + }
1.153 +
1.154 + }; //class ReducedCostMap
1.155 +
1.156 + private:
1.157 +
1.158 + // The digraph the algorithm runs on
1.159 + const Digraph &_graph;
1.160 + // The original lower bound map
1.161 + const LowerMap *_lower;
1.162 + // The modified capacity map
1.163 + CapacityArcMap _capacity;
1.164 + // The original cost map
1.165 + const CostMap &_orig_cost;
1.166 + // The scaled cost map
1.167 + LargeCostMap _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 cost map
1.180 + ResidualCostMap<LargeCostMap> _res_cost;
1.181 + // The residual digraph
1.182 + ResDigraph *_res_graph;
1.183 + // The reduced cost map
1.184 + ReducedCostMap *_red_cost;
1.185 + // The excess map
1.186 + SupplyNodeMap _excess;
1.187 + // The epsilon parameter used for cost scaling
1.188 + LCost _epsilon;
1.189 + // The scaling factor
1.190 + int _alpha;
1.191 +
1.192 + public:
1.193 +
1.194 + /// \brief General constructor (with lower bounds).
1.195 + ///
1.196 + /// General constructor (with lower bounds).
1.197 + ///
1.198 + /// \param digraph The digraph the algorithm runs on.
1.199 + /// \param lower The lower bounds of the arcs.
1.200 + /// \param capacity The capacities (upper bounds) of the arcs.
1.201 + /// \param cost The cost (length) values of the arcs.
1.202 + /// \param supply The supply values of the nodes (signed).
1.203 + CostScaling( const Digraph &digraph,
1.204 + const LowerMap &lower,
1.205 + const CapacityMap &capacity,
1.206 + const CostMap &cost,
1.207 + const SupplyMap &supply ) :
1.208 + _graph(digraph), _lower(&lower), _capacity(digraph), _orig_cost(cost),
1.209 + _cost(digraph), _supply(digraph), _flow(NULL), _local_flow(false),
1.210 + _potential(NULL), _local_potential(false), _res_cost(_cost),
1.211 + _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
1.212 + {
1.213 + // Check the sum of supply values
1.214 + Supply sum = 0;
1.215 + for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
1.216 + _valid_supply = sum == 0;
1.217 +
1.218 + for (ArcIt e(_graph); e != INVALID; ++e) _capacity[e] = capacity[e];
1.219 + for (NodeIt n(_graph); n != INVALID; ++n) _supply[n] = supply[n];
1.220 +
1.221 + // Remove non-zero lower bounds
1.222 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.223 + if (lower[e] != 0) {
1.224 + _capacity[e] -= lower[e];
1.225 + _supply[_graph.source(e)] -= lower[e];
1.226 + _supply[_graph.target(e)] += lower[e];
1.227 + }
1.228 + }
1.229 + }
1.230 +/*
1.231 + /// \brief General constructor (without lower bounds).
1.232 + ///
1.233 + /// General constructor (without lower bounds).
1.234 + ///
1.235 + /// \param digraph The digraph the algorithm runs on.
1.236 + /// \param capacity The capacities (upper bounds) of the arcs.
1.237 + /// \param cost The cost (length) values of the arcs.
1.238 + /// \param supply The supply values of the nodes (signed).
1.239 + CostScaling( const Digraph &digraph,
1.240 + const CapacityMap &capacity,
1.241 + const CostMap &cost,
1.242 + const SupplyMap &supply ) :
1.243 + _graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
1.244 + _cost(digraph), _supply(supply), _flow(NULL), _local_flow(false),
1.245 + _potential(NULL), _local_potential(false), _res_cost(_cost),
1.246 + _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
1.247 + {
1.248 + // Check the sum of supply values
1.249 + Supply sum = 0;
1.250 + for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
1.251 + _valid_supply = sum == 0;
1.252 + }
1.253 +
1.254 + /// \brief Simple constructor (with lower bounds).
1.255 + ///
1.256 + /// Simple constructor (with lower bounds).
1.257 + ///
1.258 + /// \param digraph The digraph the algorithm runs on.
1.259 + /// \param lower The lower bounds of the arcs.
1.260 + /// \param capacity The capacities (upper bounds) of the arcs.
1.261 + /// \param cost The cost (length) values of the arcs.
1.262 + /// \param s The source node.
1.263 + /// \param t The target node.
1.264 + /// \param flow_value The required amount of flow from node \c s
1.265 + /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.266 + CostScaling( const Digraph &digraph,
1.267 + const LowerMap &lower,
1.268 + const CapacityMap &capacity,
1.269 + const CostMap &cost,
1.270 + Node s, Node t,
1.271 + Supply flow_value ) :
1.272 + _graph(digraph), _lower(&lower), _capacity(capacity), _orig_cost(cost),
1.273 + _cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false),
1.274 + _potential(NULL), _local_potential(false), _res_cost(_cost),
1.275 + _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
1.276 + {
1.277 + // Remove non-zero lower bounds
1.278 + _supply[s] = flow_value;
1.279 + _supply[t] = -flow_value;
1.280 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.281 + if (lower[e] != 0) {
1.282 + _capacity[e] -= lower[e];
1.283 + _supply[_graph.source(e)] -= lower[e];
1.284 + _supply[_graph.target(e)] += lower[e];
1.285 + }
1.286 + }
1.287 + _valid_supply = true;
1.288 + }
1.289 +
1.290 + /// \brief Simple constructor (without lower bounds).
1.291 + ///
1.292 + /// Simple constructor (without lower bounds).
1.293 + ///
1.294 + /// \param digraph The digraph the algorithm runs on.
1.295 + /// \param capacity The capacities (upper bounds) of the arcs.
1.296 + /// \param cost The cost (length) values of the arcs.
1.297 + /// \param s The source node.
1.298 + /// \param t The target node.
1.299 + /// \param flow_value The required amount of flow from node \c s
1.300 + /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.301 + CostScaling( const Digraph &digraph,
1.302 + const CapacityMap &capacity,
1.303 + const CostMap &cost,
1.304 + Node s, Node t,
1.305 + Supply flow_value ) :
1.306 + _graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
1.307 + _cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false),
1.308 + _potential(NULL), _local_potential(false), _res_cost(_cost),
1.309 + _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
1.310 + {
1.311 + _supply[s] = flow_value;
1.312 + _supply[t] = -flow_value;
1.313 + _valid_supply = true;
1.314 + }
1.315 +*/
1.316 + /// Destructor.
1.317 + ~CostScaling() {
1.318 + if (_local_flow) delete _flow;
1.319 + if (_local_potential) delete _potential;
1.320 + delete _res_graph;
1.321 + delete _red_cost;
1.322 + }
1.323 +
1.324 + /// \brief Set the flow map.
1.325 + ///
1.326 + /// Set the flow map.
1.327 + ///
1.328 + /// \return \c (*this)
1.329 + CostScaling& flowMap(FlowMap &map) {
1.330 + if (_local_flow) {
1.331 + delete _flow;
1.332 + _local_flow = false;
1.333 + }
1.334 + _flow = ↦
1.335 + return *this;
1.336 + }
1.337 +
1.338 + /// \brief Set the potential map.
1.339 + ///
1.340 + /// Set the potential map.
1.341 + ///
1.342 + /// \return \c (*this)
1.343 + CostScaling& potentialMap(PotentialMap &map) {
1.344 + if (_local_potential) {
1.345 + delete _potential;
1.346 + _local_potential = false;
1.347 + }
1.348 + _potential = ↦
1.349 + return *this;
1.350 + }
1.351 +
1.352 + /// \name Execution control
1.353 +
1.354 + /// @{
1.355 +
1.356 + /// \brief Run the algorithm.
1.357 + ///
1.358 + /// Run the algorithm.
1.359 + ///
1.360 + /// \param partial_augment By default the algorithm performs
1.361 + /// partial augment and relabel operations in the cost scaling
1.362 + /// phases. Set this parameter to \c false for using local push and
1.363 + /// relabel operations instead.
1.364 + ///
1.365 + /// \return \c true if a feasible flow can be found.
1.366 + bool run(bool partial_augment = true) {
1.367 + if (partial_augment) {
1.368 + return init() && startPartialAugment();
1.369 + } else {
1.370 + return init() && startPushRelabel();
1.371 + }
1.372 + }
1.373 +
1.374 + /// @}
1.375 +
1.376 + /// \name Query Functions
1.377 + /// The result of the algorithm can be obtained using these
1.378 + /// functions.\n
1.379 + /// \ref lemon::CostScaling::run() "run()" must be called before
1.380 + /// using them.
1.381 +
1.382 + /// @{
1.383 +
1.384 + /// \brief Return a const reference to the arc map storing the
1.385 + /// found flow.
1.386 + ///
1.387 + /// Return a const reference to the arc map storing the found flow.
1.388 + ///
1.389 + /// \pre \ref run() must be called before using this function.
1.390 + const FlowMap& flowMap() const {
1.391 + return *_flow;
1.392 + }
1.393 +
1.394 + /// \brief Return a const reference to the node map storing the
1.395 + /// found potentials (the dual solution).
1.396 + ///
1.397 + /// Return a const reference to the node map storing the found
1.398 + /// potentials (the dual solution).
1.399 + ///
1.400 + /// \pre \ref run() must be called before using this function.
1.401 + const PotentialMap& potentialMap() const {
1.402 + return *_potential;
1.403 + }
1.404 +
1.405 + /// \brief Return the flow on the given arc.
1.406 + ///
1.407 + /// Return the flow on the given arc.
1.408 + ///
1.409 + /// \pre \ref run() must be called before using this function.
1.410 + Capacity flow(const Arc& arc) const {
1.411 + return (*_flow)[arc];
1.412 + }
1.413 +
1.414 + /// \brief Return the potential of the given node.
1.415 + ///
1.416 + /// Return the potential of the given node.
1.417 + ///
1.418 + /// \pre \ref run() must be called before using this function.
1.419 + Cost potential(const Node& node) const {
1.420 + return (*_potential)[node];
1.421 + }
1.422 +
1.423 + /// \brief Return the total cost of the found flow.
1.424 + ///
1.425 + /// Return the total cost of the found flow. The complexity of the
1.426 + /// function is \f$ O(e) \f$.
1.427 + ///
1.428 + /// \pre \ref run() must be called before using this function.
1.429 + Cost totalCost() const {
1.430 + Cost c = 0;
1.431 + for (ArcIt e(_graph); e != INVALID; ++e)
1.432 + c += (*_flow)[e] * _orig_cost[e];
1.433 + return c;
1.434 + }
1.435 +
1.436 + /// @}
1.437 +
1.438 + private:
1.439 +
1.440 + /// Initialize the algorithm.
1.441 + bool init() {
1.442 + if (!_valid_supply) return false;
1.443 + // The scaling factor
1.444 + _alpha = 8;
1.445 +
1.446 + // Initialize flow and potential maps
1.447 + if (!_flow) {
1.448 + _flow = new FlowMap(_graph);
1.449 + _local_flow = true;
1.450 + }
1.451 + if (!_potential) {
1.452 + _potential = new PotentialMap(_graph);
1.453 + _local_potential = true;
1.454 + }
1.455 +
1.456 + _red_cost = new ReducedCostMap(_graph, _cost, *_potential);
1.457 + _res_graph = new ResDigraph(_graph, _capacity, *_flow);
1.458 +
1.459 + // Initialize the scaled cost map and the epsilon parameter
1.460 + Cost max_cost = 0;
1.461 + int node_num = countNodes(_graph);
1.462 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.463 + _cost[e] = LCost(_orig_cost[e]) * node_num * _alpha;
1.464 + if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e];
1.465 + }
1.466 + _epsilon = max_cost * node_num;
1.467 +
1.468 + // Find a feasible flow using Circulation
1.469 + Circulation< Digraph, ConstMap<Arc, Capacity>, CapacityArcMap,
1.470 + SupplyMap >
1.471 + circulation( _graph, constMap<Arc>(Capacity(0)), _capacity,
1.472 + _supply );
1.473 + return circulation.flowMap(*_flow).run();
1.474 + }
1.475 +
1.476 + /// Execute the algorithm performing partial augmentation and
1.477 + /// relabel operations.
1.478 + bool startPartialAugment() {
1.479 + // Paramters for heuristics
1.480 +// const int BF_HEURISTIC_EPSILON_BOUND = 1000;
1.481 +// const int BF_HEURISTIC_BOUND_FACTOR = 3;
1.482 + // Maximum augment path length
1.483 + const int MAX_PATH_LENGTH = 4;
1.484 +
1.485 + // Variables
1.486 + typename Digraph::template NodeMap<Arc> pred_arc(_graph);
1.487 + typename Digraph::template NodeMap<bool> forward(_graph);
1.488 + typename Digraph::template NodeMap<OutArcIt> next_out(_graph);
1.489 + typename Digraph::template NodeMap<InArcIt> next_in(_graph);
1.490 + typename Digraph::template NodeMap<bool> next_dir(_graph);
1.491 + std::deque<Node> active_nodes;
1.492 + std::vector<Node> path_nodes;
1.493 +
1.494 +// int node_num = countNodes(_graph);
1.495 + for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.496 + 1 : _epsilon / _alpha )
1.497 + {
1.498 +/*
1.499 + // "Early Termination" heuristic: use Bellman-Ford algorithm
1.500 + // to check if the current flow is optimal
1.501 + if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
1.502 + typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
1.503 + ShiftCostMap shift_cost(_res_cost, 1);
1.504 + BellmanFord<ResDigraph, ShiftCostMap> bf(*_res_graph, shift_cost);
1.505 + bf.init(0);
1.506 + bool done = false;
1.507 + int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
1.508 + for (int i = 0; i < K && !done; ++i)
1.509 + done = bf.processNextWeakRound();
1.510 + if (done) break;
1.511 + }
1.512 +*/
1.513 + // Saturate arcs not satisfying the optimality condition
1.514 + Capacity delta;
1.515 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.516 + if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.517 + delta = _capacity[e] - (*_flow)[e];
1.518 + _excess[_graph.source(e)] -= delta;
1.519 + _excess[_graph.target(e)] += delta;
1.520 + (*_flow)[e] = _capacity[e];
1.521 + }
1.522 + if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
1.523 + _excess[_graph.target(e)] -= (*_flow)[e];
1.524 + _excess[_graph.source(e)] += (*_flow)[e];
1.525 + (*_flow)[e] = 0;
1.526 + }
1.527 + }
1.528 +
1.529 + // Find active nodes (i.e. nodes with positive excess)
1.530 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.531 + if (_excess[n] > 0) active_nodes.push_back(n);
1.532 + }
1.533 +
1.534 + // Initialize the next arc maps
1.535 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.536 + next_out[n] = OutArcIt(_graph, n);
1.537 + next_in[n] = InArcIt(_graph, n);
1.538 + next_dir[n] = true;
1.539 + }
1.540 +
1.541 + // Perform partial augment and relabel operations
1.542 + while (active_nodes.size() > 0) {
1.543 + // Select an active node (FIFO selection)
1.544 + if (_excess[active_nodes[0]] <= 0) {
1.545 + active_nodes.pop_front();
1.546 + continue;
1.547 + }
1.548 + Node start = active_nodes[0];
1.549 + path_nodes.clear();
1.550 + path_nodes.push_back(start);
1.551 +
1.552 + // Find an augmenting path from the start node
1.553 + Node u, tip = start;
1.554 + LCost min_red_cost;
1.555 + while ( _excess[tip] >= 0 &&
1.556 + int(path_nodes.size()) <= MAX_PATH_LENGTH )
1.557 + {
1.558 + if (next_dir[tip]) {
1.559 + for (OutArcIt e = next_out[tip]; e != INVALID; ++e) {
1.560 + if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.561 + u = _graph.target(e);
1.562 + pred_arc[u] = e;
1.563 + forward[u] = true;
1.564 + next_out[tip] = e;
1.565 + tip = u;
1.566 + path_nodes.push_back(tip);
1.567 + goto next_step;
1.568 + }
1.569 + }
1.570 + next_dir[tip] = false;
1.571 + }
1.572 + for (InArcIt e = next_in[tip]; e != INVALID; ++e) {
1.573 + if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
1.574 + u = _graph.source(e);
1.575 + pred_arc[u] = e;
1.576 + forward[u] = false;
1.577 + next_in[tip] = e;
1.578 + tip = u;
1.579 + path_nodes.push_back(tip);
1.580 + goto next_step;
1.581 + }
1.582 + }
1.583 +
1.584 + // Relabel tip node
1.585 + min_red_cost = std::numeric_limits<LCost>::max() / 2;
1.586 + for (OutArcIt oe(_graph, tip); oe != INVALID; ++oe) {
1.587 + if ( _capacity[oe] - (*_flow)[oe] > 0 &&
1.588 + (*_red_cost)[oe] < min_red_cost )
1.589 + min_red_cost = (*_red_cost)[oe];
1.590 + }
1.591 + for (InArcIt ie(_graph, tip); ie != INVALID; ++ie) {
1.592 + if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
1.593 + min_red_cost = -(*_red_cost)[ie];
1.594 + }
1.595 + (*_potential)[tip] -= min_red_cost + _epsilon;
1.596 +
1.597 + // Reset the next arc maps
1.598 + next_out[tip] = OutArcIt(_graph, tip);
1.599 + next_in[tip] = InArcIt(_graph, tip);
1.600 + next_dir[tip] = true;
1.601 +
1.602 + // Step back
1.603 + if (tip != start) {
1.604 + path_nodes.pop_back();
1.605 + tip = path_nodes[path_nodes.size()-1];
1.606 + }
1.607 +
1.608 + next_step:
1.609 + continue;
1.610 + }
1.611 +
1.612 + // Augment along the found path (as much flow as possible)
1.613 + Capacity delta;
1.614 + for (int i = 1; i < int(path_nodes.size()); ++i) {
1.615 + u = path_nodes[i];
1.616 + delta = forward[u] ?
1.617 + _capacity[pred_arc[u]] - (*_flow)[pred_arc[u]] :
1.618 + (*_flow)[pred_arc[u]];
1.619 + delta = std::min(delta, _excess[path_nodes[i-1]]);
1.620 + (*_flow)[pred_arc[u]] += forward[u] ? delta : -delta;
1.621 + _excess[path_nodes[i-1]] -= delta;
1.622 + _excess[u] += delta;
1.623 + if (_excess[u] > 0 && _excess[u] <= delta) active_nodes.push_back(u);
1.624 + }
1.625 + }
1.626 + }
1.627 +
1.628 + // Compute node potentials for the original costs
1.629 + ResidualCostMap<CostMap> res_cost(_orig_cost);
1.630 + BellmanFord< ResDigraph, ResidualCostMap<CostMap> >
1.631 + bf(*_res_graph, res_cost);
1.632 + bf.init(0); bf.start();
1.633 + for (NodeIt n(_graph); n != INVALID; ++n)
1.634 + (*_potential)[n] = bf.dist(n);
1.635 +
1.636 + // Handle non-zero lower bounds
1.637 + if (_lower) {
1.638 + for (ArcIt e(_graph); e != INVALID; ++e)
1.639 + (*_flow)[e] += (*_lower)[e];
1.640 + }
1.641 + return true;
1.642 + }
1.643 +
1.644 + /// Execute the algorithm performing push and relabel operations.
1.645 + bool startPushRelabel() {
1.646 + // Paramters for heuristics
1.647 +// const int BF_HEURISTIC_EPSILON_BOUND = 1000;
1.648 +// const int BF_HEURISTIC_BOUND_FACTOR = 3;
1.649 +
1.650 + typename Digraph::template NodeMap<bool> hyper(_graph, false);
1.651 + typename Digraph::template NodeMap<Arc> pred_arc(_graph);
1.652 + typename Digraph::template NodeMap<bool> forward(_graph);
1.653 + typename Digraph::template NodeMap<OutArcIt> next_out(_graph);
1.654 + typename Digraph::template NodeMap<InArcIt> next_in(_graph);
1.655 + typename Digraph::template NodeMap<bool> next_dir(_graph);
1.656 + std::deque<Node> active_nodes;
1.657 +
1.658 +// int node_num = countNodes(_graph);
1.659 + for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.660 + 1 : _epsilon / _alpha )
1.661 + {
1.662 +/*
1.663 + // "Early Termination" heuristic: use Bellman-Ford algorithm
1.664 + // to check if the current flow is optimal
1.665 + if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
1.666 + typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
1.667 + ShiftCostMap shift_cost(_res_cost, 1);
1.668 + BellmanFord<ResDigraph, ShiftCostMap> bf(*_res_graph, shift_cost);
1.669 + bf.init(0);
1.670 + bool done = false;
1.671 + int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
1.672 + for (int i = 0; i < K && !done; ++i)
1.673 + done = bf.processNextWeakRound();
1.674 + if (done) break;
1.675 + }
1.676 +*/
1.677 +
1.678 + // Saturate arcs not satisfying the optimality condition
1.679 + Capacity delta;
1.680 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.681 + if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.682 + delta = _capacity[e] - (*_flow)[e];
1.683 + _excess[_graph.source(e)] -= delta;
1.684 + _excess[_graph.target(e)] += delta;
1.685 + (*_flow)[e] = _capacity[e];
1.686 + }
1.687 + if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
1.688 + _excess[_graph.target(e)] -= (*_flow)[e];
1.689 + _excess[_graph.source(e)] += (*_flow)[e];
1.690 + (*_flow)[e] = 0;
1.691 + }
1.692 + }
1.693 +
1.694 + // Find active nodes (i.e. nodes with positive excess)
1.695 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.696 + if (_excess[n] > 0) active_nodes.push_back(n);
1.697 + }
1.698 +
1.699 + // Initialize the next arc maps
1.700 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.701 + next_out[n] = OutArcIt(_graph, n);
1.702 + next_in[n] = InArcIt(_graph, n);
1.703 + next_dir[n] = true;
1.704 + }
1.705 +
1.706 + // Perform push and relabel operations
1.707 + while (active_nodes.size() > 0) {
1.708 + // Select an active node (FIFO selection)
1.709 + Node n = active_nodes[0], t;
1.710 + bool relabel_enabled = true;
1.711 +
1.712 + // Perform push operations if there are admissible arcs
1.713 + if (_excess[n] > 0 && next_dir[n]) {
1.714 + OutArcIt e = next_out[n];
1.715 + for ( ; e != INVALID; ++e) {
1.716 + if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.717 + delta = std::min(_capacity[e] - (*_flow)[e], _excess[n]);
1.718 + t = _graph.target(e);
1.719 +
1.720 + // Push-look-ahead heuristic
1.721 + Capacity ahead = -_excess[t];
1.722 + for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) {
1.723 + if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
1.724 + ahead += _capacity[oe] - (*_flow)[oe];
1.725 + }
1.726 + for (InArcIt ie(_graph, t); ie != INVALID; ++ie) {
1.727 + if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
1.728 + ahead += (*_flow)[ie];
1.729 + }
1.730 + if (ahead < 0) ahead = 0;
1.731 +
1.732 + // Push flow along the arc
1.733 + if (ahead < delta) {
1.734 + (*_flow)[e] += ahead;
1.735 + _excess[n] -= ahead;
1.736 + _excess[t] += ahead;
1.737 + active_nodes.push_front(t);
1.738 + hyper[t] = true;
1.739 + relabel_enabled = false;
1.740 + break;
1.741 + } else {
1.742 + (*_flow)[e] += delta;
1.743 + _excess[n] -= delta;
1.744 + _excess[t] += delta;
1.745 + if (_excess[t] > 0 && _excess[t] <= delta)
1.746 + active_nodes.push_back(t);
1.747 + }
1.748 +
1.749 + if (_excess[n] == 0) break;
1.750 + }
1.751 + }
1.752 + if (e != INVALID) {
1.753 + next_out[n] = e;
1.754 + } else {
1.755 + next_dir[n] = false;
1.756 + }
1.757 + }
1.758 +
1.759 + if (_excess[n] > 0 && !next_dir[n]) {
1.760 + InArcIt e = next_in[n];
1.761 + for ( ; e != INVALID; ++e) {
1.762 + if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
1.763 + delta = std::min((*_flow)[e], _excess[n]);
1.764 + t = _graph.source(e);
1.765 +
1.766 + // Push-look-ahead heuristic
1.767 + Capacity ahead = -_excess[t];
1.768 + for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) {
1.769 + if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
1.770 + ahead += _capacity[oe] - (*_flow)[oe];
1.771 + }
1.772 + for (InArcIt ie(_graph, t); ie != INVALID; ++ie) {
1.773 + if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
1.774 + ahead += (*_flow)[ie];
1.775 + }
1.776 + if (ahead < 0) ahead = 0;
1.777 +
1.778 + // Push flow along the arc
1.779 + if (ahead < delta) {
1.780 + (*_flow)[e] -= ahead;
1.781 + _excess[n] -= ahead;
1.782 + _excess[t] += ahead;
1.783 + active_nodes.push_front(t);
1.784 + hyper[t] = true;
1.785 + relabel_enabled = false;
1.786 + break;
1.787 + } else {
1.788 + (*_flow)[e] -= delta;
1.789 + _excess[n] -= delta;
1.790 + _excess[t] += delta;
1.791 + if (_excess[t] > 0 && _excess[t] <= delta)
1.792 + active_nodes.push_back(t);
1.793 + }
1.794 +
1.795 + if (_excess[n] == 0) break;
1.796 + }
1.797 + }
1.798 + next_in[n] = e;
1.799 + }
1.800 +
1.801 + // Relabel the node if it is still active (or hyper)
1.802 + if (relabel_enabled && (_excess[n] > 0 || hyper[n])) {
1.803 + LCost min_red_cost = std::numeric_limits<LCost>::max() / 2;
1.804 + for (OutArcIt oe(_graph, n); oe != INVALID; ++oe) {
1.805 + if ( _capacity[oe] - (*_flow)[oe] > 0 &&
1.806 + (*_red_cost)[oe] < min_red_cost )
1.807 + min_red_cost = (*_red_cost)[oe];
1.808 + }
1.809 + for (InArcIt ie(_graph, n); ie != INVALID; ++ie) {
1.810 + if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
1.811 + min_red_cost = -(*_red_cost)[ie];
1.812 + }
1.813 + (*_potential)[n] -= min_red_cost + _epsilon;
1.814 + hyper[n] = false;
1.815 +
1.816 + // Reset the next arc maps
1.817 + next_out[n] = OutArcIt(_graph, n);
1.818 + next_in[n] = InArcIt(_graph, n);
1.819 + next_dir[n] = true;
1.820 + }
1.821 +
1.822 + // Remove nodes that are not active nor hyper
1.823 + while ( active_nodes.size() > 0 &&
1.824 + _excess[active_nodes[0]] <= 0 &&
1.825 + !hyper[active_nodes[0]] ) {
1.826 + active_nodes.pop_front();
1.827 + }
1.828 + }
1.829 + }
1.830 +
1.831 + // Compute node potentials for the original costs
1.832 + ResidualCostMap<CostMap> res_cost(_orig_cost);
1.833 + BellmanFord< ResDigraph, ResidualCostMap<CostMap> >
1.834 + bf(*_res_graph, res_cost);
1.835 + bf.init(0); bf.start();
1.836 + for (NodeIt n(_graph); n != INVALID; ++n)
1.837 + (*_potential)[n] = bf.dist(n);
1.838 +
1.839 + // Handle non-zero lower bounds
1.840 + if (_lower) {
1.841 + for (ArcIt e(_graph); e != INVALID; ++e)
1.842 + (*_flow)[e] += (*_lower)[e];
1.843 + }
1.844 + return true;
1.845 + }
1.846 +
1.847 + }; //class CostScaling
1.848 +
1.849 + ///@}
1.850 +
1.851 +} //namespace lemon
1.852 +
1.853 +#endif //LEMON_COST_SCALING_H