1.1 --- a/lemon/cost_scaling.h Thu Nov 12 23:29:42 2009 +0100
1.2 +++ b/lemon/cost_scaling.h Thu Nov 12 23:30:45 2009 +0100
1.3 @@ -30,548 +30,912 @@
1.4 #include <lemon/core.h>
1.5 #include <lemon/maps.h>
1.6 #include <lemon/math.h>
1.7 -#include <lemon/adaptors.h>
1.8 +#include <lemon/static_graph.h>
1.9 #include <lemon/circulation.h>
1.10 #include <lemon/bellman_ford.h>
1.11
1.12 namespace lemon {
1.13
1.14 + /// \brief Default traits class of CostScaling algorithm.
1.15 + ///
1.16 + /// Default traits class of CostScaling algorithm.
1.17 + /// \tparam GR Digraph type.
1.18 + /// \tparam V The value type used for flow amounts, capacity bounds
1.19 + /// and supply values. By default it is \c int.
1.20 + /// \tparam C The value type used for costs and potentials.
1.21 + /// By default it is the same as \c V.
1.22 +#ifdef DOXYGEN
1.23 + template <typename GR, typename V = int, typename C = V>
1.24 +#else
1.25 + template < typename GR, typename V = int, typename C = V,
1.26 + bool integer = std::numeric_limits<C>::is_integer >
1.27 +#endif
1.28 + struct CostScalingDefaultTraits
1.29 + {
1.30 + /// The type of the digraph
1.31 + typedef GR Digraph;
1.32 + /// The type of the flow amounts, capacity bounds and supply values
1.33 + typedef V Value;
1.34 + /// The type of the arc costs
1.35 + typedef C Cost;
1.36 +
1.37 + /// \brief The large cost type used for internal computations
1.38 + ///
1.39 + /// The large cost type used for internal computations.
1.40 + /// It is \c long \c long if the \c Cost type is integer,
1.41 + /// otherwise it is \c double.
1.42 + /// \c Cost must be convertible to \c LargeCost.
1.43 + typedef double LargeCost;
1.44 + };
1.45 +
1.46 + // Default traits class for integer cost types
1.47 + template <typename GR, typename V, typename C>
1.48 + struct CostScalingDefaultTraits<GR, V, C, true>
1.49 + {
1.50 + typedef GR Digraph;
1.51 + typedef V Value;
1.52 + typedef C Cost;
1.53 +#ifdef LEMON_HAVE_LONG_LONG
1.54 + typedef long long LargeCost;
1.55 +#else
1.56 + typedef long LargeCost;
1.57 +#endif
1.58 + };
1.59 +
1.60 +
1.61 /// \addtogroup min_cost_flow_algs
1.62 /// @{
1.63
1.64 - /// \brief Implementation of the cost scaling algorithm for finding a
1.65 - /// minimum cost flow.
1.66 + /// \brief Implementation of the Cost Scaling algorithm for
1.67 + /// finding a \ref min_cost_flow "minimum cost flow".
1.68 ///
1.69 - /// \ref CostScaling implements the cost scaling algorithm performing
1.70 - /// augment/push and relabel operations for finding a minimum cost
1.71 - /// flow.
1.72 + /// \ref CostScaling implements a cost scaling algorithm that performs
1.73 + /// push/augment and relabel operations for finding a minimum cost
1.74 + /// flow. It is an efficient primal-dual solution method, which
1.75 + /// can be viewed as the generalization of the \ref Preflow
1.76 + /// "preflow push-relabel" algorithm for the maximum flow problem.
1.77 ///
1.78 - /// \tparam Digraph The digraph type the algorithm runs on.
1.79 - /// \tparam LowerMap The type of the lower bound map.
1.80 - /// \tparam CapacityMap The type of the capacity (upper bound) map.
1.81 - /// \tparam CostMap The type of the cost (length) map.
1.82 - /// \tparam SupplyMap The type of the supply map.
1.83 + /// Most of the parameters of the problem (except for the digraph)
1.84 + /// can be given using separate functions, and the algorithm can be
1.85 + /// executed using the \ref run() function. If some parameters are not
1.86 + /// specified, then default values will be used.
1.87 ///
1.88 - /// \warning
1.89 - /// - Arc capacities and costs should be \e non-negative \e integers.
1.90 - /// - Supply values should be \e signed \e integers.
1.91 - /// - The value types of the maps should be convertible to each other.
1.92 - /// - \c CostMap::Value must be signed type.
1.93 + /// \tparam GR The digraph type the algorithm runs on.
1.94 + /// \tparam V The value type used for flow amounts, capacity bounds
1.95 + /// and supply values in the algorithm. By default it is \c int.
1.96 + /// \tparam C The value type used for costs and potentials in the
1.97 + /// algorithm. By default it is the same as \c V.
1.98 ///
1.99 - /// \note Arc costs are multiplied with the number of nodes during
1.100 - /// the algorithm so overflow problems may arise more easily than with
1.101 - /// other minimum cost flow algorithms.
1.102 - /// If it is available, <tt>long long int</tt> type is used instead of
1.103 - /// <tt>long int</tt> in the inside computations.
1.104 - ///
1.105 - /// \author Peter Kovacs
1.106 - template < typename Digraph,
1.107 - typename LowerMap = typename Digraph::template ArcMap<int>,
1.108 - typename CapacityMap = typename Digraph::template ArcMap<int>,
1.109 - typename CostMap = typename Digraph::template ArcMap<int>,
1.110 - typename SupplyMap = typename Digraph::template NodeMap<int> >
1.111 + /// \warning Both value types must be signed and all input data must
1.112 + /// be integer.
1.113 + /// \warning This algorithm does not support negative costs for such
1.114 + /// arcs that have infinite upper bound.
1.115 +#ifdef DOXYGEN
1.116 + template <typename GR, typename V, typename C, typename TR>
1.117 +#else
1.118 + template < typename GR, typename V = int, typename C = V,
1.119 + typename TR = CostScalingDefaultTraits<GR, V, C> >
1.120 +#endif
1.121 class CostScaling
1.122 {
1.123 - TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
1.124 + public:
1.125
1.126 - typedef typename CapacityMap::Value Capacity;
1.127 - typedef typename CostMap::Value Cost;
1.128 - typedef typename SupplyMap::Value Supply;
1.129 - typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
1.130 - typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
1.131 + /// The type of the digraph
1.132 + typedef typename TR::Digraph Digraph;
1.133 + /// The type of the flow amounts, capacity bounds and supply values
1.134 + typedef typename TR::Value Value;
1.135 + /// The type of the arc costs
1.136 + typedef typename TR::Cost Cost;
1.137
1.138 - typedef ResidualDigraph< const Digraph,
1.139 - CapacityArcMap, CapacityArcMap > ResDigraph;
1.140 - typedef typename ResDigraph::Arc ResArc;
1.141 + /// \brief The large cost type
1.142 + ///
1.143 + /// The large cost type used for internal computations.
1.144 + /// Using the \ref CostScalingDefaultTraits "default traits class",
1.145 + /// it is \c long \c long if the \c Cost type is integer,
1.146 + /// otherwise it is \c double.
1.147 + typedef typename TR::LargeCost LargeCost;
1.148
1.149 -#if defined __GNUC__ && !defined __STRICT_ANSI__
1.150 - typedef long long int LCost;
1.151 -#else
1.152 - typedef long int LCost;
1.153 -#endif
1.154 - typedef typename Digraph::template ArcMap<LCost> LargeCostMap;
1.155 + /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
1.156 + typedef TR Traits;
1.157
1.158 public:
1.159
1.160 - /// The type of the flow map.
1.161 - typedef typename Digraph::template ArcMap<Capacity> FlowMap;
1.162 - /// The type of the potential map.
1.163 - typedef typename Digraph::template NodeMap<LCost> PotentialMap;
1.164 + /// \brief Problem type constants for the \c run() function.
1.165 + ///
1.166 + /// Enum type containing the problem type constants that can be
1.167 + /// returned by the \ref run() function of the algorithm.
1.168 + enum ProblemType {
1.169 + /// The problem has no feasible solution (flow).
1.170 + INFEASIBLE,
1.171 + /// The problem has optimal solution (i.e. it is feasible and
1.172 + /// bounded), and the algorithm has found optimal flow and node
1.173 + /// potentials (primal and dual solutions).
1.174 + OPTIMAL,
1.175 + /// The digraph contains an arc of negative cost and infinite
1.176 + /// upper bound. It means that the objective function is unbounded
1.177 + /// on that arc, however note that it could actually be bounded
1.178 + /// over the feasible flows, but this algroithm cannot handle
1.179 + /// these cases.
1.180 + UNBOUNDED
1.181 + };
1.182
1.183 private:
1.184
1.185 - /// \brief Map adaptor class for handling residual arc costs.
1.186 - ///
1.187 - /// Map adaptor class for handling residual arc costs.
1.188 - template <typename Map>
1.189 - class ResidualCostMap : public MapBase<ResArc, typename Map::Value>
1.190 - {
1.191 - private:
1.192 + TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.193
1.194 - const Map &_cost_map;
1.195 + typedef std::vector<int> IntVector;
1.196 + typedef std::vector<char> BoolVector;
1.197 + typedef std::vector<Value> ValueVector;
1.198 + typedef std::vector<Cost> CostVector;
1.199 + typedef std::vector<LargeCost> LargeCostVector;
1.200
1.201 + private:
1.202 +
1.203 + template <typename KT, typename VT>
1.204 + class VectorMap {
1.205 public:
1.206 -
1.207 - ///\e
1.208 - ResidualCostMap(const Map &cost_map) :
1.209 - _cost_map(cost_map) {}
1.210 -
1.211 - ///\e
1.212 - inline typename Map::Value operator[](const ResArc &e) const {
1.213 - return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
1.214 + typedef KT Key;
1.215 + typedef VT Value;
1.216 +
1.217 + VectorMap(std::vector<Value>& v) : _v(v) {}
1.218 +
1.219 + const Value& operator[](const Key& key) const {
1.220 + return _v[StaticDigraph::id(key)];
1.221 }
1.222
1.223 - }; //class ResidualCostMap
1.224 -
1.225 - /// \brief Map adaptor class for handling reduced arc costs.
1.226 - ///
1.227 - /// Map adaptor class for handling reduced arc costs.
1.228 - class ReducedCostMap : public MapBase<Arc, LCost>
1.229 - {
1.230 - private:
1.231 -
1.232 - const Digraph &_gr;
1.233 - const LargeCostMap &_cost_map;
1.234 - const PotentialMap &_pot_map;
1.235 -
1.236 - public:
1.237 -
1.238 - ///\e
1.239 - ReducedCostMap( const Digraph &gr,
1.240 - const LargeCostMap &cost_map,
1.241 - const PotentialMap &pot_map ) :
1.242 - _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
1.243 -
1.244 - ///\e
1.245 - inline LCost operator[](const Arc &e) const {
1.246 - return _cost_map[e] + _pot_map[_gr.source(e)]
1.247 - - _pot_map[_gr.target(e)];
1.248 + Value& operator[](const Key& key) {
1.249 + return _v[StaticDigraph::id(key)];
1.250 + }
1.251 +
1.252 + void set(const Key& key, const Value& val) {
1.253 + _v[StaticDigraph::id(key)] = val;
1.254 }
1.255
1.256 - }; //class ReducedCostMap
1.257 + private:
1.258 + std::vector<Value>& _v;
1.259 + };
1.260 +
1.261 + typedef VectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
1.262 + typedef VectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
1.263
1.264 private:
1.265
1.266 - // The digraph the algorithm runs on
1.267 - const Digraph &_graph;
1.268 - // The original lower bound map
1.269 - const LowerMap *_lower;
1.270 - // The modified capacity map
1.271 - CapacityArcMap _capacity;
1.272 - // The original cost map
1.273 - const CostMap &_orig_cost;
1.274 - // The scaled cost map
1.275 - LargeCostMap _cost;
1.276 - // The modified supply map
1.277 - SupplyNodeMap _supply;
1.278 - bool _valid_supply;
1.279 + // Data related to the underlying digraph
1.280 + const GR &_graph;
1.281 + int _node_num;
1.282 + int _arc_num;
1.283 + int _res_node_num;
1.284 + int _res_arc_num;
1.285 + int _root;
1.286
1.287 - // Arc map of the current flow
1.288 - FlowMap *_flow;
1.289 - bool _local_flow;
1.290 - // Node map of the current potentials
1.291 - PotentialMap *_potential;
1.292 - bool _local_potential;
1.293 + // Parameters of the problem
1.294 + bool _have_lower;
1.295 + Value _sum_supply;
1.296
1.297 - // The residual cost map
1.298 - ResidualCostMap<LargeCostMap> _res_cost;
1.299 - // The residual digraph
1.300 - ResDigraph *_res_graph;
1.301 - // The reduced cost map
1.302 - ReducedCostMap *_red_cost;
1.303 - // The excess map
1.304 - SupplyNodeMap _excess;
1.305 - // The epsilon parameter used for cost scaling
1.306 - LCost _epsilon;
1.307 - // The scaling factor
1.308 + // Data structures for storing the digraph
1.309 + IntNodeMap _node_id;
1.310 + IntArcMap _arc_idf;
1.311 + IntArcMap _arc_idb;
1.312 + IntVector _first_out;
1.313 + BoolVector _forward;
1.314 + IntVector _source;
1.315 + IntVector _target;
1.316 + IntVector _reverse;
1.317 +
1.318 + // Node and arc data
1.319 + ValueVector _lower;
1.320 + ValueVector _upper;
1.321 + CostVector _scost;
1.322 + ValueVector _supply;
1.323 +
1.324 + ValueVector _res_cap;
1.325 + LargeCostVector _cost;
1.326 + LargeCostVector _pi;
1.327 + ValueVector _excess;
1.328 + IntVector _next_out;
1.329 + std::deque<int> _active_nodes;
1.330 +
1.331 + // Data for scaling
1.332 + LargeCost _epsilon;
1.333 int _alpha;
1.334
1.335 + // Data for a StaticDigraph structure
1.336 + typedef std::pair<int, int> IntPair;
1.337 + StaticDigraph _sgr;
1.338 + std::vector<IntPair> _arc_vec;
1.339 + std::vector<LargeCost> _cost_vec;
1.340 + LargeCostArcMap _cost_map;
1.341 + LargeCostNodeMap _pi_map;
1.342 +
1.343 + public:
1.344 +
1.345 + /// \brief Constant for infinite upper bounds (capacities).
1.346 + ///
1.347 + /// Constant for infinite upper bounds (capacities).
1.348 + /// It is \c std::numeric_limits<Value>::infinity() if available,
1.349 + /// \c std::numeric_limits<Value>::max() otherwise.
1.350 + const Value INF;
1.351 +
1.352 public:
1.353
1.354 - /// \brief General constructor (with lower bounds).
1.355 + /// \name Named Template Parameters
1.356 + /// @{
1.357 +
1.358 + template <typename T>
1.359 + struct SetLargeCostTraits : public Traits {
1.360 + typedef T LargeCost;
1.361 + };
1.362 +
1.363 + /// \brief \ref named-templ-param "Named parameter" for setting
1.364 + /// \c LargeCost type.
1.365 ///
1.366 - /// General constructor (with lower bounds).
1.367 + /// \ref named-templ-param "Named parameter" for setting \c LargeCost
1.368 + /// type, which is used for internal computations in the algorithm.
1.369 + /// \c Cost must be convertible to \c LargeCost.
1.370 + template <typename T>
1.371 + struct SetLargeCost
1.372 + : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {
1.373 + typedef CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;
1.374 + };
1.375 +
1.376 + /// @}
1.377 +
1.378 + public:
1.379 +
1.380 + /// \brief Constructor.
1.381 ///
1.382 - /// \param digraph The digraph the algorithm runs on.
1.383 - /// \param lower The lower bounds of the arcs.
1.384 - /// \param capacity The capacities (upper bounds) of the arcs.
1.385 - /// \param cost The cost (length) values of the arcs.
1.386 - /// \param supply The supply values of the nodes (signed).
1.387 - CostScaling( const Digraph &digraph,
1.388 - const LowerMap &lower,
1.389 - const CapacityMap &capacity,
1.390 - const CostMap &cost,
1.391 - const SupplyMap &supply ) :
1.392 - _graph(digraph), _lower(&lower), _capacity(digraph), _orig_cost(cost),
1.393 - _cost(digraph), _supply(digraph), _flow(NULL), _local_flow(false),
1.394 - _potential(NULL), _local_potential(false), _res_cost(_cost),
1.395 - _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
1.396 + /// The constructor of the class.
1.397 + ///
1.398 + /// \param graph The digraph the algorithm runs on.
1.399 + CostScaling(const GR& graph) :
1.400 + _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
1.401 + _cost_map(_cost_vec), _pi_map(_pi),
1.402 + INF(std::numeric_limits<Value>::has_infinity ?
1.403 + std::numeric_limits<Value>::infinity() :
1.404 + std::numeric_limits<Value>::max())
1.405 {
1.406 - // Check the sum of supply values
1.407 - Supply sum = 0;
1.408 - for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
1.409 - _valid_supply = sum == 0;
1.410 + // Check the value types
1.411 + LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
1.412 + "The flow type of CostScaling must be signed");
1.413 + LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
1.414 + "The cost type of CostScaling must be signed");
1.415 +
1.416 + // Resize vectors
1.417 + _node_num = countNodes(_graph);
1.418 + _arc_num = countArcs(_graph);
1.419 + _res_node_num = _node_num + 1;
1.420 + _res_arc_num = 2 * (_arc_num + _node_num);
1.421 + _root = _node_num;
1.422 +
1.423 + _first_out.resize(_res_node_num + 1);
1.424 + _forward.resize(_res_arc_num);
1.425 + _source.resize(_res_arc_num);
1.426 + _target.resize(_res_arc_num);
1.427 + _reverse.resize(_res_arc_num);
1.428 +
1.429 + _lower.resize(_res_arc_num);
1.430 + _upper.resize(_res_arc_num);
1.431 + _scost.resize(_res_arc_num);
1.432 + _supply.resize(_res_node_num);
1.433
1.434 - for (ArcIt e(_graph); e != INVALID; ++e) _capacity[e] = capacity[e];
1.435 - for (NodeIt n(_graph); n != INVALID; ++n) _supply[n] = supply[n];
1.436 + _res_cap.resize(_res_arc_num);
1.437 + _cost.resize(_res_arc_num);
1.438 + _pi.resize(_res_node_num);
1.439 + _excess.resize(_res_node_num);
1.440 + _next_out.resize(_res_node_num);
1.441
1.442 - // Remove non-zero lower bounds
1.443 - for (ArcIt e(_graph); e != INVALID; ++e) {
1.444 - if (lower[e] != 0) {
1.445 - _capacity[e] -= lower[e];
1.446 - _supply[_graph.source(e)] -= lower[e];
1.447 - _supply[_graph.target(e)] += lower[e];
1.448 + _arc_vec.reserve(_res_arc_num);
1.449 + _cost_vec.reserve(_res_arc_num);
1.450 +
1.451 + // Copy the graph
1.452 + int i = 0, j = 0, k = 2 * _arc_num + _node_num;
1.453 + for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.454 + _node_id[n] = i;
1.455 + }
1.456 + i = 0;
1.457 + for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.458 + _first_out[i] = j;
1.459 + for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
1.460 + _arc_idf[a] = j;
1.461 + _forward[j] = true;
1.462 + _source[j] = i;
1.463 + _target[j] = _node_id[_graph.runningNode(a)];
1.464 }
1.465 + for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
1.466 + _arc_idb[a] = j;
1.467 + _forward[j] = false;
1.468 + _source[j] = i;
1.469 + _target[j] = _node_id[_graph.runningNode(a)];
1.470 + }
1.471 + _forward[j] = false;
1.472 + _source[j] = i;
1.473 + _target[j] = _root;
1.474 + _reverse[j] = k;
1.475 + _forward[k] = true;
1.476 + _source[k] = _root;
1.477 + _target[k] = i;
1.478 + _reverse[k] = j;
1.479 + ++j; ++k;
1.480 }
1.481 - }
1.482 -/*
1.483 - /// \brief General constructor (without lower bounds).
1.484 - ///
1.485 - /// General constructor (without lower bounds).
1.486 - ///
1.487 - /// \param digraph The digraph the algorithm runs on.
1.488 - /// \param capacity The capacities (upper bounds) of the arcs.
1.489 - /// \param cost The cost (length) values of the arcs.
1.490 - /// \param supply The supply values of the nodes (signed).
1.491 - CostScaling( const Digraph &digraph,
1.492 - const CapacityMap &capacity,
1.493 - const CostMap &cost,
1.494 - const SupplyMap &supply ) :
1.495 - _graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
1.496 - _cost(digraph), _supply(supply), _flow(NULL), _local_flow(false),
1.497 - _potential(NULL), _local_potential(false), _res_cost(_cost),
1.498 - _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
1.499 - {
1.500 - // Check the sum of supply values
1.501 - Supply sum = 0;
1.502 - for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
1.503 - _valid_supply = sum == 0;
1.504 + _first_out[i] = j;
1.505 + _first_out[_res_node_num] = k;
1.506 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.507 + int fi = _arc_idf[a];
1.508 + int bi = _arc_idb[a];
1.509 + _reverse[fi] = bi;
1.510 + _reverse[bi] = fi;
1.511 + }
1.512 +
1.513 + // Reset parameters
1.514 + reset();
1.515 }
1.516
1.517 - /// \brief Simple constructor (with lower bounds).
1.518 + /// \name Parameters
1.519 + /// The parameters of the algorithm can be specified using these
1.520 + /// functions.
1.521 +
1.522 + /// @{
1.523 +
1.524 + /// \brief Set the lower bounds on the arcs.
1.525 ///
1.526 - /// Simple constructor (with lower bounds).
1.527 + /// This function sets the lower bounds on the arcs.
1.528 + /// If it is not used before calling \ref run(), the lower bounds
1.529 + /// will be set to zero on all arcs.
1.530 ///
1.531 - /// \param digraph The digraph the algorithm runs on.
1.532 - /// \param lower The lower bounds of the arcs.
1.533 - /// \param capacity The capacities (upper bounds) of the arcs.
1.534 - /// \param cost The cost (length) values of the arcs.
1.535 - /// \param s The source node.
1.536 - /// \param t The target node.
1.537 - /// \param flow_value The required amount of flow from node \c s
1.538 - /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.539 - CostScaling( const Digraph &digraph,
1.540 - const LowerMap &lower,
1.541 - const CapacityMap &capacity,
1.542 - const CostMap &cost,
1.543 - Node s, Node t,
1.544 - Supply flow_value ) :
1.545 - _graph(digraph), _lower(&lower), _capacity(capacity), _orig_cost(cost),
1.546 - _cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false),
1.547 - _potential(NULL), _local_potential(false), _res_cost(_cost),
1.548 - _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
1.549 - {
1.550 - // Remove non-zero lower bounds
1.551 - _supply[s] = flow_value;
1.552 - _supply[t] = -flow_value;
1.553 - for (ArcIt e(_graph); e != INVALID; ++e) {
1.554 - if (lower[e] != 0) {
1.555 - _capacity[e] -= lower[e];
1.556 - _supply[_graph.source(e)] -= lower[e];
1.557 - _supply[_graph.target(e)] += lower[e];
1.558 - }
1.559 + /// \param map An arc map storing the lower bounds.
1.560 + /// Its \c Value type must be convertible to the \c Value type
1.561 + /// of the algorithm.
1.562 + ///
1.563 + /// \return <tt>(*this)</tt>
1.564 + template <typename LowerMap>
1.565 + CostScaling& lowerMap(const LowerMap& map) {
1.566 + _have_lower = true;
1.567 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.568 + _lower[_arc_idf[a]] = map[a];
1.569 + _lower[_arc_idb[a]] = map[a];
1.570 }
1.571 - _valid_supply = true;
1.572 - }
1.573 -
1.574 - /// \brief Simple constructor (without lower bounds).
1.575 - ///
1.576 - /// Simple constructor (without lower bounds).
1.577 - ///
1.578 - /// \param digraph The digraph the algorithm runs on.
1.579 - /// \param capacity The capacities (upper bounds) of the arcs.
1.580 - /// \param cost The cost (length) values of the arcs.
1.581 - /// \param s The source node.
1.582 - /// \param t The target node.
1.583 - /// \param flow_value The required amount of flow from node \c s
1.584 - /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.585 - CostScaling( const Digraph &digraph,
1.586 - const CapacityMap &capacity,
1.587 - const CostMap &cost,
1.588 - Node s, Node t,
1.589 - Supply flow_value ) :
1.590 - _graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
1.591 - _cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false),
1.592 - _potential(NULL), _local_potential(false), _res_cost(_cost),
1.593 - _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
1.594 - {
1.595 - _supply[s] = flow_value;
1.596 - _supply[t] = -flow_value;
1.597 - _valid_supply = true;
1.598 - }
1.599 -*/
1.600 - /// Destructor.
1.601 - ~CostScaling() {
1.602 - if (_local_flow) delete _flow;
1.603 - if (_local_potential) delete _potential;
1.604 - delete _res_graph;
1.605 - delete _red_cost;
1.606 - }
1.607 -
1.608 - /// \brief Set the flow map.
1.609 - ///
1.610 - /// Set the flow map.
1.611 - ///
1.612 - /// \return \c (*this)
1.613 - CostScaling& flowMap(FlowMap &map) {
1.614 - if (_local_flow) {
1.615 - delete _flow;
1.616 - _local_flow = false;
1.617 - }
1.618 - _flow = ↦
1.619 return *this;
1.620 }
1.621
1.622 - /// \brief Set the potential map.
1.623 + /// \brief Set the upper bounds (capacities) on the arcs.
1.624 ///
1.625 - /// Set the potential map.
1.626 + /// This function sets the upper bounds (capacities) on the arcs.
1.627 + /// If it is not used before calling \ref run(), the upper bounds
1.628 + /// will be set to \ref INF on all arcs (i.e. the flow value will be
1.629 + /// unbounded from above on each arc).
1.630 ///
1.631 - /// \return \c (*this)
1.632 - CostScaling& potentialMap(PotentialMap &map) {
1.633 - if (_local_potential) {
1.634 - delete _potential;
1.635 - _local_potential = false;
1.636 + /// \param map An arc map storing the upper bounds.
1.637 + /// Its \c Value type must be convertible to the \c Value type
1.638 + /// of the algorithm.
1.639 + ///
1.640 + /// \return <tt>(*this)</tt>
1.641 + template<typename UpperMap>
1.642 + CostScaling& upperMap(const UpperMap& map) {
1.643 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.644 + _upper[_arc_idf[a]] = map[a];
1.645 }
1.646 - _potential = ↦
1.647 return *this;
1.648 }
1.649
1.650 + /// \brief Set the costs of the arcs.
1.651 + ///
1.652 + /// This function sets the costs of the arcs.
1.653 + /// If it is not used before calling \ref run(), the costs
1.654 + /// will be set to \c 1 on all arcs.
1.655 + ///
1.656 + /// \param map An arc map storing the costs.
1.657 + /// Its \c Value type must be convertible to the \c Cost type
1.658 + /// of the algorithm.
1.659 + ///
1.660 + /// \return <tt>(*this)</tt>
1.661 + template<typename CostMap>
1.662 + CostScaling& costMap(const CostMap& map) {
1.663 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.664 + _scost[_arc_idf[a]] = map[a];
1.665 + _scost[_arc_idb[a]] = -map[a];
1.666 + }
1.667 + return *this;
1.668 + }
1.669 +
1.670 + /// \brief Set the supply values of the nodes.
1.671 + ///
1.672 + /// This function sets the supply values of the nodes.
1.673 + /// If neither this function nor \ref stSupply() is used before
1.674 + /// calling \ref run(), the supply of each node will be set to zero.
1.675 + ///
1.676 + /// \param map A node map storing the supply values.
1.677 + /// Its \c Value type must be convertible to the \c Value type
1.678 + /// of the algorithm.
1.679 + ///
1.680 + /// \return <tt>(*this)</tt>
1.681 + template<typename SupplyMap>
1.682 + CostScaling& supplyMap(const SupplyMap& map) {
1.683 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.684 + _supply[_node_id[n]] = map[n];
1.685 + }
1.686 + return *this;
1.687 + }
1.688 +
1.689 + /// \brief Set single source and target nodes and a supply value.
1.690 + ///
1.691 + /// This function sets a single source node and a single target node
1.692 + /// and the required flow value.
1.693 + /// If neither this function nor \ref supplyMap() is used before
1.694 + /// calling \ref run(), the supply of each node will be set to zero.
1.695 + ///
1.696 + /// Using this function has the same effect as using \ref supplyMap()
1.697 + /// with such a map in which \c k is assigned to \c s, \c -k is
1.698 + /// assigned to \c t and all other nodes have zero supply value.
1.699 + ///
1.700 + /// \param s The source node.
1.701 + /// \param t The target node.
1.702 + /// \param k The required amount of flow from node \c s to node \c t
1.703 + /// (i.e. the supply of \c s and the demand of \c t).
1.704 + ///
1.705 + /// \return <tt>(*this)</tt>
1.706 + CostScaling& stSupply(const Node& s, const Node& t, Value k) {
1.707 + for (int i = 0; i != _res_node_num; ++i) {
1.708 + _supply[i] = 0;
1.709 + }
1.710 + _supply[_node_id[s]] = k;
1.711 + _supply[_node_id[t]] = -k;
1.712 + return *this;
1.713 + }
1.714 +
1.715 + /// @}
1.716 +
1.717 /// \name Execution control
1.718 + /// The algorithm can be executed using \ref run().
1.719
1.720 /// @{
1.721
1.722 /// \brief Run the algorithm.
1.723 ///
1.724 - /// Run the algorithm.
1.725 + /// This function runs the algorithm.
1.726 + /// The paramters can be specified using functions \ref lowerMap(),
1.727 + /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
1.728 + /// For example,
1.729 + /// \code
1.730 + /// CostScaling<ListDigraph> cs(graph);
1.731 + /// cs.lowerMap(lower).upperMap(upper).costMap(cost)
1.732 + /// .supplyMap(sup).run();
1.733 + /// \endcode
1.734 + ///
1.735 + /// This function can be called more than once. All the parameters
1.736 + /// that have been given are kept for the next call, unless
1.737 + /// \ref reset() is called, thus only the modified parameters
1.738 + /// have to be set again. See \ref reset() for examples.
1.739 + /// However the underlying digraph must not be modified after this
1.740 + /// class have been constructed, since it copies the digraph.
1.741 ///
1.742 /// \param partial_augment By default the algorithm performs
1.743 /// partial augment and relabel operations in the cost scaling
1.744 /// phases. Set this parameter to \c false for using local push and
1.745 /// relabel operations instead.
1.746 ///
1.747 - /// \return \c true if a feasible flow can be found.
1.748 - bool run(bool partial_augment = true) {
1.749 - if (partial_augment) {
1.750 - return init() && startPartialAugment();
1.751 - } else {
1.752 - return init() && startPushRelabel();
1.753 + /// \return \c INFEASIBLE if no feasible flow exists,
1.754 + /// \n \c OPTIMAL if the problem has optimal solution
1.755 + /// (i.e. it is feasible and bounded), and the algorithm has found
1.756 + /// optimal flow and node potentials (primal and dual solutions),
1.757 + /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
1.758 + /// and infinite upper bound. It means that the objective function
1.759 + /// is unbounded on that arc, however note that it could actually be
1.760 + /// bounded over the feasible flows, but this algroithm cannot handle
1.761 + /// these cases.
1.762 + ///
1.763 + /// \see ProblemType
1.764 + ProblemType run(bool partial_augment = true) {
1.765 + ProblemType pt = init();
1.766 + if (pt != OPTIMAL) return pt;
1.767 + start(partial_augment);
1.768 + return OPTIMAL;
1.769 + }
1.770 +
1.771 + /// \brief Reset all the parameters that have been given before.
1.772 + ///
1.773 + /// This function resets all the paramaters that have been given
1.774 + /// before using functions \ref lowerMap(), \ref upperMap(),
1.775 + /// \ref costMap(), \ref supplyMap(), \ref stSupply().
1.776 + ///
1.777 + /// It is useful for multiple run() calls. If this function is not
1.778 + /// used, all the parameters given before are kept for the next
1.779 + /// \ref run() call.
1.780 + /// However the underlying digraph must not be modified after this
1.781 + /// class have been constructed, since it copies and extends the graph.
1.782 + ///
1.783 + /// For example,
1.784 + /// \code
1.785 + /// CostScaling<ListDigraph> cs(graph);
1.786 + ///
1.787 + /// // First run
1.788 + /// cs.lowerMap(lower).upperMap(upper).costMap(cost)
1.789 + /// .supplyMap(sup).run();
1.790 + ///
1.791 + /// // Run again with modified cost map (reset() is not called,
1.792 + /// // so only the cost map have to be set again)
1.793 + /// cost[e] += 100;
1.794 + /// cs.costMap(cost).run();
1.795 + ///
1.796 + /// // Run again from scratch using reset()
1.797 + /// // (the lower bounds will be set to zero on all arcs)
1.798 + /// cs.reset();
1.799 + /// cs.upperMap(capacity).costMap(cost)
1.800 + /// .supplyMap(sup).run();
1.801 + /// \endcode
1.802 + ///
1.803 + /// \return <tt>(*this)</tt>
1.804 + CostScaling& reset() {
1.805 + for (int i = 0; i != _res_node_num; ++i) {
1.806 + _supply[i] = 0;
1.807 }
1.808 + int limit = _first_out[_root];
1.809 + for (int j = 0; j != limit; ++j) {
1.810 + _lower[j] = 0;
1.811 + _upper[j] = INF;
1.812 + _scost[j] = _forward[j] ? 1 : -1;
1.813 + }
1.814 + for (int j = limit; j != _res_arc_num; ++j) {
1.815 + _lower[j] = 0;
1.816 + _upper[j] = INF;
1.817 + _scost[j] = 0;
1.818 + _scost[_reverse[j]] = 0;
1.819 + }
1.820 + _have_lower = false;
1.821 + return *this;
1.822 }
1.823
1.824 /// @}
1.825
1.826 /// \name Query Functions
1.827 - /// The result of the algorithm can be obtained using these
1.828 + /// The results of the algorithm can be obtained using these
1.829 /// functions.\n
1.830 - /// \ref lemon::CostScaling::run() "run()" must be called before
1.831 - /// using them.
1.832 + /// The \ref run() function must be called before using them.
1.833
1.834 /// @{
1.835
1.836 - /// \brief Return a const reference to the arc map storing the
1.837 - /// found flow.
1.838 + /// \brief Return the total cost of the found flow.
1.839 ///
1.840 - /// Return a const reference to the arc map storing the found flow.
1.841 + /// This function returns the total cost of the found flow.
1.842 + /// Its complexity is O(e).
1.843 + ///
1.844 + /// \note The return type of the function can be specified as a
1.845 + /// template parameter. For example,
1.846 + /// \code
1.847 + /// cs.totalCost<double>();
1.848 + /// \endcode
1.849 + /// It is useful if the total cost cannot be stored in the \c Cost
1.850 + /// type of the algorithm, which is the default return type of the
1.851 + /// function.
1.852 ///
1.853 /// \pre \ref run() must be called before using this function.
1.854 - const FlowMap& flowMap() const {
1.855 - return *_flow;
1.856 + template <typename Number>
1.857 + Number totalCost() const {
1.858 + Number c = 0;
1.859 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.860 + int i = _arc_idb[a];
1.861 + c += static_cast<Number>(_res_cap[i]) *
1.862 + (-static_cast<Number>(_scost[i]));
1.863 + }
1.864 + return c;
1.865 }
1.866
1.867 - /// \brief Return a const reference to the node map storing the
1.868 - /// found potentials (the dual solution).
1.869 - ///
1.870 - /// Return a const reference to the node map storing the found
1.871 - /// potentials (the dual solution).
1.872 - ///
1.873 - /// \pre \ref run() must be called before using this function.
1.874 - const PotentialMap& potentialMap() const {
1.875 - return *_potential;
1.876 +#ifndef DOXYGEN
1.877 + Cost totalCost() const {
1.878 + return totalCost<Cost>();
1.879 }
1.880 +#endif
1.881
1.882 /// \brief Return the flow on the given arc.
1.883 ///
1.884 - /// Return the flow on the given arc.
1.885 + /// This function returns the flow on the given arc.
1.886 ///
1.887 /// \pre \ref run() must be called before using this function.
1.888 - Capacity flow(const Arc& arc) const {
1.889 - return (*_flow)[arc];
1.890 + Value flow(const Arc& a) const {
1.891 + return _res_cap[_arc_idb[a]];
1.892 }
1.893
1.894 - /// \brief Return the potential of the given node.
1.895 + /// \brief Return the flow map (the primal solution).
1.896 ///
1.897 - /// Return the potential of the given node.
1.898 + /// This function copies the flow value on each arc into the given
1.899 + /// map. The \c Value type of the algorithm must be convertible to
1.900 + /// the \c Value type of the map.
1.901 ///
1.902 /// \pre \ref run() must be called before using this function.
1.903 - Cost potential(const Node& node) const {
1.904 - return (*_potential)[node];
1.905 + template <typename FlowMap>
1.906 + void flowMap(FlowMap &map) const {
1.907 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.908 + map.set(a, _res_cap[_arc_idb[a]]);
1.909 + }
1.910 }
1.911
1.912 - /// \brief Return the total cost of the found flow.
1.913 + /// \brief Return the potential (dual value) of the given node.
1.914 ///
1.915 - /// Return the total cost of the found flow. The complexity of the
1.916 - /// function is \f$ O(e) \f$.
1.917 + /// This function returns the potential (dual value) of the
1.918 + /// given node.
1.919 ///
1.920 /// \pre \ref run() must be called before using this function.
1.921 - Cost totalCost() const {
1.922 - Cost c = 0;
1.923 - for (ArcIt e(_graph); e != INVALID; ++e)
1.924 - c += (*_flow)[e] * _orig_cost[e];
1.925 - return c;
1.926 + Cost potential(const Node& n) const {
1.927 + return static_cast<Cost>(_pi[_node_id[n]]);
1.928 + }
1.929 +
1.930 + /// \brief Return the potential map (the dual solution).
1.931 + ///
1.932 + /// This function copies the potential (dual value) of each node
1.933 + /// into the given map.
1.934 + /// The \c Cost type of the algorithm must be convertible to the
1.935 + /// \c Value type of the map.
1.936 + ///
1.937 + /// \pre \ref run() must be called before using this function.
1.938 + template <typename PotentialMap>
1.939 + void potentialMap(PotentialMap &map) const {
1.940 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.941 + map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
1.942 + }
1.943 }
1.944
1.945 /// @}
1.946
1.947 private:
1.948
1.949 - /// Initialize the algorithm.
1.950 - bool init() {
1.951 - if (!_valid_supply) return false;
1.952 - // The scaling factor
1.953 + // Initialize the algorithm
1.954 + ProblemType init() {
1.955 + if (_res_node_num == 0) return INFEASIBLE;
1.956 +
1.957 + // Scaling factor
1.958 _alpha = 8;
1.959
1.960 - // Initialize flow and potential maps
1.961 - if (!_flow) {
1.962 - _flow = new FlowMap(_graph);
1.963 - _local_flow = true;
1.964 + // Check the sum of supply values
1.965 + _sum_supply = 0;
1.966 + for (int i = 0; i != _root; ++i) {
1.967 + _sum_supply += _supply[i];
1.968 }
1.969 - if (!_potential) {
1.970 - _potential = new PotentialMap(_graph);
1.971 - _local_potential = true;
1.972 + if (_sum_supply > 0) return INFEASIBLE;
1.973 +
1.974 +
1.975 + // Initialize vectors
1.976 + for (int i = 0; i != _res_node_num; ++i) {
1.977 + _pi[i] = 0;
1.978 + _excess[i] = _supply[i];
1.979 + }
1.980 +
1.981 + // Remove infinite upper bounds and check negative arcs
1.982 + const Value MAX = std::numeric_limits<Value>::max();
1.983 + int last_out;
1.984 + if (_have_lower) {
1.985 + for (int i = 0; i != _root; ++i) {
1.986 + last_out = _first_out[i+1];
1.987 + for (int j = _first_out[i]; j != last_out; ++j) {
1.988 + if (_forward[j]) {
1.989 + Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
1.990 + if (c >= MAX) return UNBOUNDED;
1.991 + _excess[i] -= c;
1.992 + _excess[_target[j]] += c;
1.993 + }
1.994 + }
1.995 + }
1.996 + } else {
1.997 + for (int i = 0; i != _root; ++i) {
1.998 + last_out = _first_out[i+1];
1.999 + for (int j = _first_out[i]; j != last_out; ++j) {
1.1000 + if (_forward[j] && _scost[j] < 0) {
1.1001 + Value c = _upper[j];
1.1002 + if (c >= MAX) return UNBOUNDED;
1.1003 + _excess[i] -= c;
1.1004 + _excess[_target[j]] += c;
1.1005 + }
1.1006 + }
1.1007 + }
1.1008 + }
1.1009 + Value ex, max_cap = 0;
1.1010 + for (int i = 0; i != _res_node_num; ++i) {
1.1011 + ex = _excess[i];
1.1012 + _excess[i] = 0;
1.1013 + if (ex < 0) max_cap -= ex;
1.1014 + }
1.1015 + for (int j = 0; j != _res_arc_num; ++j) {
1.1016 + if (_upper[j] >= MAX) _upper[j] = max_cap;
1.1017 }
1.1018
1.1019 - _red_cost = new ReducedCostMap(_graph, _cost, *_potential);
1.1020 - _res_graph = new ResDigraph(_graph, _capacity, *_flow);
1.1021 + // Initialize the large cost vector and the epsilon parameter
1.1022 + _epsilon = 0;
1.1023 + LargeCost lc;
1.1024 + for (int i = 0; i != _root; ++i) {
1.1025 + last_out = _first_out[i+1];
1.1026 + for (int j = _first_out[i]; j != last_out; ++j) {
1.1027 + lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha;
1.1028 + _cost[j] = lc;
1.1029 + if (lc > _epsilon) _epsilon = lc;
1.1030 + }
1.1031 + }
1.1032 + _epsilon /= _alpha;
1.1033
1.1034 - // Initialize the scaled cost map and the epsilon parameter
1.1035 - Cost max_cost = 0;
1.1036 - int node_num = countNodes(_graph);
1.1037 - for (ArcIt e(_graph); e != INVALID; ++e) {
1.1038 - _cost[e] = LCost(_orig_cost[e]) * node_num * _alpha;
1.1039 - if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e];
1.1040 + // Initialize maps for Circulation and remove non-zero lower bounds
1.1041 + ConstMap<Arc, Value> low(0);
1.1042 + typedef typename Digraph::template ArcMap<Value> ValueArcMap;
1.1043 + typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
1.1044 + ValueArcMap cap(_graph), flow(_graph);
1.1045 + ValueNodeMap sup(_graph);
1.1046 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1047 + sup[n] = _supply[_node_id[n]];
1.1048 }
1.1049 - _epsilon = max_cost * node_num;
1.1050 + if (_have_lower) {
1.1051 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.1052 + int j = _arc_idf[a];
1.1053 + Value c = _lower[j];
1.1054 + cap[a] = _upper[j] - c;
1.1055 + sup[_graph.source(a)] -= c;
1.1056 + sup[_graph.target(a)] += c;
1.1057 + }
1.1058 + } else {
1.1059 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.1060 + cap[a] = _upper[_arc_idf[a]];
1.1061 + }
1.1062 + }
1.1063
1.1064 // Find a feasible flow using Circulation
1.1065 - Circulation< Digraph, ConstMap<Arc, Capacity>, CapacityArcMap,
1.1066 - SupplyMap >
1.1067 - circulation( _graph, constMap<Arc>(Capacity(0)), _capacity,
1.1068 - _supply );
1.1069 - return circulation.flowMap(*_flow).run();
1.1070 + Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
1.1071 + circ(_graph, low, cap, sup);
1.1072 + if (!circ.flowMap(flow).run()) return INFEASIBLE;
1.1073 +
1.1074 + // Set residual capacities and handle GEQ supply type
1.1075 + if (_sum_supply < 0) {
1.1076 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.1077 + Value fa = flow[a];
1.1078 + _res_cap[_arc_idf[a]] = cap[a] - fa;
1.1079 + _res_cap[_arc_idb[a]] = fa;
1.1080 + sup[_graph.source(a)] -= fa;
1.1081 + sup[_graph.target(a)] += fa;
1.1082 + }
1.1083 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1084 + _excess[_node_id[n]] = sup[n];
1.1085 + }
1.1086 + for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
1.1087 + int u = _target[a];
1.1088 + int ra = _reverse[a];
1.1089 + _res_cap[a] = -_sum_supply + 1;
1.1090 + _res_cap[ra] = -_excess[u];
1.1091 + _cost[a] = 0;
1.1092 + _cost[ra] = 0;
1.1093 + _excess[u] = 0;
1.1094 + }
1.1095 + } else {
1.1096 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.1097 + Value fa = flow[a];
1.1098 + _res_cap[_arc_idf[a]] = cap[a] - fa;
1.1099 + _res_cap[_arc_idb[a]] = fa;
1.1100 + }
1.1101 + for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
1.1102 + int ra = _reverse[a];
1.1103 + _res_cap[a] = 1;
1.1104 + _res_cap[ra] = 0;
1.1105 + _cost[a] = 0;
1.1106 + _cost[ra] = 0;
1.1107 + }
1.1108 + }
1.1109 +
1.1110 + return OPTIMAL;
1.1111 + }
1.1112 +
1.1113 + // Execute the algorithm and transform the results
1.1114 + void start(bool partial_augment) {
1.1115 + // Execute the algorithm
1.1116 + if (partial_augment) {
1.1117 + startPartialAugment();
1.1118 + } else {
1.1119 + startPushRelabel();
1.1120 + }
1.1121 +
1.1122 + // Compute node potentials for the original costs
1.1123 + _arc_vec.clear();
1.1124 + _cost_vec.clear();
1.1125 + for (int j = 0; j != _res_arc_num; ++j) {
1.1126 + if (_res_cap[j] > 0) {
1.1127 + _arc_vec.push_back(IntPair(_source[j], _target[j]));
1.1128 + _cost_vec.push_back(_scost[j]);
1.1129 + }
1.1130 + }
1.1131 + _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1.1132 +
1.1133 + typename BellmanFord<StaticDigraph, LargeCostArcMap>
1.1134 + ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);
1.1135 + bf.distMap(_pi_map);
1.1136 + bf.init(0);
1.1137 + bf.start();
1.1138 +
1.1139 + // Handle non-zero lower bounds
1.1140 + if (_have_lower) {
1.1141 + int limit = _first_out[_root];
1.1142 + for (int j = 0; j != limit; ++j) {
1.1143 + if (!_forward[j]) _res_cap[j] += _lower[j];
1.1144 + }
1.1145 + }
1.1146 }
1.1147
1.1148 /// Execute the algorithm performing partial augmentation and
1.1149 - /// relabel operations.
1.1150 - bool startPartialAugment() {
1.1151 + /// relabel operations
1.1152 + void startPartialAugment() {
1.1153 // Paramters for heuristics
1.1154 -// const int BF_HEURISTIC_EPSILON_BOUND = 1000;
1.1155 -// const int BF_HEURISTIC_BOUND_FACTOR = 3;
1.1156 + const int BF_HEURISTIC_EPSILON_BOUND = 1000;
1.1157 + const int BF_HEURISTIC_BOUND_FACTOR = 3;
1.1158 // Maximum augment path length
1.1159 const int MAX_PATH_LENGTH = 4;
1.1160
1.1161 - // Variables
1.1162 - typename Digraph::template NodeMap<Arc> pred_arc(_graph);
1.1163 - typename Digraph::template NodeMap<bool> forward(_graph);
1.1164 - typename Digraph::template NodeMap<OutArcIt> next_out(_graph);
1.1165 - typename Digraph::template NodeMap<InArcIt> next_in(_graph);
1.1166 - typename Digraph::template NodeMap<bool> next_dir(_graph);
1.1167 - std::deque<Node> active_nodes;
1.1168 - std::vector<Node> path_nodes;
1.1169 -
1.1170 -// int node_num = countNodes(_graph);
1.1171 + // Perform cost scaling phases
1.1172 + IntVector pred_arc(_res_node_num);
1.1173 + std::vector<int> path_nodes;
1.1174 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.1175 1 : _epsilon / _alpha )
1.1176 {
1.1177 -/*
1.1178 // "Early Termination" heuristic: use Bellman-Ford algorithm
1.1179 // to check if the current flow is optimal
1.1180 if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
1.1181 - typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
1.1182 - ShiftCostMap shift_cost(_res_cost, 1);
1.1183 - BellmanFord<ResDigraph, ShiftCostMap> bf(*_res_graph, shift_cost);
1.1184 + _arc_vec.clear();
1.1185 + _cost_vec.clear();
1.1186 + for (int j = 0; j != _res_arc_num; ++j) {
1.1187 + if (_res_cap[j] > 0) {
1.1188 + _arc_vec.push_back(IntPair(_source[j], _target[j]));
1.1189 + _cost_vec.push_back(_cost[j] + 1);
1.1190 + }
1.1191 + }
1.1192 + _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1.1193 +
1.1194 + BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
1.1195 bf.init(0);
1.1196 bool done = false;
1.1197 - int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
1.1198 + int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
1.1199 for (int i = 0; i < K && !done; ++i)
1.1200 done = bf.processNextWeakRound();
1.1201 if (done) break;
1.1202 }
1.1203 -*/
1.1204 +
1.1205 // Saturate arcs not satisfying the optimality condition
1.1206 - Capacity delta;
1.1207 - for (ArcIt e(_graph); e != INVALID; ++e) {
1.1208 - if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.1209 - delta = _capacity[e] - (*_flow)[e];
1.1210 - _excess[_graph.source(e)] -= delta;
1.1211 - _excess[_graph.target(e)] += delta;
1.1212 - (*_flow)[e] = _capacity[e];
1.1213 - }
1.1214 - if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
1.1215 - _excess[_graph.target(e)] -= (*_flow)[e];
1.1216 - _excess[_graph.source(e)] += (*_flow)[e];
1.1217 - (*_flow)[e] = 0;
1.1218 + for (int a = 0; a != _res_arc_num; ++a) {
1.1219 + if (_res_cap[a] > 0 &&
1.1220 + _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
1.1221 + Value delta = _res_cap[a];
1.1222 + _excess[_source[a]] -= delta;
1.1223 + _excess[_target[a]] += delta;
1.1224 + _res_cap[a] = 0;
1.1225 + _res_cap[_reverse[a]] += delta;
1.1226 }
1.1227 }
1.1228 -
1.1229 +
1.1230 // Find active nodes (i.e. nodes with positive excess)
1.1231 - for (NodeIt n(_graph); n != INVALID; ++n) {
1.1232 - if (_excess[n] > 0) active_nodes.push_back(n);
1.1233 + for (int u = 0; u != _res_node_num; ++u) {
1.1234 + if (_excess[u] > 0) _active_nodes.push_back(u);
1.1235 }
1.1236
1.1237 - // Initialize the next arc maps
1.1238 - for (NodeIt n(_graph); n != INVALID; ++n) {
1.1239 - next_out[n] = OutArcIt(_graph, n);
1.1240 - next_in[n] = InArcIt(_graph, n);
1.1241 - next_dir[n] = true;
1.1242 + // Initialize the next arcs
1.1243 + for (int u = 0; u != _res_node_num; ++u) {
1.1244 + _next_out[u] = _first_out[u];
1.1245 }
1.1246
1.1247 // Perform partial augment and relabel operations
1.1248 - while (active_nodes.size() > 0) {
1.1249 + while (true) {
1.1250 // Select an active node (FIFO selection)
1.1251 - if (_excess[active_nodes[0]] <= 0) {
1.1252 - active_nodes.pop_front();
1.1253 - continue;
1.1254 + while (_active_nodes.size() > 0 &&
1.1255 + _excess[_active_nodes.front()] <= 0) {
1.1256 + _active_nodes.pop_front();
1.1257 }
1.1258 - Node start = active_nodes[0];
1.1259 + if (_active_nodes.size() == 0) break;
1.1260 + int start = _active_nodes.front();
1.1261 path_nodes.clear();
1.1262 path_nodes.push_back(start);
1.1263
1.1264 // Find an augmenting path from the start node
1.1265 - Node u, tip = start;
1.1266 - LCost min_red_cost;
1.1267 - while ( _excess[tip] >= 0 &&
1.1268 - int(path_nodes.size()) <= MAX_PATH_LENGTH )
1.1269 - {
1.1270 - if (next_dir[tip]) {
1.1271 - for (OutArcIt e = next_out[tip]; e != INVALID; ++e) {
1.1272 - if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.1273 - u = _graph.target(e);
1.1274 - pred_arc[u] = e;
1.1275 - forward[u] = true;
1.1276 - next_out[tip] = e;
1.1277 - tip = u;
1.1278 - path_nodes.push_back(tip);
1.1279 - goto next_step;
1.1280 - }
1.1281 - }
1.1282 - next_dir[tip] = false;
1.1283 - }
1.1284 - for (InArcIt e = next_in[tip]; e != INVALID; ++e) {
1.1285 - if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
1.1286 - u = _graph.source(e);
1.1287 - pred_arc[u] = e;
1.1288 - forward[u] = false;
1.1289 - next_in[tip] = e;
1.1290 + int tip = start;
1.1291 + while (_excess[tip] >= 0 &&
1.1292 + int(path_nodes.size()) <= MAX_PATH_LENGTH) {
1.1293 + int u;
1.1294 + LargeCost min_red_cost, rc;
1.1295 + int last_out = _sum_supply < 0 ?
1.1296 + _first_out[tip+1] : _first_out[tip+1] - 1;
1.1297 + for (int a = _next_out[tip]; a != last_out; ++a) {
1.1298 + if (_res_cap[a] > 0 &&
1.1299 + _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
1.1300 + u = _target[a];
1.1301 + pred_arc[u] = a;
1.1302 + _next_out[tip] = a;
1.1303 tip = u;
1.1304 path_nodes.push_back(tip);
1.1305 goto next_step;
1.1306 @@ -579,266 +943,186 @@
1.1307 }
1.1308
1.1309 // Relabel tip node
1.1310 - min_red_cost = std::numeric_limits<LCost>::max() / 2;
1.1311 - for (OutArcIt oe(_graph, tip); oe != INVALID; ++oe) {
1.1312 - if ( _capacity[oe] - (*_flow)[oe] > 0 &&
1.1313 - (*_red_cost)[oe] < min_red_cost )
1.1314 - min_red_cost = (*_red_cost)[oe];
1.1315 + min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
1.1316 + for (int a = _first_out[tip]; a != last_out; ++a) {
1.1317 + rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
1.1318 + if (_res_cap[a] > 0 && rc < min_red_cost) {
1.1319 + min_red_cost = rc;
1.1320 + }
1.1321 }
1.1322 - for (InArcIt ie(_graph, tip); ie != INVALID; ++ie) {
1.1323 - if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
1.1324 - min_red_cost = -(*_red_cost)[ie];
1.1325 - }
1.1326 - (*_potential)[tip] -= min_red_cost + _epsilon;
1.1327 + _pi[tip] -= min_red_cost + _epsilon;
1.1328
1.1329 - // Reset the next arc maps
1.1330 - next_out[tip] = OutArcIt(_graph, tip);
1.1331 - next_in[tip] = InArcIt(_graph, tip);
1.1332 - next_dir[tip] = true;
1.1333 + // Reset the next arc of tip
1.1334 + _next_out[tip] = _first_out[tip];
1.1335
1.1336 // Step back
1.1337 if (tip != start) {
1.1338 path_nodes.pop_back();
1.1339 - tip = path_nodes[path_nodes.size()-1];
1.1340 + tip = path_nodes.back();
1.1341 }
1.1342
1.1343 - next_step:
1.1344 - continue;
1.1345 + next_step: ;
1.1346 }
1.1347
1.1348 // Augment along the found path (as much flow as possible)
1.1349 - Capacity delta;
1.1350 + Value delta;
1.1351 + int u, v = path_nodes.front(), pa;
1.1352 for (int i = 1; i < int(path_nodes.size()); ++i) {
1.1353 - u = path_nodes[i];
1.1354 - delta = forward[u] ?
1.1355 - _capacity[pred_arc[u]] - (*_flow)[pred_arc[u]] :
1.1356 - (*_flow)[pred_arc[u]];
1.1357 - delta = std::min(delta, _excess[path_nodes[i-1]]);
1.1358 - (*_flow)[pred_arc[u]] += forward[u] ? delta : -delta;
1.1359 - _excess[path_nodes[i-1]] -= delta;
1.1360 - _excess[u] += delta;
1.1361 - if (_excess[u] > 0 && _excess[u] <= delta) active_nodes.push_back(u);
1.1362 + u = v;
1.1363 + v = path_nodes[i];
1.1364 + pa = pred_arc[v];
1.1365 + delta = std::min(_res_cap[pa], _excess[u]);
1.1366 + _res_cap[pa] -= delta;
1.1367 + _res_cap[_reverse[pa]] += delta;
1.1368 + _excess[u] -= delta;
1.1369 + _excess[v] += delta;
1.1370 + if (_excess[v] > 0 && _excess[v] <= delta)
1.1371 + _active_nodes.push_back(v);
1.1372 }
1.1373 }
1.1374 }
1.1375 -
1.1376 - // Compute node potentials for the original costs
1.1377 - ResidualCostMap<CostMap> res_cost(_orig_cost);
1.1378 - BellmanFord< ResDigraph, ResidualCostMap<CostMap> >
1.1379 - bf(*_res_graph, res_cost);
1.1380 - bf.init(0); bf.start();
1.1381 - for (NodeIt n(_graph); n != INVALID; ++n)
1.1382 - (*_potential)[n] = bf.dist(n);
1.1383 -
1.1384 - // Handle non-zero lower bounds
1.1385 - if (_lower) {
1.1386 - for (ArcIt e(_graph); e != INVALID; ++e)
1.1387 - (*_flow)[e] += (*_lower)[e];
1.1388 - }
1.1389 - return true;
1.1390 }
1.1391
1.1392 - /// Execute the algorithm performing push and relabel operations.
1.1393 - bool startPushRelabel() {
1.1394 + /// Execute the algorithm performing push and relabel operations
1.1395 + void startPushRelabel() {
1.1396 // Paramters for heuristics
1.1397 -// const int BF_HEURISTIC_EPSILON_BOUND = 1000;
1.1398 -// const int BF_HEURISTIC_BOUND_FACTOR = 3;
1.1399 + const int BF_HEURISTIC_EPSILON_BOUND = 1000;
1.1400 + const int BF_HEURISTIC_BOUND_FACTOR = 3;
1.1401
1.1402 - typename Digraph::template NodeMap<bool> hyper(_graph, false);
1.1403 - typename Digraph::template NodeMap<Arc> pred_arc(_graph);
1.1404 - typename Digraph::template NodeMap<bool> forward(_graph);
1.1405 - typename Digraph::template NodeMap<OutArcIt> next_out(_graph);
1.1406 - typename Digraph::template NodeMap<InArcIt> next_in(_graph);
1.1407 - typename Digraph::template NodeMap<bool> next_dir(_graph);
1.1408 - std::deque<Node> active_nodes;
1.1409 -
1.1410 -// int node_num = countNodes(_graph);
1.1411 + // Perform cost scaling phases
1.1412 + BoolVector hyper(_res_node_num, false);
1.1413 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.1414 1 : _epsilon / _alpha )
1.1415 {
1.1416 -/*
1.1417 // "Early Termination" heuristic: use Bellman-Ford algorithm
1.1418 // to check if the current flow is optimal
1.1419 if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
1.1420 - typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
1.1421 - ShiftCostMap shift_cost(_res_cost, 1);
1.1422 - BellmanFord<ResDigraph, ShiftCostMap> bf(*_res_graph, shift_cost);
1.1423 + _arc_vec.clear();
1.1424 + _cost_vec.clear();
1.1425 + for (int j = 0; j != _res_arc_num; ++j) {
1.1426 + if (_res_cap[j] > 0) {
1.1427 + _arc_vec.push_back(IntPair(_source[j], _target[j]));
1.1428 + _cost_vec.push_back(_cost[j] + 1);
1.1429 + }
1.1430 + }
1.1431 + _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1.1432 +
1.1433 + BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
1.1434 bf.init(0);
1.1435 bool done = false;
1.1436 - int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
1.1437 + int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
1.1438 for (int i = 0; i < K && !done; ++i)
1.1439 done = bf.processNextWeakRound();
1.1440 if (done) break;
1.1441 }
1.1442 -*/
1.1443
1.1444 // Saturate arcs not satisfying the optimality condition
1.1445 - Capacity delta;
1.1446 - for (ArcIt e(_graph); e != INVALID; ++e) {
1.1447 - if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.1448 - delta = _capacity[e] - (*_flow)[e];
1.1449 - _excess[_graph.source(e)] -= delta;
1.1450 - _excess[_graph.target(e)] += delta;
1.1451 - (*_flow)[e] = _capacity[e];
1.1452 - }
1.1453 - if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
1.1454 - _excess[_graph.target(e)] -= (*_flow)[e];
1.1455 - _excess[_graph.source(e)] += (*_flow)[e];
1.1456 - (*_flow)[e] = 0;
1.1457 + for (int a = 0; a != _res_arc_num; ++a) {
1.1458 + if (_res_cap[a] > 0 &&
1.1459 + _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
1.1460 + Value delta = _res_cap[a];
1.1461 + _excess[_source[a]] -= delta;
1.1462 + _excess[_target[a]] += delta;
1.1463 + _res_cap[a] = 0;
1.1464 + _res_cap[_reverse[a]] += delta;
1.1465 }
1.1466 }
1.1467
1.1468 // Find active nodes (i.e. nodes with positive excess)
1.1469 - for (NodeIt n(_graph); n != INVALID; ++n) {
1.1470 - if (_excess[n] > 0) active_nodes.push_back(n);
1.1471 + for (int u = 0; u != _res_node_num; ++u) {
1.1472 + if (_excess[u] > 0) _active_nodes.push_back(u);
1.1473 }
1.1474
1.1475 - // Initialize the next arc maps
1.1476 - for (NodeIt n(_graph); n != INVALID; ++n) {
1.1477 - next_out[n] = OutArcIt(_graph, n);
1.1478 - next_in[n] = InArcIt(_graph, n);
1.1479 - next_dir[n] = true;
1.1480 + // Initialize the next arcs
1.1481 + for (int u = 0; u != _res_node_num; ++u) {
1.1482 + _next_out[u] = _first_out[u];
1.1483 }
1.1484
1.1485 // Perform push and relabel operations
1.1486 - while (active_nodes.size() > 0) {
1.1487 + while (_active_nodes.size() > 0) {
1.1488 + LargeCost min_red_cost, rc;
1.1489 + Value delta;
1.1490 + int n, t, a, last_out = _res_arc_num;
1.1491 +
1.1492 // Select an active node (FIFO selection)
1.1493 - Node n = active_nodes[0], t;
1.1494 - bool relabel_enabled = true;
1.1495 + next_node:
1.1496 + n = _active_nodes.front();
1.1497 + last_out = _sum_supply < 0 ?
1.1498 + _first_out[n+1] : _first_out[n+1] - 1;
1.1499
1.1500 // Perform push operations if there are admissible arcs
1.1501 - if (_excess[n] > 0 && next_dir[n]) {
1.1502 - OutArcIt e = next_out[n];
1.1503 - for ( ; e != INVALID; ++e) {
1.1504 - if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.1505 - delta = std::min(_capacity[e] - (*_flow)[e], _excess[n]);
1.1506 - t = _graph.target(e);
1.1507 + if (_excess[n] > 0) {
1.1508 + for (a = _next_out[n]; a != last_out; ++a) {
1.1509 + if (_res_cap[a] > 0 &&
1.1510 + _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
1.1511 + delta = std::min(_res_cap[a], _excess[n]);
1.1512 + t = _target[a];
1.1513
1.1514 // Push-look-ahead heuristic
1.1515 - Capacity ahead = -_excess[t];
1.1516 - for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) {
1.1517 - if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
1.1518 - ahead += _capacity[oe] - (*_flow)[oe];
1.1519 - }
1.1520 - for (InArcIt ie(_graph, t); ie != INVALID; ++ie) {
1.1521 - if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
1.1522 - ahead += (*_flow)[ie];
1.1523 + Value ahead = -_excess[t];
1.1524 + int last_out_t = _sum_supply < 0 ?
1.1525 + _first_out[t+1] : _first_out[t+1] - 1;
1.1526 + for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
1.1527 + if (_res_cap[ta] > 0 &&
1.1528 + _cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0)
1.1529 + ahead += _res_cap[ta];
1.1530 + if (ahead >= delta) break;
1.1531 }
1.1532 if (ahead < 0) ahead = 0;
1.1533
1.1534 // Push flow along the arc
1.1535 if (ahead < delta) {
1.1536 - (*_flow)[e] += ahead;
1.1537 + _res_cap[a] -= ahead;
1.1538 + _res_cap[_reverse[a]] += ahead;
1.1539 _excess[n] -= ahead;
1.1540 _excess[t] += ahead;
1.1541 - active_nodes.push_front(t);
1.1542 + _active_nodes.push_front(t);
1.1543 hyper[t] = true;
1.1544 - relabel_enabled = false;
1.1545 - break;
1.1546 + _next_out[n] = a;
1.1547 + goto next_node;
1.1548 } else {
1.1549 - (*_flow)[e] += delta;
1.1550 + _res_cap[a] -= delta;
1.1551 + _res_cap[_reverse[a]] += delta;
1.1552 _excess[n] -= delta;
1.1553 _excess[t] += delta;
1.1554 if (_excess[t] > 0 && _excess[t] <= delta)
1.1555 - active_nodes.push_back(t);
1.1556 + _active_nodes.push_back(t);
1.1557 }
1.1558
1.1559 - if (_excess[n] == 0) break;
1.1560 + if (_excess[n] == 0) {
1.1561 + _next_out[n] = a;
1.1562 + goto remove_nodes;
1.1563 + }
1.1564 }
1.1565 }
1.1566 - if (e != INVALID) {
1.1567 - next_out[n] = e;
1.1568 - } else {
1.1569 - next_dir[n] = false;
1.1570 - }
1.1571 - }
1.1572 -
1.1573 - if (_excess[n] > 0 && !next_dir[n]) {
1.1574 - InArcIt e = next_in[n];
1.1575 - for ( ; e != INVALID; ++e) {
1.1576 - if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
1.1577 - delta = std::min((*_flow)[e], _excess[n]);
1.1578 - t = _graph.source(e);
1.1579 -
1.1580 - // Push-look-ahead heuristic
1.1581 - Capacity ahead = -_excess[t];
1.1582 - for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) {
1.1583 - if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
1.1584 - ahead += _capacity[oe] - (*_flow)[oe];
1.1585 - }
1.1586 - for (InArcIt ie(_graph, t); ie != INVALID; ++ie) {
1.1587 - if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
1.1588 - ahead += (*_flow)[ie];
1.1589 - }
1.1590 - if (ahead < 0) ahead = 0;
1.1591 -
1.1592 - // Push flow along the arc
1.1593 - if (ahead < delta) {
1.1594 - (*_flow)[e] -= ahead;
1.1595 - _excess[n] -= ahead;
1.1596 - _excess[t] += ahead;
1.1597 - active_nodes.push_front(t);
1.1598 - hyper[t] = true;
1.1599 - relabel_enabled = false;
1.1600 - break;
1.1601 - } else {
1.1602 - (*_flow)[e] -= delta;
1.1603 - _excess[n] -= delta;
1.1604 - _excess[t] += delta;
1.1605 - if (_excess[t] > 0 && _excess[t] <= delta)
1.1606 - active_nodes.push_back(t);
1.1607 - }
1.1608 -
1.1609 - if (_excess[n] == 0) break;
1.1610 - }
1.1611 - }
1.1612 - next_in[n] = e;
1.1613 + _next_out[n] = a;
1.1614 }
1.1615
1.1616 // Relabel the node if it is still active (or hyper)
1.1617 - if (relabel_enabled && (_excess[n] > 0 || hyper[n])) {
1.1618 - LCost min_red_cost = std::numeric_limits<LCost>::max() / 2;
1.1619 - for (OutArcIt oe(_graph, n); oe != INVALID; ++oe) {
1.1620 - if ( _capacity[oe] - (*_flow)[oe] > 0 &&
1.1621 - (*_red_cost)[oe] < min_red_cost )
1.1622 - min_red_cost = (*_red_cost)[oe];
1.1623 + if (_excess[n] > 0 || hyper[n]) {
1.1624 + min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
1.1625 + for (int a = _first_out[n]; a != last_out; ++a) {
1.1626 + rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
1.1627 + if (_res_cap[a] > 0 && rc < min_red_cost) {
1.1628 + min_red_cost = rc;
1.1629 + }
1.1630 }
1.1631 - for (InArcIt ie(_graph, n); ie != INVALID; ++ie) {
1.1632 - if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
1.1633 - min_red_cost = -(*_red_cost)[ie];
1.1634 - }
1.1635 - (*_potential)[n] -= min_red_cost + _epsilon;
1.1636 + _pi[n] -= min_red_cost + _epsilon;
1.1637 hyper[n] = false;
1.1638
1.1639 - // Reset the next arc maps
1.1640 - next_out[n] = OutArcIt(_graph, n);
1.1641 - next_in[n] = InArcIt(_graph, n);
1.1642 - next_dir[n] = true;
1.1643 + // Reset the next arc
1.1644 + _next_out[n] = _first_out[n];
1.1645 }
1.1646 -
1.1647 +
1.1648 // Remove nodes that are not active nor hyper
1.1649 - while ( active_nodes.size() > 0 &&
1.1650 - _excess[active_nodes[0]] <= 0 &&
1.1651 - !hyper[active_nodes[0]] ) {
1.1652 - active_nodes.pop_front();
1.1653 + remove_nodes:
1.1654 + while ( _active_nodes.size() > 0 &&
1.1655 + _excess[_active_nodes.front()] <= 0 &&
1.1656 + !hyper[_active_nodes.front()] ) {
1.1657 + _active_nodes.pop_front();
1.1658 }
1.1659 }
1.1660 }
1.1661 -
1.1662 - // Compute node potentials for the original costs
1.1663 - ResidualCostMap<CostMap> res_cost(_orig_cost);
1.1664 - BellmanFord< ResDigraph, ResidualCostMap<CostMap> >
1.1665 - bf(*_res_graph, res_cost);
1.1666 - bf.init(0); bf.start();
1.1667 - for (NodeIt n(_graph); n != INVALID; ++n)
1.1668 - (*_potential)[n] = bf.dist(n);
1.1669 -
1.1670 - // Handle non-zero lower bounds
1.1671 - if (_lower) {
1.1672 - for (ArcIt e(_graph); e != INVALID; ++e)
1.1673 - (*_flow)[e] += (*_lower)[e];
1.1674 - }
1.1675 - return true;
1.1676 }
1.1677
1.1678 }; //class CostScaling