1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/cost_scaling.h Sun Aug 11 15:28:12 2013 +0200
1.3 @@ -0,0 +1,1316 @@
1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library.
1.7 + *
1.8 + * Copyright (C) 2003-2010
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/static_graph.h>
1.37 +#include <lemon/circulation.h>
1.38 +#include <lemon/bellman_ford.h>
1.39 +
1.40 +namespace lemon {
1.41 +
1.42 + /// \brief Default traits class of CostScaling algorithm.
1.43 + ///
1.44 + /// Default traits class of CostScaling algorithm.
1.45 + /// \tparam GR Digraph type.
1.46 + /// \tparam V The number type used for flow amounts, capacity bounds
1.47 + /// and supply values. By default it is \c int.
1.48 + /// \tparam C The number type used for costs and potentials.
1.49 + /// By default it is the same as \c V.
1.50 +#ifdef DOXYGEN
1.51 + template <typename GR, typename V = int, typename C = V>
1.52 +#else
1.53 + template < typename GR, typename V = int, typename C = V,
1.54 + bool integer = std::numeric_limits<C>::is_integer >
1.55 +#endif
1.56 + struct CostScalingDefaultTraits
1.57 + {
1.58 + /// The type of the digraph
1.59 + typedef GR Digraph;
1.60 + /// The type of the flow amounts, capacity bounds and supply values
1.61 + typedef V Value;
1.62 + /// The type of the arc costs
1.63 + typedef C Cost;
1.64 +
1.65 + /// \brief The large cost type used for internal computations
1.66 + ///
1.67 + /// The large cost type used for internal computations.
1.68 + /// It is \c long \c long if the \c Cost type is integer,
1.69 + /// otherwise it is \c double.
1.70 + /// \c Cost must be convertible to \c LargeCost.
1.71 + typedef double LargeCost;
1.72 + };
1.73 +
1.74 + // Default traits class for integer cost types
1.75 + template <typename GR, typename V, typename C>
1.76 + struct CostScalingDefaultTraits<GR, V, C, true>
1.77 + {
1.78 + typedef GR Digraph;
1.79 + typedef V Value;
1.80 + typedef C Cost;
1.81 +#ifdef LEMON_HAVE_LONG_LONG
1.82 + typedef long long LargeCost;
1.83 +#else
1.84 + typedef long LargeCost;
1.85 +#endif
1.86 + };
1.87 +
1.88 +
1.89 + /// \addtogroup min_cost_flow_algs
1.90 + /// @{
1.91 +
1.92 + /// \brief Implementation of the Cost Scaling algorithm for
1.93 + /// finding a \ref min_cost_flow "minimum cost flow".
1.94 + ///
1.95 + /// \ref CostScaling implements a cost scaling algorithm that performs
1.96 + /// push/augment and relabel operations for finding a \ref min_cost_flow
1.97 + /// "minimum cost flow" \ref amo93networkflows, \ref goldberg90approximation,
1.98 + /// \ref goldberg97efficient, \ref bunnagel98efficient.
1.99 + /// It is a highly efficient primal-dual solution method, which
1.100 + /// can be viewed as the generalization of the \ref Preflow
1.101 + /// "preflow push-relabel" algorithm for the maximum flow problem.
1.102 + ///
1.103 + /// Most of the parameters of the problem (except for the digraph)
1.104 + /// can be given using separate functions, and the algorithm can be
1.105 + /// executed using the \ref run() function. If some parameters are not
1.106 + /// specified, then default values will be used.
1.107 + ///
1.108 + /// \tparam GR The digraph type the algorithm runs on.
1.109 + /// \tparam V The number type used for flow amounts, capacity bounds
1.110 + /// and supply values in the algorithm. By default, it is \c int.
1.111 + /// \tparam C The number type used for costs and potentials in the
1.112 + /// algorithm. By default, it is the same as \c V.
1.113 + /// \tparam TR The traits class that defines various types used by the
1.114 + /// algorithm. By default, it is \ref CostScalingDefaultTraits
1.115 + /// "CostScalingDefaultTraits<GR, V, C>".
1.116 + /// In most cases, this parameter should not be set directly,
1.117 + /// consider to use the named template parameters instead.
1.118 + ///
1.119 + /// \warning Both number types must be signed and all input data must
1.120 + /// be integer.
1.121 + /// \warning This algorithm does not support negative costs for such
1.122 + /// arcs that have infinite upper bound.
1.123 + ///
1.124 + /// \note %CostScaling provides three different internal methods,
1.125 + /// from which the most efficient one is used by default.
1.126 + /// For more information, see \ref Method.
1.127 +#ifdef DOXYGEN
1.128 + template <typename GR, typename V, typename C, typename TR>
1.129 +#else
1.130 + template < typename GR, typename V = int, typename C = V,
1.131 + typename TR = CostScalingDefaultTraits<GR, V, C> >
1.132 +#endif
1.133 + class CostScaling
1.134 + {
1.135 + public:
1.136 +
1.137 + /// The type of the digraph
1.138 + typedef typename TR::Digraph Digraph;
1.139 + /// The type of the flow amounts, capacity bounds and supply values
1.140 + typedef typename TR::Value Value;
1.141 + /// The type of the arc costs
1.142 + typedef typename TR::Cost Cost;
1.143 +
1.144 + /// \brief The large cost type
1.145 + ///
1.146 + /// The large cost type used for internal computations.
1.147 + /// By default, it is \c long \c long if the \c Cost type is integer,
1.148 + /// otherwise it is \c double.
1.149 + typedef typename TR::LargeCost LargeCost;
1.150 +
1.151 + /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
1.152 + typedef TR Traits;
1.153 +
1.154 + public:
1.155 +
1.156 + /// \brief Problem type constants for the \c run() function.
1.157 + ///
1.158 + /// Enum type containing the problem type constants that can be
1.159 + /// returned by the \ref run() function of the algorithm.
1.160 + enum ProblemType {
1.161 + /// The problem has no feasible solution (flow).
1.162 + INFEASIBLE,
1.163 + /// The problem has optimal solution (i.e. it is feasible and
1.164 + /// bounded), and the algorithm has found optimal flow and node
1.165 + /// potentials (primal and dual solutions).
1.166 + OPTIMAL,
1.167 + /// The digraph contains an arc of negative cost and infinite
1.168 + /// upper bound. It means that the objective function is unbounded
1.169 + /// on that arc, however, note that it could actually be bounded
1.170 + /// over the feasible flows, but this algroithm cannot handle
1.171 + /// these cases.
1.172 + UNBOUNDED
1.173 + };
1.174 +
1.175 + /// \brief Constants for selecting the internal method.
1.176 + ///
1.177 + /// Enum type containing constants for selecting the internal method
1.178 + /// for the \ref run() function.
1.179 + ///
1.180 + /// \ref CostScaling provides three internal methods that differ mainly
1.181 + /// in their base operations, which are used in conjunction with the
1.182 + /// relabel operation.
1.183 + /// By default, the so called \ref PARTIAL_AUGMENT
1.184 + /// "Partial Augment-Relabel" method is used, which proved to be
1.185 + /// the most efficient and the most robust on various test inputs.
1.186 + /// However, the other methods can be selected using the \ref run()
1.187 + /// function with the proper parameter.
1.188 + enum Method {
1.189 + /// Local push operations are used, i.e. flow is moved only on one
1.190 + /// admissible arc at once.
1.191 + PUSH,
1.192 + /// Augment operations are used, i.e. flow is moved on admissible
1.193 + /// paths from a node with excess to a node with deficit.
1.194 + AUGMENT,
1.195 + /// Partial augment operations are used, i.e. flow is moved on
1.196 + /// admissible paths started from a node with excess, but the
1.197 + /// lengths of these paths are limited. This method can be viewed
1.198 + /// as a combined version of the previous two operations.
1.199 + PARTIAL_AUGMENT
1.200 + };
1.201 +
1.202 + private:
1.203 +
1.204 + TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.205 +
1.206 + typedef std::vector<int> IntVector;
1.207 + typedef std::vector<Value> ValueVector;
1.208 + typedef std::vector<Cost> CostVector;
1.209 + typedef std::vector<LargeCost> LargeCostVector;
1.210 + typedef std::vector<char> BoolVector;
1.211 + // Note: vector<char> is used instead of vector<bool> for efficiency reasons
1.212 +
1.213 + private:
1.214 +
1.215 + template <typename KT, typename VT>
1.216 + class StaticVectorMap {
1.217 + public:
1.218 + typedef KT Key;
1.219 + typedef VT Value;
1.220 +
1.221 + StaticVectorMap(std::vector<Value>& v) : _v(v) {}
1.222 +
1.223 + const Value& operator[](const Key& key) const {
1.224 + return _v[StaticDigraph::id(key)];
1.225 + }
1.226 +
1.227 + Value& operator[](const Key& key) {
1.228 + return _v[StaticDigraph::id(key)];
1.229 + }
1.230 +
1.231 + void set(const Key& key, const Value& val) {
1.232 + _v[StaticDigraph::id(key)] = val;
1.233 + }
1.234 +
1.235 + private:
1.236 + std::vector<Value>& _v;
1.237 + };
1.238 +
1.239 + typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
1.240 + typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
1.241 +
1.242 + private:
1.243 +
1.244 + // Data related to the underlying digraph
1.245 + const GR &_graph;
1.246 + int _node_num;
1.247 + int _arc_num;
1.248 + int _res_node_num;
1.249 + int _res_arc_num;
1.250 + int _root;
1.251 +
1.252 + // Parameters of the problem
1.253 + bool _have_lower;
1.254 + Value _sum_supply;
1.255 + int _sup_node_num;
1.256 +
1.257 + // Data structures for storing the digraph
1.258 + IntNodeMap _node_id;
1.259 + IntArcMap _arc_idf;
1.260 + IntArcMap _arc_idb;
1.261 + IntVector _first_out;
1.262 + BoolVector _forward;
1.263 + IntVector _source;
1.264 + IntVector _target;
1.265 + IntVector _reverse;
1.266 +
1.267 + // Node and arc data
1.268 + ValueVector _lower;
1.269 + ValueVector _upper;
1.270 + CostVector _scost;
1.271 + ValueVector _supply;
1.272 +
1.273 + ValueVector _res_cap;
1.274 + LargeCostVector _cost;
1.275 + LargeCostVector _pi;
1.276 + ValueVector _excess;
1.277 + IntVector _next_out;
1.278 + std::deque<int> _active_nodes;
1.279 +
1.280 + // Data for scaling
1.281 + LargeCost _epsilon;
1.282 + int _alpha;
1.283 +
1.284 + IntVector _buckets;
1.285 + IntVector _bucket_next;
1.286 + IntVector _bucket_prev;
1.287 + IntVector _rank;
1.288 + int _max_rank;
1.289 +
1.290 + // Data for a StaticDigraph structure
1.291 + typedef std::pair<int, int> IntPair;
1.292 + StaticDigraph _sgr;
1.293 + std::vector<IntPair> _arc_vec;
1.294 + std::vector<LargeCost> _cost_vec;
1.295 + LargeCostArcMap _cost_map;
1.296 + LargeCostNodeMap _pi_map;
1.297 +
1.298 + public:
1.299 +
1.300 + /// \brief Constant for infinite upper bounds (capacities).
1.301 + ///
1.302 + /// Constant for infinite upper bounds (capacities).
1.303 + /// It is \c std::numeric_limits<Value>::infinity() if available,
1.304 + /// \c std::numeric_limits<Value>::max() otherwise.
1.305 + const Value INF;
1.306 +
1.307 + public:
1.308 +
1.309 + /// \name Named Template Parameters
1.310 + /// @{
1.311 +
1.312 + template <typename T>
1.313 + struct SetLargeCostTraits : public Traits {
1.314 + typedef T LargeCost;
1.315 + };
1.316 +
1.317 + /// \brief \ref named-templ-param "Named parameter" for setting
1.318 + /// \c LargeCost type.
1.319 + ///
1.320 + /// \ref named-templ-param "Named parameter" for setting \c LargeCost
1.321 + /// type, which is used for internal computations in the algorithm.
1.322 + /// \c Cost must be convertible to \c LargeCost.
1.323 + template <typename T>
1.324 + struct SetLargeCost
1.325 + : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {
1.326 + typedef CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;
1.327 + };
1.328 +
1.329 + /// @}
1.330 +
1.331 + protected:
1.332 +
1.333 + CostScaling() {}
1.334 +
1.335 + public:
1.336 +
1.337 + /// \brief Constructor.
1.338 + ///
1.339 + /// The constructor of the class.
1.340 + ///
1.341 + /// \param graph The digraph the algorithm runs on.
1.342 + CostScaling(const GR& graph) :
1.343 + _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
1.344 + _cost_map(_cost_vec), _pi_map(_pi),
1.345 + INF(std::numeric_limits<Value>::has_infinity ?
1.346 + std::numeric_limits<Value>::infinity() :
1.347 + std::numeric_limits<Value>::max())
1.348 + {
1.349 + // Check the number types
1.350 + LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
1.351 + "The flow type of CostScaling must be signed");
1.352 + LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
1.353 + "The cost type of CostScaling must be signed");
1.354 +
1.355 + // Reset data structures
1.356 + reset();
1.357 + }
1.358 +
1.359 + /// \name Parameters
1.360 + /// The parameters of the algorithm can be specified using these
1.361 + /// functions.
1.362 +
1.363 + /// @{
1.364 +
1.365 + /// \brief Set the lower bounds on the arcs.
1.366 + ///
1.367 + /// This function sets the lower bounds on the arcs.
1.368 + /// If it is not used before calling \ref run(), the lower bounds
1.369 + /// will be set to zero on all arcs.
1.370 + ///
1.371 + /// \param map An arc map storing the lower bounds.
1.372 + /// Its \c Value type must be convertible to the \c Value type
1.373 + /// of the algorithm.
1.374 + ///
1.375 + /// \return <tt>(*this)</tt>
1.376 + template <typename LowerMap>
1.377 + CostScaling& lowerMap(const LowerMap& map) {
1.378 + _have_lower = true;
1.379 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.380 + _lower[_arc_idf[a]] = map[a];
1.381 + _lower[_arc_idb[a]] = map[a];
1.382 + }
1.383 + return *this;
1.384 + }
1.385 +
1.386 + /// \brief Set the upper bounds (capacities) on the arcs.
1.387 + ///
1.388 + /// This function sets the upper bounds (capacities) on the arcs.
1.389 + /// If it is not used before calling \ref run(), the upper bounds
1.390 + /// will be set to \ref INF on all arcs (i.e. the flow value will be
1.391 + /// unbounded from above).
1.392 + ///
1.393 + /// \param map An arc map storing the upper bounds.
1.394 + /// Its \c Value type must be convertible to the \c Value type
1.395 + /// of the algorithm.
1.396 + ///
1.397 + /// \return <tt>(*this)</tt>
1.398 + template<typename UpperMap>
1.399 + CostScaling& upperMap(const UpperMap& map) {
1.400 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.401 + _upper[_arc_idf[a]] = map[a];
1.402 + }
1.403 + return *this;
1.404 + }
1.405 +
1.406 + /// \brief Set the costs of the arcs.
1.407 + ///
1.408 + /// This function sets the costs of the arcs.
1.409 + /// If it is not used before calling \ref run(), the costs
1.410 + /// will be set to \c 1 on all arcs.
1.411 + ///
1.412 + /// \param map An arc map storing the costs.
1.413 + /// Its \c Value type must be convertible to the \c Cost type
1.414 + /// of the algorithm.
1.415 + ///
1.416 + /// \return <tt>(*this)</tt>
1.417 + template<typename CostMap>
1.418 + CostScaling& costMap(const CostMap& map) {
1.419 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.420 + _scost[_arc_idf[a]] = map[a];
1.421 + _scost[_arc_idb[a]] = -map[a];
1.422 + }
1.423 + return *this;
1.424 + }
1.425 +
1.426 + /// \brief Set the supply values of the nodes.
1.427 + ///
1.428 + /// This function sets the supply values of the nodes.
1.429 + /// If neither this function nor \ref stSupply() is used before
1.430 + /// calling \ref run(), the supply of each node will be set to zero.
1.431 + ///
1.432 + /// \param map A node map storing the supply values.
1.433 + /// Its \c Value type must be convertible to the \c Value type
1.434 + /// of the algorithm.
1.435 + ///
1.436 + /// \return <tt>(*this)</tt>
1.437 + template<typename SupplyMap>
1.438 + CostScaling& supplyMap(const SupplyMap& map) {
1.439 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.440 + _supply[_node_id[n]] = map[n];
1.441 + }
1.442 + return *this;
1.443 + }
1.444 +
1.445 + /// \brief Set single source and target nodes and a supply value.
1.446 + ///
1.447 + /// This function sets a single source node and a single target node
1.448 + /// and the required flow value.
1.449 + /// If neither this function nor \ref supplyMap() is used before
1.450 + /// calling \ref run(), the supply of each node will be set to zero.
1.451 + ///
1.452 + /// Using this function has the same effect as using \ref supplyMap()
1.453 + /// with such a map in which \c k is assigned to \c s, \c -k is
1.454 + /// assigned to \c t and all other nodes have zero supply value.
1.455 + ///
1.456 + /// \param s The source node.
1.457 + /// \param t The target node.
1.458 + /// \param k The required amount of flow from node \c s to node \c t
1.459 + /// (i.e. the supply of \c s and the demand of \c t).
1.460 + ///
1.461 + /// \return <tt>(*this)</tt>
1.462 + CostScaling& stSupply(const Node& s, const Node& t, Value k) {
1.463 + for (int i = 0; i != _res_node_num; ++i) {
1.464 + _supply[i] = 0;
1.465 + }
1.466 + _supply[_node_id[s]] = k;
1.467 + _supply[_node_id[t]] = -k;
1.468 + return *this;
1.469 + }
1.470 +
1.471 + /// @}
1.472 +
1.473 + /// \name Execution control
1.474 + /// The algorithm can be executed using \ref run().
1.475 +
1.476 + /// @{
1.477 +
1.478 + /// \brief Run the algorithm.
1.479 + ///
1.480 + /// This function runs the algorithm.
1.481 + /// The paramters can be specified using functions \ref lowerMap(),
1.482 + /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
1.483 + /// For example,
1.484 + /// \code
1.485 + /// CostScaling<ListDigraph> cs(graph);
1.486 + /// cs.lowerMap(lower).upperMap(upper).costMap(cost)
1.487 + /// .supplyMap(sup).run();
1.488 + /// \endcode
1.489 + ///
1.490 + /// This function can be called more than once. All the given parameters
1.491 + /// are kept for the next call, unless \ref resetParams() or \ref reset()
1.492 + /// is used, thus only the modified parameters have to be set again.
1.493 + /// If the underlying digraph was also modified after the construction
1.494 + /// of the class (or the last \ref reset() call), then the \ref reset()
1.495 + /// function must be called.
1.496 + ///
1.497 + /// \param method The internal method that will be used in the
1.498 + /// algorithm. For more information, see \ref Method.
1.499 + /// \param factor The cost scaling factor. It must be larger than one.
1.500 + ///
1.501 + /// \return \c INFEASIBLE if no feasible flow exists,
1.502 + /// \n \c OPTIMAL if the problem has optimal solution
1.503 + /// (i.e. it is feasible and bounded), and the algorithm has found
1.504 + /// optimal flow and node potentials (primal and dual solutions),
1.505 + /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
1.506 + /// and infinite upper bound. It means that the objective function
1.507 + /// is unbounded on that arc, however, note that it could actually be
1.508 + /// bounded over the feasible flows, but this algroithm cannot handle
1.509 + /// these cases.
1.510 + ///
1.511 + /// \see ProblemType, Method
1.512 + /// \see resetParams(), reset()
1.513 + ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) {
1.514 + _alpha = factor;
1.515 + ProblemType pt = init();
1.516 + if (pt != OPTIMAL) return pt;
1.517 + start(method);
1.518 + return OPTIMAL;
1.519 + }
1.520 +
1.521 + /// \brief Reset all the parameters that have been given before.
1.522 + ///
1.523 + /// This function resets all the paramaters that have been given
1.524 + /// before using functions \ref lowerMap(), \ref upperMap(),
1.525 + /// \ref costMap(), \ref supplyMap(), \ref stSupply().
1.526 + ///
1.527 + /// It is useful for multiple \ref run() calls. Basically, all the given
1.528 + /// parameters are kept for the next \ref run() call, unless
1.529 + /// \ref resetParams() or \ref reset() is used.
1.530 + /// If the underlying digraph was also modified after the construction
1.531 + /// of the class or the last \ref reset() call, then the \ref reset()
1.532 + /// function must be used, otherwise \ref resetParams() is sufficient.
1.533 + ///
1.534 + /// For example,
1.535 + /// \code
1.536 + /// CostScaling<ListDigraph> cs(graph);
1.537 + ///
1.538 + /// // First run
1.539 + /// cs.lowerMap(lower).upperMap(upper).costMap(cost)
1.540 + /// .supplyMap(sup).run();
1.541 + ///
1.542 + /// // Run again with modified cost map (resetParams() is not called,
1.543 + /// // so only the cost map have to be set again)
1.544 + /// cost[e] += 100;
1.545 + /// cs.costMap(cost).run();
1.546 + ///
1.547 + /// // Run again from scratch using resetParams()
1.548 + /// // (the lower bounds will be set to zero on all arcs)
1.549 + /// cs.resetParams();
1.550 + /// cs.upperMap(capacity).costMap(cost)
1.551 + /// .supplyMap(sup).run();
1.552 + /// \endcode
1.553 + ///
1.554 + /// \return <tt>(*this)</tt>
1.555 + ///
1.556 + /// \see reset(), run()
1.557 + CostScaling& resetParams() {
1.558 + for (int i = 0; i != _res_node_num; ++i) {
1.559 + _supply[i] = 0;
1.560 + }
1.561 + int limit = _first_out[_root];
1.562 + for (int j = 0; j != limit; ++j) {
1.563 + _lower[j] = 0;
1.564 + _upper[j] = INF;
1.565 + _scost[j] = _forward[j] ? 1 : -1;
1.566 + }
1.567 + for (int j = limit; j != _res_arc_num; ++j) {
1.568 + _lower[j] = 0;
1.569 + _upper[j] = INF;
1.570 + _scost[j] = 0;
1.571 + _scost[_reverse[j]] = 0;
1.572 + }
1.573 + _have_lower = false;
1.574 + return *this;
1.575 + }
1.576 +
1.577 + /// \brief Reset all the parameters that have been given before.
1.578 + ///
1.579 + /// This function resets all the paramaters that have been given
1.580 + /// before using functions \ref lowerMap(), \ref upperMap(),
1.581 + /// \ref costMap(), \ref supplyMap(), \ref stSupply().
1.582 + ///
1.583 + /// It is useful for multiple run() calls. If this function is not
1.584 + /// used, all the parameters given before are kept for the next
1.585 + /// \ref run() call.
1.586 + /// However, the underlying digraph must not be modified after this
1.587 + /// class have been constructed, since it copies and extends the graph.
1.588 + /// \return <tt>(*this)</tt>
1.589 + CostScaling& reset() {
1.590 + // Resize vectors
1.591 + _node_num = countNodes(_graph);
1.592 + _arc_num = countArcs(_graph);
1.593 + _res_node_num = _node_num + 1;
1.594 + _res_arc_num = 2 * (_arc_num + _node_num);
1.595 + _root = _node_num;
1.596 +
1.597 + _first_out.resize(_res_node_num + 1);
1.598 + _forward.resize(_res_arc_num);
1.599 + _source.resize(_res_arc_num);
1.600 + _target.resize(_res_arc_num);
1.601 + _reverse.resize(_res_arc_num);
1.602 +
1.603 + _lower.resize(_res_arc_num);
1.604 + _upper.resize(_res_arc_num);
1.605 + _scost.resize(_res_arc_num);
1.606 + _supply.resize(_res_node_num);
1.607 +
1.608 + _res_cap.resize(_res_arc_num);
1.609 + _cost.resize(_res_arc_num);
1.610 + _pi.resize(_res_node_num);
1.611 + _excess.resize(_res_node_num);
1.612 + _next_out.resize(_res_node_num);
1.613 +
1.614 + _arc_vec.reserve(_res_arc_num);
1.615 + _cost_vec.reserve(_res_arc_num);
1.616 +
1.617 + // Copy the graph
1.618 + int i = 0, j = 0, k = 2 * _arc_num + _node_num;
1.619 + for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.620 + _node_id[n] = i;
1.621 + }
1.622 + i = 0;
1.623 + for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.624 + _first_out[i] = j;
1.625 + for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
1.626 + _arc_idf[a] = j;
1.627 + _forward[j] = true;
1.628 + _source[j] = i;
1.629 + _target[j] = _node_id[_graph.runningNode(a)];
1.630 + }
1.631 + for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
1.632 + _arc_idb[a] = j;
1.633 + _forward[j] = false;
1.634 + _source[j] = i;
1.635 + _target[j] = _node_id[_graph.runningNode(a)];
1.636 + }
1.637 + _forward[j] = false;
1.638 + _source[j] = i;
1.639 + _target[j] = _root;
1.640 + _reverse[j] = k;
1.641 + _forward[k] = true;
1.642 + _source[k] = _root;
1.643 + _target[k] = i;
1.644 + _reverse[k] = j;
1.645 + ++j; ++k;
1.646 + }
1.647 + _first_out[i] = j;
1.648 + _first_out[_res_node_num] = k;
1.649 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.650 + int fi = _arc_idf[a];
1.651 + int bi = _arc_idb[a];
1.652 + _reverse[fi] = bi;
1.653 + _reverse[bi] = fi;
1.654 + }
1.655 +
1.656 + // Reset parameters
1.657 + resetParams();
1.658 + return *this;
1.659 + }
1.660 +
1.661 + /// @}
1.662 +
1.663 + /// \name Query Functions
1.664 + /// The results of the algorithm can be obtained using these
1.665 + /// functions.\n
1.666 + /// The \ref run() function must be called before using them.
1.667 +
1.668 + /// @{
1.669 +
1.670 + /// \brief Return the total cost of the found flow.
1.671 + ///
1.672 + /// This function returns the total cost of the found flow.
1.673 + /// Its complexity is O(e).
1.674 + ///
1.675 + /// \note The return type of the function can be specified as a
1.676 + /// template parameter. For example,
1.677 + /// \code
1.678 + /// cs.totalCost<double>();
1.679 + /// \endcode
1.680 + /// It is useful if the total cost cannot be stored in the \c Cost
1.681 + /// type of the algorithm, which is the default return type of the
1.682 + /// function.
1.683 + ///
1.684 + /// \pre \ref run() must be called before using this function.
1.685 + template <typename Number>
1.686 + Number totalCost() const {
1.687 + Number c = 0;
1.688 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.689 + int i = _arc_idb[a];
1.690 + c += static_cast<Number>(_res_cap[i]) *
1.691 + (-static_cast<Number>(_scost[i]));
1.692 + }
1.693 + return c;
1.694 + }
1.695 +
1.696 +#ifndef DOXYGEN
1.697 + Cost totalCost() const {
1.698 + return totalCost<Cost>();
1.699 + }
1.700 +#endif
1.701 +
1.702 + /// \brief Return the flow on the given arc.
1.703 + ///
1.704 + /// This function returns the flow on the given arc.
1.705 + ///
1.706 + /// \pre \ref run() must be called before using this function.
1.707 + Value flow(const Arc& a) const {
1.708 + return _res_cap[_arc_idb[a]];
1.709 + }
1.710 +
1.711 + /// \brief Return the flow map (the primal solution).
1.712 + ///
1.713 + /// This function copies the flow value on each arc into the given
1.714 + /// map. The \c Value type of the algorithm must be convertible to
1.715 + /// the \c Value type of the map.
1.716 + ///
1.717 + /// \pre \ref run() must be called before using this function.
1.718 + template <typename FlowMap>
1.719 + void flowMap(FlowMap &map) const {
1.720 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.721 + map.set(a, _res_cap[_arc_idb[a]]);
1.722 + }
1.723 + }
1.724 +
1.725 + /// \brief Return the potential (dual value) of the given node.
1.726 + ///
1.727 + /// This function returns the potential (dual value) of the
1.728 + /// given node.
1.729 + ///
1.730 + /// \pre \ref run() must be called before using this function.
1.731 + Cost potential(const Node& n) const {
1.732 + return static_cast<Cost>(_pi[_node_id[n]]);
1.733 + }
1.734 +
1.735 + /// \brief Return the potential map (the dual solution).
1.736 + ///
1.737 + /// This function copies the potential (dual value) of each node
1.738 + /// into the given map.
1.739 + /// The \c Cost type of the algorithm must be convertible to the
1.740 + /// \c Value type of the map.
1.741 + ///
1.742 + /// \pre \ref run() must be called before using this function.
1.743 + template <typename PotentialMap>
1.744 + void potentialMap(PotentialMap &map) const {
1.745 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.746 + map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
1.747 + }
1.748 + }
1.749 +
1.750 + /// @}
1.751 +
1.752 + private:
1.753 +
1.754 + // Initialize the algorithm
1.755 + ProblemType init() {
1.756 + if (_res_node_num <= 1) return INFEASIBLE;
1.757 +
1.758 + // Check the sum of supply values
1.759 + _sum_supply = 0;
1.760 + for (int i = 0; i != _root; ++i) {
1.761 + _sum_supply += _supply[i];
1.762 + }
1.763 + if (_sum_supply > 0) return INFEASIBLE;
1.764 +
1.765 +
1.766 + // Initialize vectors
1.767 + for (int i = 0; i != _res_node_num; ++i) {
1.768 + _pi[i] = 0;
1.769 + _excess[i] = _supply[i];
1.770 + }
1.771 +
1.772 + // Remove infinite upper bounds and check negative arcs
1.773 + const Value MAX = std::numeric_limits<Value>::max();
1.774 + int last_out;
1.775 + if (_have_lower) {
1.776 + for (int i = 0; i != _root; ++i) {
1.777 + last_out = _first_out[i+1];
1.778 + for (int j = _first_out[i]; j != last_out; ++j) {
1.779 + if (_forward[j]) {
1.780 + Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
1.781 + if (c >= MAX) return UNBOUNDED;
1.782 + _excess[i] -= c;
1.783 + _excess[_target[j]] += c;
1.784 + }
1.785 + }
1.786 + }
1.787 + } else {
1.788 + for (int i = 0; i != _root; ++i) {
1.789 + last_out = _first_out[i+1];
1.790 + for (int j = _first_out[i]; j != last_out; ++j) {
1.791 + if (_forward[j] && _scost[j] < 0) {
1.792 + Value c = _upper[j];
1.793 + if (c >= MAX) return UNBOUNDED;
1.794 + _excess[i] -= c;
1.795 + _excess[_target[j]] += c;
1.796 + }
1.797 + }
1.798 + }
1.799 + }
1.800 + Value ex, max_cap = 0;
1.801 + for (int i = 0; i != _res_node_num; ++i) {
1.802 + ex = _excess[i];
1.803 + _excess[i] = 0;
1.804 + if (ex < 0) max_cap -= ex;
1.805 + }
1.806 + for (int j = 0; j != _res_arc_num; ++j) {
1.807 + if (_upper[j] >= MAX) _upper[j] = max_cap;
1.808 + }
1.809 +
1.810 + // Initialize the large cost vector and the epsilon parameter
1.811 + _epsilon = 0;
1.812 + LargeCost lc;
1.813 + for (int i = 0; i != _root; ++i) {
1.814 + last_out = _first_out[i+1];
1.815 + for (int j = _first_out[i]; j != last_out; ++j) {
1.816 + lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha;
1.817 + _cost[j] = lc;
1.818 + if (lc > _epsilon) _epsilon = lc;
1.819 + }
1.820 + }
1.821 + _epsilon /= _alpha;
1.822 +
1.823 + // Initialize maps for Circulation and remove non-zero lower bounds
1.824 + ConstMap<Arc, Value> low(0);
1.825 + typedef typename Digraph::template ArcMap<Value> ValueArcMap;
1.826 + typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
1.827 + ValueArcMap cap(_graph), flow(_graph);
1.828 + ValueNodeMap sup(_graph);
1.829 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.830 + sup[n] = _supply[_node_id[n]];
1.831 + }
1.832 + if (_have_lower) {
1.833 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.834 + int j = _arc_idf[a];
1.835 + Value c = _lower[j];
1.836 + cap[a] = _upper[j] - c;
1.837 + sup[_graph.source(a)] -= c;
1.838 + sup[_graph.target(a)] += c;
1.839 + }
1.840 + } else {
1.841 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.842 + cap[a] = _upper[_arc_idf[a]];
1.843 + }
1.844 + }
1.845 +
1.846 + _sup_node_num = 0;
1.847 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.848 + if (sup[n] > 0) ++_sup_node_num;
1.849 + }
1.850 +
1.851 + // Find a feasible flow using Circulation
1.852 + Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
1.853 + circ(_graph, low, cap, sup);
1.854 + if (!circ.flowMap(flow).run()) return INFEASIBLE;
1.855 +
1.856 + // Set residual capacities and handle GEQ supply type
1.857 + if (_sum_supply < 0) {
1.858 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.859 + Value fa = flow[a];
1.860 + _res_cap[_arc_idf[a]] = cap[a] - fa;
1.861 + _res_cap[_arc_idb[a]] = fa;
1.862 + sup[_graph.source(a)] -= fa;
1.863 + sup[_graph.target(a)] += fa;
1.864 + }
1.865 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.866 + _excess[_node_id[n]] = sup[n];
1.867 + }
1.868 + for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
1.869 + int u = _target[a];
1.870 + int ra = _reverse[a];
1.871 + _res_cap[a] = -_sum_supply + 1;
1.872 + _res_cap[ra] = -_excess[u];
1.873 + _cost[a] = 0;
1.874 + _cost[ra] = 0;
1.875 + _excess[u] = 0;
1.876 + }
1.877 + } else {
1.878 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.879 + Value fa = flow[a];
1.880 + _res_cap[_arc_idf[a]] = cap[a] - fa;
1.881 + _res_cap[_arc_idb[a]] = fa;
1.882 + }
1.883 + for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
1.884 + int ra = _reverse[a];
1.885 + _res_cap[a] = 0;
1.886 + _res_cap[ra] = 0;
1.887 + _cost[a] = 0;
1.888 + _cost[ra] = 0;
1.889 + }
1.890 + }
1.891 +
1.892 + return OPTIMAL;
1.893 + }
1.894 +
1.895 + // Execute the algorithm and transform the results
1.896 + void start(Method method) {
1.897 + // Maximum path length for partial augment
1.898 + const int MAX_PATH_LENGTH = 4;
1.899 +
1.900 + // Initialize data structures for buckets
1.901 + _max_rank = _alpha * _res_node_num;
1.902 + _buckets.resize(_max_rank);
1.903 + _bucket_next.resize(_res_node_num + 1);
1.904 + _bucket_prev.resize(_res_node_num + 1);
1.905 + _rank.resize(_res_node_num + 1);
1.906 +
1.907 + // Execute the algorithm
1.908 + switch (method) {
1.909 + case PUSH:
1.910 + startPush();
1.911 + break;
1.912 + case AUGMENT:
1.913 + startAugment(_res_node_num - 1);
1.914 + break;
1.915 + case PARTIAL_AUGMENT:
1.916 + startAugment(MAX_PATH_LENGTH);
1.917 + break;
1.918 + }
1.919 +
1.920 + // Compute node potentials for the original costs
1.921 + _arc_vec.clear();
1.922 + _cost_vec.clear();
1.923 + for (int j = 0; j != _res_arc_num; ++j) {
1.924 + if (_res_cap[j] > 0) {
1.925 + _arc_vec.push_back(IntPair(_source[j], _target[j]));
1.926 + _cost_vec.push_back(_scost[j]);
1.927 + }
1.928 + }
1.929 + _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1.930 +
1.931 + typename BellmanFord<StaticDigraph, LargeCostArcMap>
1.932 + ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);
1.933 + bf.distMap(_pi_map);
1.934 + bf.init(0);
1.935 + bf.start();
1.936 +
1.937 + // Handle non-zero lower bounds
1.938 + if (_have_lower) {
1.939 + int limit = _first_out[_root];
1.940 + for (int j = 0; j != limit; ++j) {
1.941 + if (!_forward[j]) _res_cap[j] += _lower[j];
1.942 + }
1.943 + }
1.944 + }
1.945 +
1.946 + // Initialize a cost scaling phase
1.947 + void initPhase() {
1.948 + // Saturate arcs not satisfying the optimality condition
1.949 + for (int u = 0; u != _res_node_num; ++u) {
1.950 + int last_out = _first_out[u+1];
1.951 + LargeCost pi_u = _pi[u];
1.952 + for (int a = _first_out[u]; a != last_out; ++a) {
1.953 + int v = _target[a];
1.954 + if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) {
1.955 + Value delta = _res_cap[a];
1.956 + _excess[u] -= delta;
1.957 + _excess[v] += delta;
1.958 + _res_cap[a] = 0;
1.959 + _res_cap[_reverse[a]] += delta;
1.960 + }
1.961 + }
1.962 + }
1.963 +
1.964 + // Find active nodes (i.e. nodes with positive excess)
1.965 + for (int u = 0; u != _res_node_num; ++u) {
1.966 + if (_excess[u] > 0) _active_nodes.push_back(u);
1.967 + }
1.968 +
1.969 + // Initialize the next arcs
1.970 + for (int u = 0; u != _res_node_num; ++u) {
1.971 + _next_out[u] = _first_out[u];
1.972 + }
1.973 + }
1.974 +
1.975 + // Early termination heuristic
1.976 + bool earlyTermination() {
1.977 + const double EARLY_TERM_FACTOR = 3.0;
1.978 +
1.979 + // Build a static residual graph
1.980 + _arc_vec.clear();
1.981 + _cost_vec.clear();
1.982 + for (int j = 0; j != _res_arc_num; ++j) {
1.983 + if (_res_cap[j] > 0) {
1.984 + _arc_vec.push_back(IntPair(_source[j], _target[j]));
1.985 + _cost_vec.push_back(_cost[j] + 1);
1.986 + }
1.987 + }
1.988 + _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1.989 +
1.990 + // Run Bellman-Ford algorithm to check if the current flow is optimal
1.991 + BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
1.992 + bf.init(0);
1.993 + bool done = false;
1.994 + int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num)));
1.995 + for (int i = 0; i < K && !done; ++i) {
1.996 + done = bf.processNextWeakRound();
1.997 + }
1.998 + return done;
1.999 + }
1.1000 +
1.1001 + // Global potential update heuristic
1.1002 + void globalUpdate() {
1.1003 + int bucket_end = _root + 1;
1.1004 +
1.1005 + // Initialize buckets
1.1006 + for (int r = 0; r != _max_rank; ++r) {
1.1007 + _buckets[r] = bucket_end;
1.1008 + }
1.1009 + Value total_excess = 0;
1.1010 + for (int i = 0; i != _res_node_num; ++i) {
1.1011 + if (_excess[i] < 0) {
1.1012 + _rank[i] = 0;
1.1013 + _bucket_next[i] = _buckets[0];
1.1014 + _bucket_prev[_buckets[0]] = i;
1.1015 + _buckets[0] = i;
1.1016 + } else {
1.1017 + total_excess += _excess[i];
1.1018 + _rank[i] = _max_rank;
1.1019 + }
1.1020 + }
1.1021 + if (total_excess == 0) return;
1.1022 +
1.1023 + // Search the buckets
1.1024 + int r = 0;
1.1025 + for ( ; r != _max_rank; ++r) {
1.1026 + while (_buckets[r] != bucket_end) {
1.1027 + // Remove the first node from the current bucket
1.1028 + int u = _buckets[r];
1.1029 + _buckets[r] = _bucket_next[u];
1.1030 +
1.1031 + // Search the incomming arcs of u
1.1032 + LargeCost pi_u = _pi[u];
1.1033 + int last_out = _first_out[u+1];
1.1034 + for (int a = _first_out[u]; a != last_out; ++a) {
1.1035 + int ra = _reverse[a];
1.1036 + if (_res_cap[ra] > 0) {
1.1037 + int v = _source[ra];
1.1038 + int old_rank_v = _rank[v];
1.1039 + if (r < old_rank_v) {
1.1040 + // Compute the new rank of v
1.1041 + LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
1.1042 + int new_rank_v = old_rank_v;
1.1043 + if (nrc < LargeCost(_max_rank))
1.1044 + new_rank_v = r + 1 + int(nrc);
1.1045 +
1.1046 + // Change the rank of v
1.1047 + if (new_rank_v < old_rank_v) {
1.1048 + _rank[v] = new_rank_v;
1.1049 + _next_out[v] = _first_out[v];
1.1050 +
1.1051 + // Remove v from its old bucket
1.1052 + if (old_rank_v < _max_rank) {
1.1053 + if (_buckets[old_rank_v] == v) {
1.1054 + _buckets[old_rank_v] = _bucket_next[v];
1.1055 + } else {
1.1056 + _bucket_next[_bucket_prev[v]] = _bucket_next[v];
1.1057 + _bucket_prev[_bucket_next[v]] = _bucket_prev[v];
1.1058 + }
1.1059 + }
1.1060 +
1.1061 + // Insert v to its new bucket
1.1062 + _bucket_next[v] = _buckets[new_rank_v];
1.1063 + _bucket_prev[_buckets[new_rank_v]] = v;
1.1064 + _buckets[new_rank_v] = v;
1.1065 + }
1.1066 + }
1.1067 + }
1.1068 + }
1.1069 +
1.1070 + // Finish search if there are no more active nodes
1.1071 + if (_excess[u] > 0) {
1.1072 + total_excess -= _excess[u];
1.1073 + if (total_excess <= 0) break;
1.1074 + }
1.1075 + }
1.1076 + if (total_excess <= 0) break;
1.1077 + }
1.1078 +
1.1079 + // Relabel nodes
1.1080 + for (int u = 0; u != _res_node_num; ++u) {
1.1081 + int k = std::min(_rank[u], r);
1.1082 + if (k > 0) {
1.1083 + _pi[u] -= _epsilon * k;
1.1084 + _next_out[u] = _first_out[u];
1.1085 + }
1.1086 + }
1.1087 + }
1.1088 +
1.1089 + /// Execute the algorithm performing augment and relabel operations
1.1090 + void startAugment(int max_length) {
1.1091 + // Paramters for heuristics
1.1092 + const int EARLY_TERM_EPSILON_LIMIT = 1000;
1.1093 + const double GLOBAL_UPDATE_FACTOR = 3.0;
1.1094 +
1.1095 + const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
1.1096 + (_res_node_num + _sup_node_num * _sup_node_num));
1.1097 + int next_update_limit = global_update_freq;
1.1098 +
1.1099 + int relabel_cnt = 0;
1.1100 +
1.1101 + // Perform cost scaling phases
1.1102 + std::vector<int> path;
1.1103 + for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.1104 + 1 : _epsilon / _alpha )
1.1105 + {
1.1106 + // Early termination heuristic
1.1107 + if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
1.1108 + if (earlyTermination()) break;
1.1109 + }
1.1110 +
1.1111 + // Initialize current phase
1.1112 + initPhase();
1.1113 +
1.1114 + // Perform partial augment and relabel operations
1.1115 + while (true) {
1.1116 + // Select an active node (FIFO selection)
1.1117 + while (_active_nodes.size() > 0 &&
1.1118 + _excess[_active_nodes.front()] <= 0) {
1.1119 + _active_nodes.pop_front();
1.1120 + }
1.1121 + if (_active_nodes.size() == 0) break;
1.1122 + int start = _active_nodes.front();
1.1123 +
1.1124 + // Find an augmenting path from the start node
1.1125 + path.clear();
1.1126 + int tip = start;
1.1127 + while (_excess[tip] >= 0 && int(path.size()) < max_length) {
1.1128 + int u;
1.1129 + LargeCost min_red_cost, rc, pi_tip = _pi[tip];
1.1130 + int last_out = _first_out[tip+1];
1.1131 + for (int a = _next_out[tip]; a != last_out; ++a) {
1.1132 + u = _target[a];
1.1133 + if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) {
1.1134 + path.push_back(a);
1.1135 + _next_out[tip] = a;
1.1136 + tip = u;
1.1137 + goto next_step;
1.1138 + }
1.1139 + }
1.1140 +
1.1141 + // Relabel tip node
1.1142 + min_red_cost = std::numeric_limits<LargeCost>::max();
1.1143 + if (tip != start) {
1.1144 + int ra = _reverse[path.back()];
1.1145 + min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]];
1.1146 + }
1.1147 + for (int a = _first_out[tip]; a != last_out; ++a) {
1.1148 + rc = _cost[a] + pi_tip - _pi[_target[a]];
1.1149 + if (_res_cap[a] > 0 && rc < min_red_cost) {
1.1150 + min_red_cost = rc;
1.1151 + }
1.1152 + }
1.1153 + _pi[tip] -= min_red_cost + _epsilon;
1.1154 + _next_out[tip] = _first_out[tip];
1.1155 + ++relabel_cnt;
1.1156 +
1.1157 + // Step back
1.1158 + if (tip != start) {
1.1159 + tip = _source[path.back()];
1.1160 + path.pop_back();
1.1161 + }
1.1162 +
1.1163 + next_step: ;
1.1164 + }
1.1165 +
1.1166 + // Augment along the found path (as much flow as possible)
1.1167 + Value delta;
1.1168 + int pa, u, v = start;
1.1169 + for (int i = 0; i != int(path.size()); ++i) {
1.1170 + pa = path[i];
1.1171 + u = v;
1.1172 + v = _target[pa];
1.1173 + delta = std::min(_res_cap[pa], _excess[u]);
1.1174 + _res_cap[pa] -= delta;
1.1175 + _res_cap[_reverse[pa]] += delta;
1.1176 + _excess[u] -= delta;
1.1177 + _excess[v] += delta;
1.1178 + if (_excess[v] > 0 && _excess[v] <= delta)
1.1179 + _active_nodes.push_back(v);
1.1180 + }
1.1181 +
1.1182 + // Global update heuristic
1.1183 + if (relabel_cnt >= next_update_limit) {
1.1184 + globalUpdate();
1.1185 + next_update_limit += global_update_freq;
1.1186 + }
1.1187 + }
1.1188 + }
1.1189 + }
1.1190 +
1.1191 + /// Execute the algorithm performing push and relabel operations
1.1192 + void startPush() {
1.1193 + // Paramters for heuristics
1.1194 + const int EARLY_TERM_EPSILON_LIMIT = 1000;
1.1195 + const double GLOBAL_UPDATE_FACTOR = 2.0;
1.1196 +
1.1197 + const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
1.1198 + (_res_node_num + _sup_node_num * _sup_node_num));
1.1199 + int next_update_limit = global_update_freq;
1.1200 +
1.1201 + int relabel_cnt = 0;
1.1202 +
1.1203 + // Perform cost scaling phases
1.1204 + BoolVector hyper(_res_node_num, false);
1.1205 + LargeCostVector hyper_cost(_res_node_num);
1.1206 + for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.1207 + 1 : _epsilon / _alpha )
1.1208 + {
1.1209 + // Early termination heuristic
1.1210 + if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
1.1211 + if (earlyTermination()) break;
1.1212 + }
1.1213 +
1.1214 + // Initialize current phase
1.1215 + initPhase();
1.1216 +
1.1217 + // Perform push and relabel operations
1.1218 + while (_active_nodes.size() > 0) {
1.1219 + LargeCost min_red_cost, rc, pi_n;
1.1220 + Value delta;
1.1221 + int n, t, a, last_out = _res_arc_num;
1.1222 +
1.1223 + next_node:
1.1224 + // Select an active node (FIFO selection)
1.1225 + n = _active_nodes.front();
1.1226 + last_out = _first_out[n+1];
1.1227 + pi_n = _pi[n];
1.1228 +
1.1229 + // Perform push operations if there are admissible arcs
1.1230 + if (_excess[n] > 0) {
1.1231 + for (a = _next_out[n]; a != last_out; ++a) {
1.1232 + if (_res_cap[a] > 0 &&
1.1233 + _cost[a] + pi_n - _pi[_target[a]] < 0) {
1.1234 + delta = std::min(_res_cap[a], _excess[n]);
1.1235 + t = _target[a];
1.1236 +
1.1237 + // Push-look-ahead heuristic
1.1238 + Value ahead = -_excess[t];
1.1239 + int last_out_t = _first_out[t+1];
1.1240 + LargeCost pi_t = _pi[t];
1.1241 + for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
1.1242 + if (_res_cap[ta] > 0 &&
1.1243 + _cost[ta] + pi_t - _pi[_target[ta]] < 0)
1.1244 + ahead += _res_cap[ta];
1.1245 + if (ahead >= delta) break;
1.1246 + }
1.1247 + if (ahead < 0) ahead = 0;
1.1248 +
1.1249 + // Push flow along the arc
1.1250 + if (ahead < delta && !hyper[t]) {
1.1251 + _res_cap[a] -= ahead;
1.1252 + _res_cap[_reverse[a]] += ahead;
1.1253 + _excess[n] -= ahead;
1.1254 + _excess[t] += ahead;
1.1255 + _active_nodes.push_front(t);
1.1256 + hyper[t] = true;
1.1257 + hyper_cost[t] = _cost[a] + pi_n - pi_t;
1.1258 + _next_out[n] = a;
1.1259 + goto next_node;
1.1260 + } else {
1.1261 + _res_cap[a] -= delta;
1.1262 + _res_cap[_reverse[a]] += delta;
1.1263 + _excess[n] -= delta;
1.1264 + _excess[t] += delta;
1.1265 + if (_excess[t] > 0 && _excess[t] <= delta)
1.1266 + _active_nodes.push_back(t);
1.1267 + }
1.1268 +
1.1269 + if (_excess[n] == 0) {
1.1270 + _next_out[n] = a;
1.1271 + goto remove_nodes;
1.1272 + }
1.1273 + }
1.1274 + }
1.1275 + _next_out[n] = a;
1.1276 + }
1.1277 +
1.1278 + // Relabel the node if it is still active (or hyper)
1.1279 + if (_excess[n] > 0 || hyper[n]) {
1.1280 + min_red_cost = hyper[n] ? -hyper_cost[n] :
1.1281 + std::numeric_limits<LargeCost>::max();
1.1282 + for (int a = _first_out[n]; a != last_out; ++a) {
1.1283 + rc = _cost[a] + pi_n - _pi[_target[a]];
1.1284 + if (_res_cap[a] > 0 && rc < min_red_cost) {
1.1285 + min_red_cost = rc;
1.1286 + }
1.1287 + }
1.1288 + _pi[n] -= min_red_cost + _epsilon;
1.1289 + _next_out[n] = _first_out[n];
1.1290 + hyper[n] = false;
1.1291 + ++relabel_cnt;
1.1292 + }
1.1293 +
1.1294 + // Remove nodes that are not active nor hyper
1.1295 + remove_nodes:
1.1296 + while ( _active_nodes.size() > 0 &&
1.1297 + _excess[_active_nodes.front()] <= 0 &&
1.1298 + !hyper[_active_nodes.front()] ) {
1.1299 + _active_nodes.pop_front();
1.1300 + }
1.1301 +
1.1302 + // Global update heuristic
1.1303 + if (relabel_cnt >= next_update_limit) {
1.1304 + globalUpdate();
1.1305 + for (int u = 0; u != _res_node_num; ++u)
1.1306 + hyper[u] = false;
1.1307 + next_update_limit += global_update_freq;
1.1308 + }
1.1309 + }
1.1310 + }
1.1311 + }
1.1312 +
1.1313 + }; //class CostScaling
1.1314 +
1.1315 + ///@}
1.1316 +
1.1317 +} //namespace lemon
1.1318 +
1.1319 +#endif //LEMON_COST_SCALING_H