1.1 --- a/lemon/Makefile.am Thu Nov 12 23:27:21 2009 +0100
1.2 +++ b/lemon/Makefile.am Thu Nov 12 23:29:42 2009 +0100
1.3 @@ -69,8 +69,9 @@
1.4 lemon/color.h \
1.5 lemon/concept_check.h \
1.6 lemon/connectivity.h \
1.7 + lemon/core.h \
1.8 + lemon/cost_scaling.h \
1.9 lemon/counter.h \
1.10 - lemon/core.h \
1.11 lemon/cplex.h \
1.12 lemon/dfs.h \
1.13 lemon/dijkstra.h \
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/lemon/cost_scaling.h Thu Nov 12 23:29:42 2009 +0100
2.3 @@ -0,0 +1,850 @@
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_COST_SCALING_H
2.23 +#define LEMON_COST_SCALING_H
2.24 +
2.25 +/// \ingroup min_cost_flow_algs
2.26 +/// \file
2.27 +/// \brief Cost scaling algorithm for finding a minimum cost flow.
2.28 +
2.29 +#include <vector>
2.30 +#include <deque>
2.31 +#include <limits>
2.32 +
2.33 +#include <lemon/core.h>
2.34 +#include <lemon/maps.h>
2.35 +#include <lemon/math.h>
2.36 +#include <lemon/adaptors.h>
2.37 +#include <lemon/circulation.h>
2.38 +#include <lemon/bellman_ford.h>
2.39 +
2.40 +namespace lemon {
2.41 +
2.42 + /// \addtogroup min_cost_flow_algs
2.43 + /// @{
2.44 +
2.45 + /// \brief Implementation of the cost scaling algorithm for finding a
2.46 + /// minimum cost flow.
2.47 + ///
2.48 + /// \ref CostScaling implements the cost scaling algorithm performing
2.49 + /// augment/push and relabel operations for finding a minimum cost
2.50 + /// 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 + /// \note Arc costs are multiplied with the number of nodes during
2.65 + /// the algorithm so overflow problems may arise more easily than with
2.66 + /// other minimum cost flow algorithms.
2.67 + /// If it is available, <tt>long long int</tt> type is used instead of
2.68 + /// <tt>long int</tt> in the inside computations.
2.69 + ///
2.70 + /// \author Peter Kovacs
2.71 + template < typename Digraph,
2.72 + typename LowerMap = typename Digraph::template ArcMap<int>,
2.73 + typename CapacityMap = typename Digraph::template ArcMap<int>,
2.74 + typename CostMap = typename Digraph::template ArcMap<int>,
2.75 + typename SupplyMap = typename Digraph::template NodeMap<int> >
2.76 + class CostScaling
2.77 + {
2.78 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
2.79 +
2.80 + typedef typename CapacityMap::Value Capacity;
2.81 + typedef typename CostMap::Value Cost;
2.82 + typedef typename SupplyMap::Value Supply;
2.83 + typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
2.84 + typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
2.85 +
2.86 + typedef ResidualDigraph< const Digraph,
2.87 + CapacityArcMap, CapacityArcMap > ResDigraph;
2.88 + typedef typename ResDigraph::Arc ResArc;
2.89 +
2.90 +#if defined __GNUC__ && !defined __STRICT_ANSI__
2.91 + typedef long long int LCost;
2.92 +#else
2.93 + typedef long int LCost;
2.94 +#endif
2.95 + typedef typename Digraph::template ArcMap<LCost> LargeCostMap;
2.96 +
2.97 + public:
2.98 +
2.99 + /// The type of the flow map.
2.100 + typedef typename Digraph::template ArcMap<Capacity> FlowMap;
2.101 + /// The type of the potential map.
2.102 + typedef typename Digraph::template NodeMap<LCost> PotentialMap;
2.103 +
2.104 + private:
2.105 +
2.106 + /// \brief Map adaptor class for handling residual arc costs.
2.107 + ///
2.108 + /// Map adaptor class for handling residual arc costs.
2.109 + template <typename Map>
2.110 + class ResidualCostMap : public MapBase<ResArc, typename Map::Value>
2.111 + {
2.112 + private:
2.113 +
2.114 + const Map &_cost_map;
2.115 +
2.116 + public:
2.117 +
2.118 + ///\e
2.119 + ResidualCostMap(const Map &cost_map) :
2.120 + _cost_map(cost_map) {}
2.121 +
2.122 + ///\e
2.123 + inline typename Map::Value operator[](const ResArc &e) const {
2.124 + return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
2.125 + }
2.126 +
2.127 + }; //class ResidualCostMap
2.128 +
2.129 + /// \brief Map adaptor class for handling reduced arc costs.
2.130 + ///
2.131 + /// Map adaptor class for handling reduced arc costs.
2.132 + class ReducedCostMap : public MapBase<Arc, LCost>
2.133 + {
2.134 + private:
2.135 +
2.136 + const Digraph &_gr;
2.137 + const LargeCostMap &_cost_map;
2.138 + const PotentialMap &_pot_map;
2.139 +
2.140 + public:
2.141 +
2.142 + ///\e
2.143 + ReducedCostMap( const Digraph &gr,
2.144 + const LargeCostMap &cost_map,
2.145 + const PotentialMap &pot_map ) :
2.146 + _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
2.147 +
2.148 + ///\e
2.149 + inline LCost operator[](const Arc &e) const {
2.150 + return _cost_map[e] + _pot_map[_gr.source(e)]
2.151 + - _pot_map[_gr.target(e)];
2.152 + }
2.153 +
2.154 + }; //class ReducedCostMap
2.155 +
2.156 + private:
2.157 +
2.158 + // The digraph the algorithm runs on
2.159 + const Digraph &_graph;
2.160 + // The original lower bound map
2.161 + const LowerMap *_lower;
2.162 + // The modified capacity map
2.163 + CapacityArcMap _capacity;
2.164 + // The original cost map
2.165 + const CostMap &_orig_cost;
2.166 + // The scaled cost map
2.167 + LargeCostMap _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 cost map
2.180 + ResidualCostMap<LargeCostMap> _res_cost;
2.181 + // The residual digraph
2.182 + ResDigraph *_res_graph;
2.183 + // The reduced cost map
2.184 + ReducedCostMap *_red_cost;
2.185 + // The excess map
2.186 + SupplyNodeMap _excess;
2.187 + // The epsilon parameter used for cost scaling
2.188 + LCost _epsilon;
2.189 + // The scaling factor
2.190 + int _alpha;
2.191 +
2.192 + public:
2.193 +
2.194 + /// \brief General constructor (with lower bounds).
2.195 + ///
2.196 + /// General constructor (with lower bounds).
2.197 + ///
2.198 + /// \param digraph The digraph the algorithm runs on.
2.199 + /// \param lower The lower bounds of the arcs.
2.200 + /// \param capacity The capacities (upper bounds) of the arcs.
2.201 + /// \param cost The cost (length) values of the arcs.
2.202 + /// \param supply The supply values of the nodes (signed).
2.203 + CostScaling( const Digraph &digraph,
2.204 + const LowerMap &lower,
2.205 + const CapacityMap &capacity,
2.206 + const CostMap &cost,
2.207 + const SupplyMap &supply ) :
2.208 + _graph(digraph), _lower(&lower), _capacity(digraph), _orig_cost(cost),
2.209 + _cost(digraph), _supply(digraph), _flow(NULL), _local_flow(false),
2.210 + _potential(NULL), _local_potential(false), _res_cost(_cost),
2.211 + _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
2.212 + {
2.213 + // Check the sum of supply values
2.214 + Supply sum = 0;
2.215 + for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
2.216 + _valid_supply = sum == 0;
2.217 +
2.218 + for (ArcIt e(_graph); e != INVALID; ++e) _capacity[e] = capacity[e];
2.219 + for (NodeIt n(_graph); n != INVALID; ++n) _supply[n] = supply[n];
2.220 +
2.221 + // Remove non-zero lower bounds
2.222 + for (ArcIt e(_graph); e != INVALID; ++e) {
2.223 + if (lower[e] != 0) {
2.224 + _capacity[e] -= lower[e];
2.225 + _supply[_graph.source(e)] -= lower[e];
2.226 + _supply[_graph.target(e)] += lower[e];
2.227 + }
2.228 + }
2.229 + }
2.230 +/*
2.231 + /// \brief General constructor (without lower bounds).
2.232 + ///
2.233 + /// General constructor (without lower bounds).
2.234 + ///
2.235 + /// \param digraph The digraph the algorithm runs on.
2.236 + /// \param capacity The capacities (upper bounds) of the arcs.
2.237 + /// \param cost The cost (length) values of the arcs.
2.238 + /// \param supply The supply values of the nodes (signed).
2.239 + CostScaling( const Digraph &digraph,
2.240 + const CapacityMap &capacity,
2.241 + const CostMap &cost,
2.242 + const SupplyMap &supply ) :
2.243 + _graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
2.244 + _cost(digraph), _supply(supply), _flow(NULL), _local_flow(false),
2.245 + _potential(NULL), _local_potential(false), _res_cost(_cost),
2.246 + _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
2.247 + {
2.248 + // Check the sum of supply values
2.249 + Supply sum = 0;
2.250 + for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
2.251 + _valid_supply = sum == 0;
2.252 + }
2.253 +
2.254 + /// \brief Simple constructor (with lower bounds).
2.255 + ///
2.256 + /// Simple constructor (with lower bounds).
2.257 + ///
2.258 + /// \param digraph The digraph the algorithm runs on.
2.259 + /// \param lower The lower bounds of the arcs.
2.260 + /// \param capacity The capacities (upper bounds) of the arcs.
2.261 + /// \param cost The cost (length) values of the arcs.
2.262 + /// \param s The source node.
2.263 + /// \param t The target node.
2.264 + /// \param flow_value The required amount of flow from node \c s
2.265 + /// to node \c t (i.e. the supply of \c s and the demand of \c t).
2.266 + CostScaling( const Digraph &digraph,
2.267 + const LowerMap &lower,
2.268 + const CapacityMap &capacity,
2.269 + const CostMap &cost,
2.270 + Node s, Node t,
2.271 + Supply flow_value ) :
2.272 + _graph(digraph), _lower(&lower), _capacity(capacity), _orig_cost(cost),
2.273 + _cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false),
2.274 + _potential(NULL), _local_potential(false), _res_cost(_cost),
2.275 + _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
2.276 + {
2.277 + // Remove non-zero lower bounds
2.278 + _supply[s] = flow_value;
2.279 + _supply[t] = -flow_value;
2.280 + for (ArcIt e(_graph); e != INVALID; ++e) {
2.281 + if (lower[e] != 0) {
2.282 + _capacity[e] -= lower[e];
2.283 + _supply[_graph.source(e)] -= lower[e];
2.284 + _supply[_graph.target(e)] += lower[e];
2.285 + }
2.286 + }
2.287 + _valid_supply = true;
2.288 + }
2.289 +
2.290 + /// \brief Simple constructor (without lower bounds).
2.291 + ///
2.292 + /// Simple constructor (without lower bounds).
2.293 + ///
2.294 + /// \param digraph The digraph the algorithm runs on.
2.295 + /// \param capacity The capacities (upper bounds) of the arcs.
2.296 + /// \param cost The cost (length) values of the arcs.
2.297 + /// \param s The source node.
2.298 + /// \param t The target node.
2.299 + /// \param flow_value The required amount of flow from node \c s
2.300 + /// to node \c t (i.e. the supply of \c s and the demand of \c t).
2.301 + CostScaling( const Digraph &digraph,
2.302 + const CapacityMap &capacity,
2.303 + const CostMap &cost,
2.304 + Node s, Node t,
2.305 + Supply flow_value ) :
2.306 + _graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
2.307 + _cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false),
2.308 + _potential(NULL), _local_potential(false), _res_cost(_cost),
2.309 + _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
2.310 + {
2.311 + _supply[s] = flow_value;
2.312 + _supply[t] = -flow_value;
2.313 + _valid_supply = true;
2.314 + }
2.315 +*/
2.316 + /// Destructor.
2.317 + ~CostScaling() {
2.318 + if (_local_flow) delete _flow;
2.319 + if (_local_potential) delete _potential;
2.320 + delete _res_graph;
2.321 + delete _red_cost;
2.322 + }
2.323 +
2.324 + /// \brief Set the flow map.
2.325 + ///
2.326 + /// Set the flow map.
2.327 + ///
2.328 + /// \return \c (*this)
2.329 + CostScaling& flowMap(FlowMap &map) {
2.330 + if (_local_flow) {
2.331 + delete _flow;
2.332 + _local_flow = false;
2.333 + }
2.334 + _flow = ↦
2.335 + return *this;
2.336 + }
2.337 +
2.338 + /// \brief Set the potential map.
2.339 + ///
2.340 + /// Set the potential map.
2.341 + ///
2.342 + /// \return \c (*this)
2.343 + CostScaling& potentialMap(PotentialMap &map) {
2.344 + if (_local_potential) {
2.345 + delete _potential;
2.346 + _local_potential = false;
2.347 + }
2.348 + _potential = ↦
2.349 + return *this;
2.350 + }
2.351 +
2.352 + /// \name Execution control
2.353 +
2.354 + /// @{
2.355 +
2.356 + /// \brief Run the algorithm.
2.357 + ///
2.358 + /// Run the algorithm.
2.359 + ///
2.360 + /// \param partial_augment By default the algorithm performs
2.361 + /// partial augment and relabel operations in the cost scaling
2.362 + /// phases. Set this parameter to \c false for using local push and
2.363 + /// relabel operations instead.
2.364 + ///
2.365 + /// \return \c true if a feasible flow can be found.
2.366 + bool run(bool partial_augment = true) {
2.367 + if (partial_augment) {
2.368 + return init() && startPartialAugment();
2.369 + } else {
2.370 + return init() && startPushRelabel();
2.371 + }
2.372 + }
2.373 +
2.374 + /// @}
2.375 +
2.376 + /// \name Query Functions
2.377 + /// The result of the algorithm can be obtained using these
2.378 + /// functions.\n
2.379 + /// \ref lemon::CostScaling::run() "run()" must be called before
2.380 + /// using them.
2.381 +
2.382 + /// @{
2.383 +
2.384 + /// \brief Return a const reference to the arc map storing the
2.385 + /// found flow.
2.386 + ///
2.387 + /// Return a const reference to the arc map storing the found flow.
2.388 + ///
2.389 + /// \pre \ref run() must be called before using this function.
2.390 + const FlowMap& flowMap() const {
2.391 + return *_flow;
2.392 + }
2.393 +
2.394 + /// \brief Return a const reference to the node map storing the
2.395 + /// found potentials (the dual solution).
2.396 + ///
2.397 + /// Return a const reference to the node map storing the found
2.398 + /// potentials (the dual solution).
2.399 + ///
2.400 + /// \pre \ref run() must be called before using this function.
2.401 + const PotentialMap& potentialMap() const {
2.402 + return *_potential;
2.403 + }
2.404 +
2.405 + /// \brief Return the flow on the given arc.
2.406 + ///
2.407 + /// Return the flow on the given arc.
2.408 + ///
2.409 + /// \pre \ref run() must be called before using this function.
2.410 + Capacity flow(const Arc& arc) const {
2.411 + return (*_flow)[arc];
2.412 + }
2.413 +
2.414 + /// \brief Return the potential of the given node.
2.415 + ///
2.416 + /// Return the potential of the given node.
2.417 + ///
2.418 + /// \pre \ref run() must be called before using this function.
2.419 + Cost potential(const Node& node) const {
2.420 + return (*_potential)[node];
2.421 + }
2.422 +
2.423 + /// \brief Return the total cost of the found flow.
2.424 + ///
2.425 + /// Return the total cost of the found flow. The complexity of the
2.426 + /// function is \f$ O(e) \f$.
2.427 + ///
2.428 + /// \pre \ref run() must be called before using this function.
2.429 + Cost totalCost() const {
2.430 + Cost c = 0;
2.431 + for (ArcIt e(_graph); e != INVALID; ++e)
2.432 + c += (*_flow)[e] * _orig_cost[e];
2.433 + return c;
2.434 + }
2.435 +
2.436 + /// @}
2.437 +
2.438 + private:
2.439 +
2.440 + /// Initialize the algorithm.
2.441 + bool init() {
2.442 + if (!_valid_supply) return false;
2.443 + // The scaling factor
2.444 + _alpha = 8;
2.445 +
2.446 + // Initialize flow and potential maps
2.447 + if (!_flow) {
2.448 + _flow = new FlowMap(_graph);
2.449 + _local_flow = true;
2.450 + }
2.451 + if (!_potential) {
2.452 + _potential = new PotentialMap(_graph);
2.453 + _local_potential = true;
2.454 + }
2.455 +
2.456 + _red_cost = new ReducedCostMap(_graph, _cost, *_potential);
2.457 + _res_graph = new ResDigraph(_graph, _capacity, *_flow);
2.458 +
2.459 + // Initialize the scaled cost map and the epsilon parameter
2.460 + Cost max_cost = 0;
2.461 + int node_num = countNodes(_graph);
2.462 + for (ArcIt e(_graph); e != INVALID; ++e) {
2.463 + _cost[e] = LCost(_orig_cost[e]) * node_num * _alpha;
2.464 + if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e];
2.465 + }
2.466 + _epsilon = max_cost * node_num;
2.467 +
2.468 + // Find a feasible flow using Circulation
2.469 + Circulation< Digraph, ConstMap<Arc, Capacity>, CapacityArcMap,
2.470 + SupplyMap >
2.471 + circulation( _graph, constMap<Arc>(Capacity(0)), _capacity,
2.472 + _supply );
2.473 + return circulation.flowMap(*_flow).run();
2.474 + }
2.475 +
2.476 + /// Execute the algorithm performing partial augmentation and
2.477 + /// relabel operations.
2.478 + bool startPartialAugment() {
2.479 + // Paramters for heuristics
2.480 +// const int BF_HEURISTIC_EPSILON_BOUND = 1000;
2.481 +// const int BF_HEURISTIC_BOUND_FACTOR = 3;
2.482 + // Maximum augment path length
2.483 + const int MAX_PATH_LENGTH = 4;
2.484 +
2.485 + // Variables
2.486 + typename Digraph::template NodeMap<Arc> pred_arc(_graph);
2.487 + typename Digraph::template NodeMap<bool> forward(_graph);
2.488 + typename Digraph::template NodeMap<OutArcIt> next_out(_graph);
2.489 + typename Digraph::template NodeMap<InArcIt> next_in(_graph);
2.490 + typename Digraph::template NodeMap<bool> next_dir(_graph);
2.491 + std::deque<Node> active_nodes;
2.492 + std::vector<Node> path_nodes;
2.493 +
2.494 +// int node_num = countNodes(_graph);
2.495 + for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
2.496 + 1 : _epsilon / _alpha )
2.497 + {
2.498 +/*
2.499 + // "Early Termination" heuristic: use Bellman-Ford algorithm
2.500 + // to check if the current flow is optimal
2.501 + if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
2.502 + typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
2.503 + ShiftCostMap shift_cost(_res_cost, 1);
2.504 + BellmanFord<ResDigraph, ShiftCostMap> bf(*_res_graph, shift_cost);
2.505 + bf.init(0);
2.506 + bool done = false;
2.507 + int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
2.508 + for (int i = 0; i < K && !done; ++i)
2.509 + done = bf.processNextWeakRound();
2.510 + if (done) break;
2.511 + }
2.512 +*/
2.513 + // Saturate arcs not satisfying the optimality condition
2.514 + Capacity delta;
2.515 + for (ArcIt e(_graph); e != INVALID; ++e) {
2.516 + if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
2.517 + delta = _capacity[e] - (*_flow)[e];
2.518 + _excess[_graph.source(e)] -= delta;
2.519 + _excess[_graph.target(e)] += delta;
2.520 + (*_flow)[e] = _capacity[e];
2.521 + }
2.522 + if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
2.523 + _excess[_graph.target(e)] -= (*_flow)[e];
2.524 + _excess[_graph.source(e)] += (*_flow)[e];
2.525 + (*_flow)[e] = 0;
2.526 + }
2.527 + }
2.528 +
2.529 + // Find active nodes (i.e. nodes with positive excess)
2.530 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.531 + if (_excess[n] > 0) active_nodes.push_back(n);
2.532 + }
2.533 +
2.534 + // Initialize the next arc maps
2.535 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.536 + next_out[n] = OutArcIt(_graph, n);
2.537 + next_in[n] = InArcIt(_graph, n);
2.538 + next_dir[n] = true;
2.539 + }
2.540 +
2.541 + // Perform partial augment and relabel operations
2.542 + while (active_nodes.size() > 0) {
2.543 + // Select an active node (FIFO selection)
2.544 + if (_excess[active_nodes[0]] <= 0) {
2.545 + active_nodes.pop_front();
2.546 + continue;
2.547 + }
2.548 + Node start = active_nodes[0];
2.549 + path_nodes.clear();
2.550 + path_nodes.push_back(start);
2.551 +
2.552 + // Find an augmenting path from the start node
2.553 + Node u, tip = start;
2.554 + LCost min_red_cost;
2.555 + while ( _excess[tip] >= 0 &&
2.556 + int(path_nodes.size()) <= MAX_PATH_LENGTH )
2.557 + {
2.558 + if (next_dir[tip]) {
2.559 + for (OutArcIt e = next_out[tip]; e != INVALID; ++e) {
2.560 + if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
2.561 + u = _graph.target(e);
2.562 + pred_arc[u] = e;
2.563 + forward[u] = true;
2.564 + next_out[tip] = e;
2.565 + tip = u;
2.566 + path_nodes.push_back(tip);
2.567 + goto next_step;
2.568 + }
2.569 + }
2.570 + next_dir[tip] = false;
2.571 + }
2.572 + for (InArcIt e = next_in[tip]; e != INVALID; ++e) {
2.573 + if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
2.574 + u = _graph.source(e);
2.575 + pred_arc[u] = e;
2.576 + forward[u] = false;
2.577 + next_in[tip] = e;
2.578 + tip = u;
2.579 + path_nodes.push_back(tip);
2.580 + goto next_step;
2.581 + }
2.582 + }
2.583 +
2.584 + // Relabel tip node
2.585 + min_red_cost = std::numeric_limits<LCost>::max() / 2;
2.586 + for (OutArcIt oe(_graph, tip); oe != INVALID; ++oe) {
2.587 + if ( _capacity[oe] - (*_flow)[oe] > 0 &&
2.588 + (*_red_cost)[oe] < min_red_cost )
2.589 + min_red_cost = (*_red_cost)[oe];
2.590 + }
2.591 + for (InArcIt ie(_graph, tip); ie != INVALID; ++ie) {
2.592 + if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
2.593 + min_red_cost = -(*_red_cost)[ie];
2.594 + }
2.595 + (*_potential)[tip] -= min_red_cost + _epsilon;
2.596 +
2.597 + // Reset the next arc maps
2.598 + next_out[tip] = OutArcIt(_graph, tip);
2.599 + next_in[tip] = InArcIt(_graph, tip);
2.600 + next_dir[tip] = true;
2.601 +
2.602 + // Step back
2.603 + if (tip != start) {
2.604 + path_nodes.pop_back();
2.605 + tip = path_nodes[path_nodes.size()-1];
2.606 + }
2.607 +
2.608 + next_step:
2.609 + continue;
2.610 + }
2.611 +
2.612 + // Augment along the found path (as much flow as possible)
2.613 + Capacity delta;
2.614 + for (int i = 1; i < int(path_nodes.size()); ++i) {
2.615 + u = path_nodes[i];
2.616 + delta = forward[u] ?
2.617 + _capacity[pred_arc[u]] - (*_flow)[pred_arc[u]] :
2.618 + (*_flow)[pred_arc[u]];
2.619 + delta = std::min(delta, _excess[path_nodes[i-1]]);
2.620 + (*_flow)[pred_arc[u]] += forward[u] ? delta : -delta;
2.621 + _excess[path_nodes[i-1]] -= delta;
2.622 + _excess[u] += delta;
2.623 + if (_excess[u] > 0 && _excess[u] <= delta) active_nodes.push_back(u);
2.624 + }
2.625 + }
2.626 + }
2.627 +
2.628 + // Compute node potentials for the original costs
2.629 + ResidualCostMap<CostMap> res_cost(_orig_cost);
2.630 + BellmanFord< ResDigraph, ResidualCostMap<CostMap> >
2.631 + bf(*_res_graph, res_cost);
2.632 + bf.init(0); bf.start();
2.633 + for (NodeIt n(_graph); n != INVALID; ++n)
2.634 + (*_potential)[n] = bf.dist(n);
2.635 +
2.636 + // Handle non-zero lower bounds
2.637 + if (_lower) {
2.638 + for (ArcIt e(_graph); e != INVALID; ++e)
2.639 + (*_flow)[e] += (*_lower)[e];
2.640 + }
2.641 + return true;
2.642 + }
2.643 +
2.644 + /// Execute the algorithm performing push and relabel operations.
2.645 + bool startPushRelabel() {
2.646 + // Paramters for heuristics
2.647 +// const int BF_HEURISTIC_EPSILON_BOUND = 1000;
2.648 +// const int BF_HEURISTIC_BOUND_FACTOR = 3;
2.649 +
2.650 + typename Digraph::template NodeMap<bool> hyper(_graph, false);
2.651 + typename Digraph::template NodeMap<Arc> pred_arc(_graph);
2.652 + typename Digraph::template NodeMap<bool> forward(_graph);
2.653 + typename Digraph::template NodeMap<OutArcIt> next_out(_graph);
2.654 + typename Digraph::template NodeMap<InArcIt> next_in(_graph);
2.655 + typename Digraph::template NodeMap<bool> next_dir(_graph);
2.656 + std::deque<Node> active_nodes;
2.657 +
2.658 +// int node_num = countNodes(_graph);
2.659 + for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
2.660 + 1 : _epsilon / _alpha )
2.661 + {
2.662 +/*
2.663 + // "Early Termination" heuristic: use Bellman-Ford algorithm
2.664 + // to check if the current flow is optimal
2.665 + if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
2.666 + typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
2.667 + ShiftCostMap shift_cost(_res_cost, 1);
2.668 + BellmanFord<ResDigraph, ShiftCostMap> bf(*_res_graph, shift_cost);
2.669 + bf.init(0);
2.670 + bool done = false;
2.671 + int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
2.672 + for (int i = 0; i < K && !done; ++i)
2.673 + done = bf.processNextWeakRound();
2.674 + if (done) break;
2.675 + }
2.676 +*/
2.677 +
2.678 + // Saturate arcs not satisfying the optimality condition
2.679 + Capacity delta;
2.680 + for (ArcIt e(_graph); e != INVALID; ++e) {
2.681 + if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
2.682 + delta = _capacity[e] - (*_flow)[e];
2.683 + _excess[_graph.source(e)] -= delta;
2.684 + _excess[_graph.target(e)] += delta;
2.685 + (*_flow)[e] = _capacity[e];
2.686 + }
2.687 + if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
2.688 + _excess[_graph.target(e)] -= (*_flow)[e];
2.689 + _excess[_graph.source(e)] += (*_flow)[e];
2.690 + (*_flow)[e] = 0;
2.691 + }
2.692 + }
2.693 +
2.694 + // Find active nodes (i.e. nodes with positive excess)
2.695 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.696 + if (_excess[n] > 0) active_nodes.push_back(n);
2.697 + }
2.698 +
2.699 + // Initialize the next arc maps
2.700 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.701 + next_out[n] = OutArcIt(_graph, n);
2.702 + next_in[n] = InArcIt(_graph, n);
2.703 + next_dir[n] = true;
2.704 + }
2.705 +
2.706 + // Perform push and relabel operations
2.707 + while (active_nodes.size() > 0) {
2.708 + // Select an active node (FIFO selection)
2.709 + Node n = active_nodes[0], t;
2.710 + bool relabel_enabled = true;
2.711 +
2.712 + // Perform push operations if there are admissible arcs
2.713 + if (_excess[n] > 0 && next_dir[n]) {
2.714 + OutArcIt e = next_out[n];
2.715 + for ( ; e != INVALID; ++e) {
2.716 + if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
2.717 + delta = std::min(_capacity[e] - (*_flow)[e], _excess[n]);
2.718 + t = _graph.target(e);
2.719 +
2.720 + // Push-look-ahead heuristic
2.721 + Capacity ahead = -_excess[t];
2.722 + for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) {
2.723 + if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
2.724 + ahead += _capacity[oe] - (*_flow)[oe];
2.725 + }
2.726 + for (InArcIt ie(_graph, t); ie != INVALID; ++ie) {
2.727 + if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
2.728 + ahead += (*_flow)[ie];
2.729 + }
2.730 + if (ahead < 0) ahead = 0;
2.731 +
2.732 + // Push flow along the arc
2.733 + if (ahead < delta) {
2.734 + (*_flow)[e] += ahead;
2.735 + _excess[n] -= ahead;
2.736 + _excess[t] += ahead;
2.737 + active_nodes.push_front(t);
2.738 + hyper[t] = true;
2.739 + relabel_enabled = false;
2.740 + break;
2.741 + } else {
2.742 + (*_flow)[e] += delta;
2.743 + _excess[n] -= delta;
2.744 + _excess[t] += delta;
2.745 + if (_excess[t] > 0 && _excess[t] <= delta)
2.746 + active_nodes.push_back(t);
2.747 + }
2.748 +
2.749 + if (_excess[n] == 0) break;
2.750 + }
2.751 + }
2.752 + if (e != INVALID) {
2.753 + next_out[n] = e;
2.754 + } else {
2.755 + next_dir[n] = false;
2.756 + }
2.757 + }
2.758 +
2.759 + if (_excess[n] > 0 && !next_dir[n]) {
2.760 + InArcIt e = next_in[n];
2.761 + for ( ; e != INVALID; ++e) {
2.762 + if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
2.763 + delta = std::min((*_flow)[e], _excess[n]);
2.764 + t = _graph.source(e);
2.765 +
2.766 + // Push-look-ahead heuristic
2.767 + Capacity ahead = -_excess[t];
2.768 + for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) {
2.769 + if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
2.770 + ahead += _capacity[oe] - (*_flow)[oe];
2.771 + }
2.772 + for (InArcIt ie(_graph, t); ie != INVALID; ++ie) {
2.773 + if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
2.774 + ahead += (*_flow)[ie];
2.775 + }
2.776 + if (ahead < 0) ahead = 0;
2.777 +
2.778 + // Push flow along the arc
2.779 + if (ahead < delta) {
2.780 + (*_flow)[e] -= ahead;
2.781 + _excess[n] -= ahead;
2.782 + _excess[t] += ahead;
2.783 + active_nodes.push_front(t);
2.784 + hyper[t] = true;
2.785 + relabel_enabled = false;
2.786 + break;
2.787 + } else {
2.788 + (*_flow)[e] -= delta;
2.789 + _excess[n] -= delta;
2.790 + _excess[t] += delta;
2.791 + if (_excess[t] > 0 && _excess[t] <= delta)
2.792 + active_nodes.push_back(t);
2.793 + }
2.794 +
2.795 + if (_excess[n] == 0) break;
2.796 + }
2.797 + }
2.798 + next_in[n] = e;
2.799 + }
2.800 +
2.801 + // Relabel the node if it is still active (or hyper)
2.802 + if (relabel_enabled && (_excess[n] > 0 || hyper[n])) {
2.803 + LCost min_red_cost = std::numeric_limits<LCost>::max() / 2;
2.804 + for (OutArcIt oe(_graph, n); oe != INVALID; ++oe) {
2.805 + if ( _capacity[oe] - (*_flow)[oe] > 0 &&
2.806 + (*_red_cost)[oe] < min_red_cost )
2.807 + min_red_cost = (*_red_cost)[oe];
2.808 + }
2.809 + for (InArcIt ie(_graph, n); ie != INVALID; ++ie) {
2.810 + if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
2.811 + min_red_cost = -(*_red_cost)[ie];
2.812 + }
2.813 + (*_potential)[n] -= min_red_cost + _epsilon;
2.814 + hyper[n] = false;
2.815 +
2.816 + // Reset the next arc maps
2.817 + next_out[n] = OutArcIt(_graph, n);
2.818 + next_in[n] = InArcIt(_graph, n);
2.819 + next_dir[n] = true;
2.820 + }
2.821 +
2.822 + // Remove nodes that are not active nor hyper
2.823 + while ( active_nodes.size() > 0 &&
2.824 + _excess[active_nodes[0]] <= 0 &&
2.825 + !hyper[active_nodes[0]] ) {
2.826 + active_nodes.pop_front();
2.827 + }
2.828 + }
2.829 + }
2.830 +
2.831 + // Compute node potentials for the original costs
2.832 + ResidualCostMap<CostMap> res_cost(_orig_cost);
2.833 + BellmanFord< ResDigraph, ResidualCostMap<CostMap> >
2.834 + bf(*_res_graph, res_cost);
2.835 + bf.init(0); bf.start();
2.836 + for (NodeIt n(_graph); n != INVALID; ++n)
2.837 + (*_potential)[n] = bf.dist(n);
2.838 +
2.839 + // Handle non-zero lower bounds
2.840 + if (_lower) {
2.841 + for (ArcIt e(_graph); e != INVALID; ++e)
2.842 + (*_flow)[e] += (*_lower)[e];
2.843 + }
2.844 + return true;
2.845 + }
2.846 +
2.847 + }; //class CostScaling
2.848 +
2.849 + ///@}
2.850 +
2.851 +} //namespace lemon
2.852 +
2.853 +#endif //LEMON_COST_SCALING_H