kpeter@2577: /* -*- C++ -*- kpeter@2577: * kpeter@2577: * This file is a part of LEMON, a generic C++ optimization library kpeter@2577: * kpeter@2577: * Copyright (C) 2003-2008 kpeter@2577: * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport kpeter@2577: * (Egervary Research Group on Combinatorial Optimization, EGRES). kpeter@2577: * kpeter@2577: * Permission to use, modify and distribute this software is granted kpeter@2577: * provided that this copyright notice appears in all copies. For kpeter@2577: * precise terms see the accompanying LICENSE file. kpeter@2577: * kpeter@2577: * This software is provided "AS IS" with no warranty of any kind, kpeter@2577: * express or implied, and with no claim as to its suitability for any kpeter@2577: * purpose. kpeter@2577: * kpeter@2577: */ kpeter@2577: kpeter@2577: #ifndef LEMON_COST_SCALING_H kpeter@2577: #define LEMON_COST_SCALING_H kpeter@2577: kpeter@2577: /// \ingroup min_cost_flow kpeter@2577: /// kpeter@2577: /// \file kpeter@2577: /// \brief Cost scaling algorithm for finding a minimum cost flow. kpeter@2577: kpeter@2577: #include kpeter@2577: #include kpeter@2577: #include kpeter@2577: #include kpeter@2577: #include kpeter@2577: kpeter@2577: #include kpeter@2577: #include kpeter@2577: kpeter@2577: namespace lemon { kpeter@2577: kpeter@2577: /// \addtogroup min_cost_flow kpeter@2577: /// @{ kpeter@2577: kpeter@2577: /// \brief Implementation of the cost scaling algorithm for finding a kpeter@2577: /// minimum cost flow. kpeter@2577: /// kpeter@2577: /// \ref CostScaling implements the cost scaling algorithm performing kpeter@2577: /// generalized push-relabel operations for finding a minimum cost kpeter@2577: /// flow. kpeter@2577: /// kpeter@2577: /// \tparam Graph The directed graph type the algorithm runs on. kpeter@2577: /// \tparam LowerMap The type of the lower bound map. kpeter@2577: /// \tparam CapacityMap The type of the capacity (upper bound) map. kpeter@2577: /// \tparam CostMap The type of the cost (length) map. kpeter@2577: /// \tparam SupplyMap The type of the supply map. kpeter@2577: /// kpeter@2577: /// \warning kpeter@2577: /// - Edge capacities and costs should be \e non-negative \e integers. kpeter@2577: /// - Supply values should be \e signed \e integers. kpeter@2581: /// - The value types of the maps should be convertible to each other. kpeter@2581: /// - \c CostMap::Value must be signed type. kpeter@2577: /// kpeter@2577: /// \note Edge costs are multiplied with the number of nodes during kpeter@2577: /// the algorithm so overflow problems may arise more easily than with kpeter@2577: /// other minimum cost flow algorithms. kpeter@2577: /// If it is available, long long int type is used instead of kpeter@2577: /// long int in the inside computations. kpeter@2577: /// kpeter@2577: /// \author Peter Kovacs kpeter@2577: kpeter@2577: template < typename Graph, kpeter@2577: typename LowerMap = typename Graph::template EdgeMap, kpeter@2577: typename CapacityMap = typename Graph::template EdgeMap, kpeter@2577: typename CostMap = typename Graph::template EdgeMap, kpeter@2577: typename SupplyMap = typename Graph::template NodeMap > kpeter@2577: class CostScaling kpeter@2577: { kpeter@2577: GRAPH_TYPEDEFS(typename Graph); kpeter@2577: kpeter@2577: typedef typename CapacityMap::Value Capacity; kpeter@2577: typedef typename CostMap::Value Cost; kpeter@2577: typedef typename SupplyMap::Value Supply; kpeter@2577: typedef typename Graph::template EdgeMap CapacityEdgeMap; kpeter@2577: typedef typename Graph::template NodeMap SupplyNodeMap; kpeter@2577: kpeter@2577: typedef ResGraphAdaptor< const Graph, Capacity, kpeter@2577: CapacityEdgeMap, CapacityEdgeMap > ResGraph; kpeter@2577: typedef typename ResGraph::Edge ResEdge; kpeter@2577: kpeter@2577: #if defined __GNUC__ && !defined __STRICT_ANSI__ kpeter@2577: typedef long long int LCost; kpeter@2577: #else kpeter@2577: typedef long int LCost; kpeter@2577: #endif kpeter@2577: typedef typename Graph::template EdgeMap LargeCostMap; kpeter@2577: kpeter@2577: public: kpeter@2577: kpeter@2577: /// The type of the flow map. kpeter@2581: typedef typename Graph::template EdgeMap FlowMap; kpeter@2577: /// The type of the potential map. kpeter@2577: typedef typename Graph::template NodeMap PotentialMap; kpeter@2577: kpeter@2577: private: kpeter@2577: kpeter@2577: /// \brief Map adaptor class for handling residual edge costs. kpeter@2577: /// kpeter@2577: /// \ref ResidualCostMap is a map adaptor class for handling kpeter@2577: /// residual edge costs. kpeter@2581: template kpeter@2581: class ResidualCostMap : public MapBase kpeter@2577: { kpeter@2577: private: kpeter@2577: kpeter@2581: const Map &_cost_map; kpeter@2577: kpeter@2577: public: kpeter@2577: kpeter@2577: ///\e kpeter@2581: ResidualCostMap(const Map &cost_map) : kpeter@2577: _cost_map(cost_map) {} kpeter@2577: kpeter@2577: ///\e kpeter@2581: typename Map::Value operator[](const ResEdge &e) const { kpeter@2577: return ResGraph::forward(e) ? _cost_map[e] : -_cost_map[e]; kpeter@2577: } kpeter@2577: kpeter@2577: }; //class ResidualCostMap kpeter@2577: kpeter@2577: /// \brief Map adaptor class for handling reduced edge costs. kpeter@2577: /// kpeter@2577: /// \ref ReducedCostMap is a map adaptor class for handling reduced kpeter@2577: /// edge costs. kpeter@2577: class ReducedCostMap : public MapBase kpeter@2577: { kpeter@2577: private: kpeter@2577: kpeter@2577: const Graph &_gr; kpeter@2577: const LargeCostMap &_cost_map; kpeter@2577: const PotentialMap &_pot_map; kpeter@2577: kpeter@2577: public: kpeter@2577: kpeter@2577: ///\e kpeter@2577: ReducedCostMap( const Graph &gr, kpeter@2577: const LargeCostMap &cost_map, kpeter@2577: const PotentialMap &pot_map ) : kpeter@2577: _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {} kpeter@2577: kpeter@2577: ///\e kpeter@2577: LCost operator[](const Edge &e) const { kpeter@2577: return _cost_map[e] + _pot_map[_gr.source(e)] kpeter@2577: - _pot_map[_gr.target(e)]; kpeter@2577: } kpeter@2577: kpeter@2577: }; //class ReducedCostMap kpeter@2577: kpeter@2577: private: kpeter@2577: kpeter@2577: // Scaling factor kpeter@2577: static const int ALPHA = 4; kpeter@2577: kpeter@2577: // Paramters for heuristics kpeter@2581: static const int BF_HEURISTIC_EPSILON_BOUND = 5000; kpeter@2581: static const int BF_HEURISTIC_BOUND_FACTOR = 3; kpeter@2577: kpeter@2577: private: kpeter@2577: kpeter@2577: // The directed graph the algorithm runs on kpeter@2577: const Graph &_graph; kpeter@2577: // The original lower bound map kpeter@2577: const LowerMap *_lower; kpeter@2577: // The modified capacity map kpeter@2577: CapacityEdgeMap _capacity; kpeter@2577: // The original cost map kpeter@2577: const CostMap &_orig_cost; kpeter@2577: // The scaled cost map kpeter@2577: LargeCostMap _cost; kpeter@2577: // The modified supply map kpeter@2577: SupplyNodeMap _supply; kpeter@2577: bool _valid_supply; kpeter@2577: kpeter@2577: // Edge map of the current flow kpeter@2581: FlowMap *_flow; kpeter@2581: bool _local_flow; kpeter@2577: // Node map of the current potentials kpeter@2581: PotentialMap *_potential; kpeter@2581: bool _local_potential; kpeter@2577: kpeter@2577: // The residual graph kpeter@2581: ResGraph *_res_graph; kpeter@2577: // The residual cost map kpeter@2581: ResidualCostMap _res_cost; kpeter@2577: // The reduced cost map kpeter@2581: ReducedCostMap *_red_cost; kpeter@2577: // The excess map kpeter@2577: SupplyNodeMap _excess; kpeter@2577: // The epsilon parameter used for cost scaling kpeter@2577: LCost _epsilon; kpeter@2577: kpeter@2577: public: kpeter@2577: kpeter@2581: /// \brief General constructor (with lower bounds). kpeter@2577: /// kpeter@2581: /// General constructor (with lower bounds). kpeter@2577: /// kpeter@2577: /// \param graph The directed graph the algorithm runs on. kpeter@2577: /// \param lower The lower bounds of the edges. kpeter@2577: /// \param capacity The capacities (upper bounds) of the edges. kpeter@2577: /// \param cost The cost (length) values of the edges. kpeter@2577: /// \param supply The supply values of the nodes (signed). kpeter@2577: CostScaling( const Graph &graph, kpeter@2577: const LowerMap &lower, kpeter@2577: const CapacityMap &capacity, kpeter@2577: const CostMap &cost, kpeter@2577: const SupplyMap &supply ) : kpeter@2577: _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost), kpeter@2581: _cost(graph), _supply(graph), _flow(0), _local_flow(false), kpeter@2581: _potential(0), _local_potential(false), _res_cost(_cost), kpeter@2581: _excess(graph, 0) kpeter@2577: { kpeter@2577: // Removing non-zero lower bounds kpeter@2577: _capacity = subMap(capacity, lower); kpeter@2577: Supply sum = 0; kpeter@2577: for (NodeIt n(_graph); n != INVALID; ++n) { kpeter@2577: Supply s = supply[n]; kpeter@2577: for (InEdgeIt e(_graph, n); e != INVALID; ++e) kpeter@2577: s += lower[e]; kpeter@2577: for (OutEdgeIt e(_graph, n); e != INVALID; ++e) kpeter@2577: s -= lower[e]; kpeter@2577: _supply[n] = s; kpeter@2577: sum += s; kpeter@2577: } kpeter@2577: _valid_supply = sum == 0; kpeter@2577: } kpeter@2577: kpeter@2581: /// \brief General constructor (without lower bounds). kpeter@2577: /// kpeter@2581: /// General constructor (without lower bounds). kpeter@2577: /// kpeter@2577: /// \param graph The directed graph the algorithm runs on. kpeter@2577: /// \param capacity The capacities (upper bounds) of the edges. kpeter@2577: /// \param cost The cost (length) values of the edges. kpeter@2577: /// \param supply The supply values of the nodes (signed). kpeter@2577: CostScaling( const Graph &graph, kpeter@2577: const CapacityMap &capacity, kpeter@2577: const CostMap &cost, kpeter@2577: const SupplyMap &supply ) : kpeter@2577: _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost), kpeter@2581: _cost(graph), _supply(supply), _flow(0), _local_flow(false), kpeter@2581: _potential(0), _local_potential(false), _res_cost(_cost), kpeter@2581: _excess(graph, 0) kpeter@2577: { kpeter@2577: // Checking the sum of supply values kpeter@2577: Supply sum = 0; kpeter@2577: for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n]; kpeter@2577: _valid_supply = sum == 0; kpeter@2577: } kpeter@2577: kpeter@2581: /// \brief Simple constructor (with lower bounds). kpeter@2577: /// kpeter@2581: /// Simple constructor (with lower bounds). kpeter@2577: /// kpeter@2577: /// \param graph The directed graph the algorithm runs on. kpeter@2577: /// \param lower The lower bounds of the edges. kpeter@2577: /// \param capacity The capacities (upper bounds) of the edges. kpeter@2577: /// \param cost The cost (length) values of the edges. kpeter@2577: /// \param s The source node. kpeter@2577: /// \param t The target node. kpeter@2577: /// \param flow_value The required amount of flow from node \c s kpeter@2577: /// to node \c t (i.e. the supply of \c s and the demand of \c t). kpeter@2577: CostScaling( const Graph &graph, kpeter@2577: const LowerMap &lower, kpeter@2577: const CapacityMap &capacity, kpeter@2577: const CostMap &cost, kpeter@2577: Node s, Node t, kpeter@2577: Supply flow_value ) : kpeter@2577: _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost), kpeter@2581: _cost(graph), _supply(graph), _flow(0), _local_flow(false), kpeter@2581: _potential(0), _local_potential(false), _res_cost(_cost), kpeter@2581: _excess(graph, 0) kpeter@2577: { kpeter@2577: // Removing nonzero lower bounds kpeter@2577: _capacity = subMap(capacity, lower); kpeter@2577: for (NodeIt n(_graph); n != INVALID; ++n) { kpeter@2577: Supply sum = 0; kpeter@2577: if (n == s) sum = flow_value; kpeter@2577: if (n == t) sum = -flow_value; kpeter@2577: for (InEdgeIt e(_graph, n); e != INVALID; ++e) kpeter@2577: sum += lower[e]; kpeter@2577: for (OutEdgeIt e(_graph, n); e != INVALID; ++e) kpeter@2577: sum -= lower[e]; kpeter@2577: _supply[n] = sum; kpeter@2577: } kpeter@2577: _valid_supply = true; kpeter@2577: } kpeter@2577: kpeter@2581: /// \brief Simple constructor (without lower bounds). kpeter@2577: /// kpeter@2581: /// Simple constructor (without lower bounds). kpeter@2577: /// kpeter@2577: /// \param graph The directed graph the algorithm runs on. kpeter@2577: /// \param capacity The capacities (upper bounds) of the edges. kpeter@2577: /// \param cost The cost (length) values of the edges. kpeter@2577: /// \param s The source node. kpeter@2577: /// \param t The target node. kpeter@2577: /// \param flow_value The required amount of flow from node \c s kpeter@2577: /// to node \c t (i.e. the supply of \c s and the demand of \c t). kpeter@2577: CostScaling( const Graph &graph, kpeter@2577: const CapacityMap &capacity, kpeter@2577: const CostMap &cost, kpeter@2577: Node s, Node t, kpeter@2577: Supply flow_value ) : kpeter@2577: _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost), kpeter@2581: _cost(graph), _supply(graph, 0), _flow(0), _local_flow(false), kpeter@2581: _potential(0), _local_potential(false), _res_cost(_cost), kpeter@2581: _excess(graph, 0) kpeter@2577: { kpeter@2577: _supply[s] = flow_value; kpeter@2577: _supply[t] = -flow_value; kpeter@2577: _valid_supply = true; kpeter@2577: } kpeter@2577: kpeter@2581: /// Destructor. kpeter@2581: ~CostScaling() { kpeter@2581: if (_local_flow) delete _flow; kpeter@2581: if (_local_potential) delete _potential; kpeter@2581: delete _res_graph; kpeter@2581: delete _red_cost; kpeter@2581: } kpeter@2581: kpeter@2581: /// \brief Sets the flow map. kpeter@2581: /// kpeter@2581: /// Sets the flow map. kpeter@2581: /// kpeter@2581: /// \return \c (*this) kpeter@2581: CostScaling& flowMap(FlowMap &map) { kpeter@2581: if (_local_flow) { kpeter@2581: delete _flow; kpeter@2581: _local_flow = false; kpeter@2581: } kpeter@2581: _flow = ↦ kpeter@2581: return *this; kpeter@2581: } kpeter@2581: kpeter@2581: /// \brief Sets the potential map. kpeter@2581: /// kpeter@2581: /// Sets the potential map. kpeter@2581: /// kpeter@2581: /// \return \c (*this) kpeter@2581: CostScaling& potentialMap(PotentialMap &map) { kpeter@2581: if (_local_potential) { kpeter@2581: delete _potential; kpeter@2581: _local_potential = false; kpeter@2581: } kpeter@2581: _potential = ↦ kpeter@2581: return *this; kpeter@2581: } kpeter@2581: kpeter@2581: /// \name Execution control kpeter@2581: /// The only way to execute the algorithm is to call the run() kpeter@2581: /// function. kpeter@2581: kpeter@2581: /// @{ kpeter@2581: kpeter@2577: /// \brief Runs the algorithm. kpeter@2577: /// kpeter@2577: /// Runs the algorithm. kpeter@2577: /// kpeter@2577: /// \return \c true if a feasible flow can be found. kpeter@2577: bool run() { kpeter@2581: return init() && start(); kpeter@2577: } kpeter@2577: kpeter@2581: /// @} kpeter@2581: kpeter@2581: /// \name Query Functions kpeter@2581: /// The result of the algorithm can be obtained using these kpeter@2581: /// functions. kpeter@2581: /// \n run() must be called before using them. kpeter@2581: kpeter@2581: /// @{ kpeter@2581: kpeter@2577: /// \brief Returns a const reference to the edge map storing the kpeter@2577: /// found flow. kpeter@2577: /// kpeter@2577: /// Returns a const reference to the edge map storing the found flow. kpeter@2577: /// kpeter@2577: /// \pre \ref run() must be called before using this function. kpeter@2577: const FlowMap& flowMap() const { kpeter@2581: return *_flow; kpeter@2577: } kpeter@2577: kpeter@2577: /// \brief Returns a const reference to the node map storing the kpeter@2577: /// found potentials (the dual solution). kpeter@2577: /// kpeter@2577: /// Returns a const reference to the node map storing the found kpeter@2577: /// potentials (the dual solution). kpeter@2577: /// kpeter@2577: /// \pre \ref run() must be called before using this function. kpeter@2577: const PotentialMap& potentialMap() const { kpeter@2581: return *_potential; kpeter@2581: } kpeter@2581: kpeter@2581: /// \brief Returns the flow on the edge. kpeter@2581: /// kpeter@2581: /// Returns the flow on the edge. kpeter@2581: /// kpeter@2581: /// \pre \ref run() must be called before using this function. kpeter@2581: Capacity flow(const Edge& edge) const { kpeter@2581: return (*_flow)[edge]; kpeter@2581: } kpeter@2581: kpeter@2581: /// \brief Returns the potential of the node. kpeter@2581: /// kpeter@2581: /// Returns the potential of the node. kpeter@2581: /// kpeter@2581: /// \pre \ref run() must be called before using this function. kpeter@2581: Cost potential(const Node& node) const { kpeter@2581: return (*_potential)[node]; kpeter@2577: } kpeter@2577: kpeter@2577: /// \brief Returns the total cost of the found flow. kpeter@2577: /// kpeter@2577: /// Returns the total cost of the found flow. The complexity of the kpeter@2577: /// function is \f$ O(e) \f$. kpeter@2577: /// kpeter@2577: /// \pre \ref run() must be called before using this function. kpeter@2577: Cost totalCost() const { kpeter@2577: Cost c = 0; kpeter@2577: for (EdgeIt e(_graph); e != INVALID; ++e) kpeter@2581: c += (*_flow)[e] * _orig_cost[e]; kpeter@2577: return c; kpeter@2577: } kpeter@2577: kpeter@2581: /// @} kpeter@2581: kpeter@2577: private: kpeter@2577: kpeter@2577: /// Initializes the algorithm. kpeter@2577: bool init() { kpeter@2577: if (!_valid_supply) return false; kpeter@2577: kpeter@2581: // Initializing flow and potential maps kpeter@2581: if (!_flow) { kpeter@2581: _flow = new FlowMap(_graph); kpeter@2581: _local_flow = true; kpeter@2581: } kpeter@2581: if (!_potential) { kpeter@2581: _potential = new PotentialMap(_graph); kpeter@2581: _local_potential = true; kpeter@2581: } kpeter@2581: kpeter@2581: _red_cost = new ReducedCostMap(_graph, _cost, *_potential); kpeter@2581: _res_graph = new ResGraph(_graph, _capacity, *_flow); kpeter@2581: kpeter@2577: // Initializing the scaled cost map and the epsilon parameter kpeter@2577: Cost max_cost = 0; kpeter@2577: int node_num = countNodes(_graph); kpeter@2577: for (EdgeIt e(_graph); e != INVALID; ++e) { kpeter@2577: _cost[e] = LCost(_orig_cost[e]) * node_num * ALPHA; kpeter@2577: if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e]; kpeter@2577: } kpeter@2577: _epsilon = max_cost * node_num; kpeter@2577: kpeter@2577: // Finding a feasible flow using Circulation kpeter@2577: Circulation< Graph, ConstMap, CapacityEdgeMap, kpeter@2577: SupplyMap > kpeter@2581: circulation( _graph, constMap(Capacity(0)), _capacity, kpeter@2577: _supply ); kpeter@2581: return circulation.flowMap(*_flow).run(); kpeter@2577: } kpeter@2577: kpeter@2577: kpeter@2577: /// Executes the algorithm. kpeter@2577: bool start() { kpeter@2577: std::deque active_nodes; kpeter@2577: typename Graph::template NodeMap hyper(_graph, false); kpeter@2577: kpeter@2577: int node_num = countNodes(_graph); kpeter@2577: for ( ; _epsilon >= 1; _epsilon = _epsilon < ALPHA && _epsilon > 1 ? kpeter@2577: 1 : _epsilon / ALPHA ) kpeter@2577: { kpeter@2577: // Performing price refinement heuristic using Bellman-Ford kpeter@2577: // algorithm kpeter@2577: if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) { kpeter@2581: typedef ShiftMap< ResidualCostMap > ShiftCostMap; kpeter@2577: ShiftCostMap shift_cost(_res_cost, _epsilon); kpeter@2581: BellmanFord bf(*_res_graph, shift_cost); kpeter@2577: bf.init(0); kpeter@2577: bool done = false; kpeter@2577: int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num)); kpeter@2577: for (int i = 0; i < K && !done; ++i) kpeter@2577: done = bf.processNextWeakRound(); kpeter@2577: if (done) { kpeter@2577: for (NodeIt n(_graph); n != INVALID; ++n) kpeter@2581: (*_potential)[n] = bf.dist(n); kpeter@2577: continue; kpeter@2577: } kpeter@2577: } kpeter@2577: kpeter@2577: // Saturating edges not satisfying the optimality condition kpeter@2577: Capacity delta; kpeter@2577: for (EdgeIt e(_graph); e != INVALID; ++e) { kpeter@2581: if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) { kpeter@2581: delta = _capacity[e] - (*_flow)[e]; kpeter@2577: _excess[_graph.source(e)] -= delta; kpeter@2577: _excess[_graph.target(e)] += delta; kpeter@2581: (*_flow)[e] = _capacity[e]; kpeter@2577: } kpeter@2581: if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) { kpeter@2581: _excess[_graph.target(e)] -= (*_flow)[e]; kpeter@2581: _excess[_graph.source(e)] += (*_flow)[e]; kpeter@2581: (*_flow)[e] = 0; kpeter@2577: } kpeter@2577: } kpeter@2577: kpeter@2577: // Finding active nodes (i.e. nodes with positive excess) kpeter@2577: for (NodeIt n(_graph); n != INVALID; ++n) kpeter@2577: if (_excess[n] > 0) active_nodes.push_back(n); kpeter@2577: kpeter@2577: // Performing push and relabel operations kpeter@2577: while (active_nodes.size() > 0) { kpeter@2577: Node n = active_nodes[0], t; kpeter@2577: bool relabel_enabled = true; kpeter@2577: kpeter@2577: // Performing push operations if there are admissible edges kpeter@2577: if (_excess[n] > 0) { kpeter@2577: for (OutEdgeIt e(_graph, n); e != INVALID; ++e) { kpeter@2581: if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) { kpeter@2581: delta = _capacity[e] - (*_flow)[e] <= _excess[n] ? kpeter@2581: _capacity[e] - (*_flow)[e] : _excess[n]; kpeter@2577: t = _graph.target(e); kpeter@2577: kpeter@2577: // Push-look-ahead heuristic kpeter@2577: Capacity ahead = -_excess[t]; kpeter@2577: for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) { kpeter@2581: if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0) kpeter@2581: ahead += _capacity[oe] - (*_flow)[oe]; kpeter@2577: } kpeter@2577: for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) { kpeter@2581: if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0) kpeter@2581: ahead += (*_flow)[ie]; kpeter@2577: } kpeter@2577: if (ahead < 0) ahead = 0; kpeter@2577: kpeter@2577: // Pushing flow along the edge kpeter@2577: if (ahead < delta) { kpeter@2581: (*_flow)[e] += ahead; kpeter@2577: _excess[n] -= ahead; kpeter@2577: _excess[t] += ahead; kpeter@2577: active_nodes.push_front(t); kpeter@2577: hyper[t] = true; kpeter@2577: relabel_enabled = false; kpeter@2577: break; kpeter@2577: } else { kpeter@2581: (*_flow)[e] += delta; kpeter@2577: _excess[n] -= delta; kpeter@2577: _excess[t] += delta; kpeter@2577: if (_excess[t] > 0 && _excess[t] <= delta) kpeter@2577: active_nodes.push_back(t); kpeter@2577: } kpeter@2577: kpeter@2577: if (_excess[n] == 0) break; kpeter@2577: } kpeter@2577: } kpeter@2577: } kpeter@2577: kpeter@2577: if (_excess[n] > 0) { kpeter@2577: for (InEdgeIt e(_graph, n); e != INVALID; ++e) { kpeter@2581: if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) { kpeter@2581: delta = (*_flow)[e] <= _excess[n] ? (*_flow)[e] : _excess[n]; kpeter@2577: t = _graph.source(e); kpeter@2577: kpeter@2577: // Push-look-ahead heuristic kpeter@2577: Capacity ahead = -_excess[t]; kpeter@2577: for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) { kpeter@2581: if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0) kpeter@2581: ahead += _capacity[oe] - (*_flow)[oe]; kpeter@2577: } kpeter@2577: for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) { kpeter@2581: if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0) kpeter@2581: ahead += (*_flow)[ie]; kpeter@2577: } kpeter@2577: if (ahead < 0) ahead = 0; kpeter@2577: kpeter@2577: // Pushing flow along the edge kpeter@2577: if (ahead < delta) { kpeter@2581: (*_flow)[e] -= ahead; kpeter@2577: _excess[n] -= ahead; kpeter@2577: _excess[t] += ahead; kpeter@2577: active_nodes.push_front(t); kpeter@2577: hyper[t] = true; kpeter@2577: relabel_enabled = false; kpeter@2577: break; kpeter@2577: } else { kpeter@2581: (*_flow)[e] -= delta; kpeter@2577: _excess[n] -= delta; kpeter@2577: _excess[t] += delta; kpeter@2577: if (_excess[t] > 0 && _excess[t] <= delta) kpeter@2577: active_nodes.push_back(t); kpeter@2577: } kpeter@2577: kpeter@2577: if (_excess[n] == 0) break; kpeter@2577: } kpeter@2577: } kpeter@2577: } kpeter@2577: kpeter@2577: if (relabel_enabled && (_excess[n] > 0 || hyper[n])) { kpeter@2577: // Performing relabel operation if the node is still active kpeter@2577: LCost min_red_cost = std::numeric_limits::max(); kpeter@2577: for (OutEdgeIt oe(_graph, n); oe != INVALID; ++oe) { kpeter@2581: if ( _capacity[oe] - (*_flow)[oe] > 0 && kpeter@2581: (*_red_cost)[oe] < min_red_cost ) kpeter@2581: min_red_cost = (*_red_cost)[oe]; kpeter@2577: } kpeter@2577: for (InEdgeIt ie(_graph, n); ie != INVALID; ++ie) { kpeter@2581: if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost) kpeter@2581: min_red_cost = -(*_red_cost)[ie]; kpeter@2577: } kpeter@2581: (*_potential)[n] -= min_red_cost + _epsilon; kpeter@2577: hyper[n] = false; kpeter@2577: } kpeter@2577: kpeter@2577: // Removing active nodes with non-positive excess kpeter@2577: while ( active_nodes.size() > 0 && kpeter@2577: _excess[active_nodes[0]] <= 0 && kpeter@2577: !hyper[active_nodes[0]] ) { kpeter@2577: active_nodes.pop_front(); kpeter@2577: } kpeter@2577: } kpeter@2577: } kpeter@2577: kpeter@2581: // Computing node potentials for the original costs kpeter@2581: ResidualCostMap res_cost(_orig_cost); kpeter@2581: BellmanFord< ResGraph, ResidualCostMap > kpeter@2581: bf(*_res_graph, res_cost); kpeter@2581: bf.init(0); bf.start(); kpeter@2581: for (NodeIt n(_graph); n != INVALID; ++n) kpeter@2581: (*_potential)[n] = bf.dist(n); kpeter@2581: kpeter@2577: // Handling non-zero lower bounds kpeter@2577: if (_lower) { kpeter@2577: for (EdgeIt e(_graph); e != INVALID; ++e) kpeter@2581: (*_flow)[e] += (*_lower)[e]; kpeter@2577: } kpeter@2577: return true; kpeter@2577: } kpeter@2577: kpeter@2577: }; //class CostScaling kpeter@2577: kpeter@2577: ///@} kpeter@2577: kpeter@2577: } //namespace lemon kpeter@2577: kpeter@2577: #endif //LEMON_COST_SCALING_H