deba@2440: /* -*- C++ -*- deba@2440: * deba@2440: * This file is a part of LEMON, a generic C++ optimization library deba@2440: * alpar@2553: * Copyright (C) 2003-2008 deba@2440: * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport deba@2440: * (Egervary Research Group on Combinatorial Optimization, EGRES). deba@2440: * deba@2440: * Permission to use, modify and distribute this software is granted deba@2440: * provided that this copyright notice appears in all copies. For deba@2440: * precise terms see the accompanying LICENSE file. deba@2440: * deba@2440: * This software is provided "AS IS" with no warranty of any kind, deba@2440: * express or implied, and with no claim as to its suitability for any deba@2440: * purpose. deba@2440: * deba@2440: */ deba@2440: deba@2440: #ifndef LEMON_CAPACITY_SCALING_H deba@2440: #define LEMON_CAPACITY_SCALING_H deba@2440: deba@2440: /// \ingroup min_cost_flow deba@2440: /// deba@2440: /// \file kpeter@2574: /// \brief Capacity scaling algorithm for finding a minimum cost flow. kpeter@2574: kpeter@2574: #include kpeter@2535: #include deba@2457: deba@2440: namespace lemon { deba@2440: deba@2440: /// \addtogroup min_cost_flow deba@2440: /// @{ deba@2440: kpeter@2574: /// \brief Implementation of the capacity scaling algorithm for kpeter@2574: /// finding a minimum cost flow. deba@2440: /// kpeter@2535: /// \ref CapacityScaling implements the capacity scaling version kpeter@2535: /// of the successive shortest path algorithm for finding a minimum kpeter@2535: /// cost flow. deba@2440: /// kpeter@2574: /// \tparam Graph The directed graph type the algorithm runs on. kpeter@2574: /// \tparam LowerMap The type of the lower bound map. kpeter@2574: /// \tparam CapacityMap The type of the capacity (upper bound) map. kpeter@2574: /// \tparam CostMap The type of the cost (length) map. kpeter@2574: /// \tparam SupplyMap The type of the supply map. deba@2440: /// deba@2440: /// \warning kpeter@2574: /// - Edge capacities and costs should be \e non-negative \e integers. kpeter@2574: /// - 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. deba@2440: /// deba@2440: /// \author Peter Kovacs kpeter@2533: template < typename Graph, kpeter@2535: typename LowerMap = typename Graph::template EdgeMap, kpeter@2574: typename CapacityMap = typename Graph::template EdgeMap, kpeter@2535: typename CostMap = typename Graph::template EdgeMap, kpeter@2574: typename SupplyMap = typename Graph::template NodeMap > deba@2440: class CapacityScaling deba@2440: { kpeter@2556: GRAPH_TYPEDEFS(typename Graph); deba@2440: deba@2440: typedef typename CapacityMap::Value Capacity; deba@2440: typedef typename CostMap::Value Cost; deba@2440: typedef typename SupplyMap::Value Supply; kpeter@2556: typedef typename Graph::template EdgeMap CapacityEdgeMap; kpeter@2556: typedef typename Graph::template NodeMap SupplyNodeMap; kpeter@2535: typedef typename Graph::template NodeMap PredMap; deba@2440: deba@2440: public: deba@2440: kpeter@2556: /// The type of the flow map. kpeter@2556: typedef typename Graph::template EdgeMap FlowMap; kpeter@2556: /// The type of the potential map. deba@2440: typedef typename Graph::template NodeMap PotentialMap; deba@2440: kpeter@2574: private: deba@2440: kpeter@2535: /// \brief Special implementation of the \ref Dijkstra algorithm kpeter@2574: /// for finding shortest paths in the residual network. kpeter@2574: /// kpeter@2574: /// \ref ResidualDijkstra is a special implementation of the kpeter@2574: /// \ref Dijkstra algorithm for finding shortest paths in the kpeter@2574: /// residual network of the graph with respect to the reduced edge kpeter@2574: /// costs and modifying the node potentials according to the kpeter@2574: /// distance of the nodes. kpeter@2535: class ResidualDijkstra deba@2440: { kpeter@2535: typedef typename Graph::template NodeMap HeapCrossRef; kpeter@2535: typedef BinHeap Heap; kpeter@2535: kpeter@2574: private: kpeter@2535: kpeter@2574: // The directed graph the algorithm runs on kpeter@2574: const Graph &_graph; kpeter@2535: kpeter@2574: // The main maps kpeter@2574: const FlowMap &_flow; kpeter@2574: const CapacityEdgeMap &_res_cap; kpeter@2574: const CostMap &_cost; kpeter@2574: const SupplyNodeMap &_excess; kpeter@2574: PotentialMap &_potential; kpeter@2535: kpeter@2574: // The distance map kpeter@2588: PotentialMap _dist; kpeter@2574: // The pred edge map kpeter@2574: PredMap &_pred; kpeter@2574: // The processed (i.e. permanently labeled) nodes kpeter@2574: std::vector _proc_nodes; deba@2440: deba@2440: public: deba@2440: kpeter@2581: /// Constructor. kpeter@2574: ResidualDijkstra( const Graph &graph, kpeter@2574: const FlowMap &flow, kpeter@2574: const CapacityEdgeMap &res_cap, kpeter@2574: const CostMap &cost, kpeter@2574: const SupplyMap &excess, kpeter@2574: PotentialMap &potential, kpeter@2574: PredMap &pred ) : kpeter@2574: _graph(graph), _flow(flow), _res_cap(res_cap), _cost(cost), kpeter@2574: _excess(excess), _potential(potential), _dist(graph), kpeter@2574: _pred(pred) kpeter@2535: {} deba@2440: kpeter@2620: /// Run the algorithm from the given source node. kpeter@2588: Node run(Node s, Capacity delta = 1) { kpeter@2574: HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP); kpeter@2535: Heap heap(heap_cross_ref); kpeter@2535: heap.push(s, 0); kpeter@2574: _pred[s] = INVALID; kpeter@2574: _proc_nodes.clear(); kpeter@2535: kpeter@2535: // Processing nodes kpeter@2574: while (!heap.empty() && _excess[heap.top()] > -delta) { kpeter@2535: Node u = heap.top(), v; kpeter@2574: Cost d = heap.prio() + _potential[u], nd; kpeter@2574: _dist[u] = heap.prio(); kpeter@2535: heap.pop(); kpeter@2574: _proc_nodes.push_back(u); kpeter@2535: kpeter@2535: // Traversing outgoing edges kpeter@2574: for (OutEdgeIt e(_graph, u); e != INVALID; ++e) { kpeter@2574: if (_res_cap[e] >= delta) { kpeter@2574: v = _graph.target(e); kpeter@2535: switch(heap.state(v)) { kpeter@2535: case Heap::PRE_HEAP: kpeter@2574: heap.push(v, d + _cost[e] - _potential[v]); kpeter@2574: _pred[v] = e; kpeter@2535: break; kpeter@2535: case Heap::IN_HEAP: kpeter@2574: nd = d + _cost[e] - _potential[v]; kpeter@2535: if (nd < heap[v]) { kpeter@2535: heap.decrease(v, nd); kpeter@2574: _pred[v] = e; kpeter@2535: } kpeter@2535: break; kpeter@2535: case Heap::POST_HEAP: kpeter@2535: break; kpeter@2535: } kpeter@2535: } kpeter@2535: } kpeter@2535: kpeter@2535: // Traversing incoming edges kpeter@2574: for (InEdgeIt e(_graph, u); e != INVALID; ++e) { kpeter@2574: if (_flow[e] >= delta) { kpeter@2574: v = _graph.source(e); kpeter@2535: switch(heap.state(v)) { kpeter@2535: case Heap::PRE_HEAP: kpeter@2574: heap.push(v, d - _cost[e] - _potential[v]); kpeter@2574: _pred[v] = e; kpeter@2535: break; kpeter@2535: case Heap::IN_HEAP: kpeter@2574: nd = d - _cost[e] - _potential[v]; kpeter@2535: if (nd < heap[v]) { kpeter@2535: heap.decrease(v, nd); kpeter@2574: _pred[v] = e; kpeter@2535: } kpeter@2535: break; kpeter@2535: case Heap::POST_HEAP: kpeter@2535: break; kpeter@2535: } kpeter@2535: } kpeter@2535: } kpeter@2535: } kpeter@2535: if (heap.empty()) return INVALID; kpeter@2535: kpeter@2535: // Updating potentials of processed nodes kpeter@2535: Node t = heap.top(); kpeter@2574: Cost t_dist = heap.prio(); kpeter@2574: for (int i = 0; i < int(_proc_nodes.size()); ++i) kpeter@2574: _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist; kpeter@2535: kpeter@2535: return t; deba@2440: } deba@2440: kpeter@2535: }; //class ResidualDijkstra deba@2440: kpeter@2574: private: deba@2440: kpeter@2574: // The directed graph the algorithm runs on kpeter@2574: const Graph &_graph; kpeter@2574: // The original lower bound map kpeter@2574: const LowerMap *_lower; kpeter@2574: // The modified capacity map kpeter@2574: CapacityEdgeMap _capacity; kpeter@2574: // The original cost map kpeter@2574: const CostMap &_cost; kpeter@2574: // The modified supply map kpeter@2574: SupplyNodeMap _supply; kpeter@2574: bool _valid_supply; deba@2440: kpeter@2574: // Edge map of the current flow kpeter@2581: FlowMap *_flow; kpeter@2581: bool _local_flow; kpeter@2574: // Node map of the current potentials kpeter@2581: PotentialMap *_potential; kpeter@2581: bool _local_potential; deba@2440: kpeter@2574: // The residual capacity map kpeter@2574: CapacityEdgeMap _res_cap; kpeter@2574: // The excess map kpeter@2574: SupplyNodeMap _excess; kpeter@2574: // The excess nodes (i.e. nodes with positive excess) kpeter@2574: std::vector _excess_nodes; kpeter@2574: // The deficit nodes (i.e. nodes with negative excess) kpeter@2574: std::vector _deficit_nodes; deba@2440: kpeter@2574: // The delta parameter used for capacity scaling kpeter@2574: Capacity _delta; kpeter@2574: // The maximum number of phases kpeter@2574: int _phase_num; deba@2440: kpeter@2574: // The pred edge map kpeter@2574: PredMap _pred; kpeter@2574: // Implementation of the Dijkstra algorithm for finding augmenting kpeter@2574: // shortest paths in the residual network kpeter@2581: ResidualDijkstra *_dijkstra; deba@2440: kpeter@2581: public: deba@2440: kpeter@2581: /// \brief General constructor (with lower bounds). deba@2440: /// kpeter@2581: /// General constructor (with lower bounds). deba@2440: /// kpeter@2574: /// \param graph The directed graph the algorithm runs on. kpeter@2574: /// \param lower The lower bounds of the edges. kpeter@2574: /// \param capacity The capacities (upper bounds) of the edges. kpeter@2574: /// \param cost The cost (length) values of the edges. kpeter@2574: /// \param supply The supply values of the nodes (signed). kpeter@2574: CapacityScaling( const Graph &graph, kpeter@2574: const LowerMap &lower, kpeter@2574: const CapacityMap &capacity, kpeter@2574: const CostMap &cost, kpeter@2574: const SupplyMap &supply ) : kpeter@2629: _graph(graph), _lower(&lower), _capacity(capacity), _cost(cost), kpeter@2629: _supply(supply), _flow(NULL), _local_flow(false), kpeter@2623: _potential(NULL), _local_potential(false), kpeter@2629: _res_cap(capacity), _excess(supply), _pred(graph), _dijkstra(NULL) deba@2440: { kpeter@2629: // Check the sum of supply values deba@2440: Supply sum = 0; kpeter@2629: for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n]; kpeter@2629: _valid_supply = sum == 0; kpeter@2629: kpeter@2629: // Remove non-zero lower bounds kpeter@2629: typename LowerMap::Value lcap; kpeter@2629: for (EdgeIt e(_graph); e != INVALID; ++e) { kpeter@2629: if ((lcap = lower[e]) != 0) { kpeter@2629: _capacity[e] -= lcap; kpeter@2629: _res_cap[e] -= lcap; kpeter@2629: _supply[_graph.source(e)] -= lcap; kpeter@2629: _supply[_graph.target(e)] += lcap; kpeter@2629: _excess[_graph.source(e)] -= lcap; kpeter@2629: _excess[_graph.target(e)] += lcap; kpeter@2629: } deba@2440: } deba@2440: } deba@2440: kpeter@2581: /// \brief General constructor (without lower bounds). deba@2440: /// kpeter@2581: /// General constructor (without lower bounds). deba@2440: /// kpeter@2574: /// \param graph The directed graph the algorithm runs on. kpeter@2574: /// \param capacity The capacities (upper bounds) of the edges. kpeter@2574: /// \param cost The cost (length) values of the edges. kpeter@2574: /// \param supply The supply values of the nodes (signed). kpeter@2574: CapacityScaling( const Graph &graph, kpeter@2574: const CapacityMap &capacity, kpeter@2574: const CostMap &cost, kpeter@2574: const SupplyMap &supply ) : kpeter@2574: _graph(graph), _lower(NULL), _capacity(capacity), _cost(cost), kpeter@2623: _supply(supply), _flow(NULL), _local_flow(false), kpeter@2623: _potential(NULL), _local_potential(false), kpeter@2629: _res_cap(capacity), _excess(supply), _pred(graph), _dijkstra(NULL) deba@2440: { kpeter@2629: // Check the sum of supply values deba@2440: Supply sum = 0; kpeter@2574: for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n]; kpeter@2574: _valid_supply = sum == 0; deba@2440: } deba@2440: kpeter@2581: /// \brief Simple constructor (with lower bounds). deba@2440: /// kpeter@2581: /// Simple constructor (with lower bounds). deba@2440: /// kpeter@2574: /// \param graph The directed graph the algorithm runs on. kpeter@2574: /// \param lower The lower bounds of the edges. kpeter@2574: /// \param capacity The capacities (upper bounds) of the edges. kpeter@2574: /// \param cost The cost (length) values of the edges. kpeter@2574: /// \param s The source node. kpeter@2574: /// \param t The target node. kpeter@2574: /// \param flow_value The required amount of flow from node \c s kpeter@2574: /// to node \c t (i.e. the supply of \c s and the demand of \c t). kpeter@2574: CapacityScaling( const Graph &graph, kpeter@2574: const LowerMap &lower, kpeter@2574: const CapacityMap &capacity, kpeter@2574: const CostMap &cost, kpeter@2574: Node s, Node t, kpeter@2574: Supply flow_value ) : kpeter@2629: _graph(graph), _lower(&lower), _capacity(capacity), _cost(cost), kpeter@2629: _supply(graph, 0), _flow(NULL), _local_flow(false), kpeter@2623: _potential(NULL), _local_potential(false), kpeter@2629: _res_cap(capacity), _excess(graph, 0), _pred(graph), _dijkstra(NULL) deba@2440: { kpeter@2629: // Remove non-zero lower bounds kpeter@2629: _supply[s] = _excess[s] = flow_value; kpeter@2629: _supply[t] = _excess[t] = -flow_value; kpeter@2629: typename LowerMap::Value lcap; kpeter@2629: for (EdgeIt e(_graph); e != INVALID; ++e) { kpeter@2629: if ((lcap = lower[e]) != 0) { kpeter@2629: _capacity[e] -= lcap; kpeter@2629: _res_cap[e] -= lcap; kpeter@2629: _supply[_graph.source(e)] -= lcap; kpeter@2629: _supply[_graph.target(e)] += lcap; kpeter@2629: _excess[_graph.source(e)] -= lcap; kpeter@2629: _excess[_graph.target(e)] += lcap; kpeter@2629: } deba@2440: } kpeter@2574: _valid_supply = true; deba@2440: } deba@2440: kpeter@2581: /// \brief Simple constructor (without lower bounds). deba@2440: /// kpeter@2581: /// Simple constructor (without lower bounds). deba@2440: /// kpeter@2574: /// \param graph The directed graph the algorithm runs on. kpeter@2574: /// \param capacity The capacities (upper bounds) of the edges. kpeter@2574: /// \param cost The cost (length) values of the edges. kpeter@2574: /// \param s The source node. kpeter@2574: /// \param t The target node. kpeter@2574: /// \param flow_value The required amount of flow from node \c s kpeter@2574: /// to node \c t (i.e. the supply of \c s and the demand of \c t). kpeter@2574: CapacityScaling( const Graph &graph, kpeter@2574: const CapacityMap &capacity, kpeter@2574: const CostMap &cost, kpeter@2574: Node s, Node t, kpeter@2574: Supply flow_value ) : kpeter@2574: _graph(graph), _lower(NULL), _capacity(capacity), _cost(cost), kpeter@2623: _supply(graph, 0), _flow(NULL), _local_flow(false), kpeter@2623: _potential(NULL), _local_potential(false), kpeter@2629: _res_cap(capacity), _excess(graph, 0), _pred(graph), _dijkstra(NULL) deba@2440: { kpeter@2629: _supply[s] = _excess[s] = flow_value; kpeter@2629: _supply[t] = _excess[t] = -flow_value; kpeter@2574: _valid_supply = true; deba@2440: } deba@2440: kpeter@2581: /// Destructor. kpeter@2581: ~CapacityScaling() { kpeter@2581: if (_local_flow) delete _flow; kpeter@2581: if (_local_potential) delete _potential; kpeter@2581: delete _dijkstra; kpeter@2581: } kpeter@2581: kpeter@2620: /// \brief Set the flow map. kpeter@2581: /// kpeter@2620: /// Set the flow map. kpeter@2581: /// kpeter@2581: /// \return \c (*this) kpeter@2581: CapacityScaling& 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@2620: /// \brief Set the potential map. kpeter@2581: /// kpeter@2620: /// Set the potential map. kpeter@2581: /// kpeter@2581: /// \return \c (*this) kpeter@2581: CapacityScaling& 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: kpeter@2581: /// @{ kpeter@2581: kpeter@2620: /// \brief Run the algorithm. kpeter@2556: /// kpeter@2620: /// This function runs the algorithm. kpeter@2556: /// kpeter@2574: /// \param scaling Enable or disable capacity scaling. kpeter@2556: /// If the maximum edge capacity and/or the amount of total supply kpeter@2574: /// is rather small, the algorithm could be slightly faster without kpeter@2556: /// scaling. kpeter@2556: /// kpeter@2556: /// \return \c true if a feasible flow can be found. kpeter@2574: bool run(bool scaling = true) { kpeter@2574: return init(scaling) && start(); kpeter@2556: } kpeter@2556: kpeter@2581: /// @} kpeter@2581: kpeter@2581: /// \name Query Functions kpeter@2620: /// The results of the algorithm can be obtained using these kpeter@2620: /// functions.\n kpeter@2620: /// \ref lemon::CapacityScaling::run() "run()" must be called before kpeter@2620: /// using them. kpeter@2581: kpeter@2581: /// @{ kpeter@2581: kpeter@2620: /// \brief Return a const reference to the edge map storing the kpeter@2574: /// found flow. deba@2440: /// kpeter@2620: /// Return a const reference to the edge map storing the found flow. deba@2440: /// deba@2440: /// \pre \ref run() must be called before using this function. deba@2440: const FlowMap& flowMap() const { kpeter@2581: return *_flow; deba@2440: } deba@2440: kpeter@2620: /// \brief Return a const reference to the node map storing the kpeter@2574: /// found potentials (the dual solution). deba@2440: /// kpeter@2620: /// Return a const reference to the node map storing the found kpeter@2574: /// potentials (the dual solution). deba@2440: /// deba@2440: /// \pre \ref run() must be called before using this function. deba@2440: const PotentialMap& potentialMap() const { kpeter@2581: return *_potential; kpeter@2581: } kpeter@2581: kpeter@2620: /// \brief Return the flow on the given edge. kpeter@2581: /// kpeter@2620: /// Return the flow on the given 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@2620: /// \brief Return the potential of the given node. kpeter@2581: /// kpeter@2620: /// Return the potential of the given 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]; deba@2440: } deba@2440: kpeter@2620: /// \brief Return the total cost of the found flow. deba@2440: /// kpeter@2620: /// Return the total cost of the found flow. The complexity of the deba@2440: /// function is \f$ O(e) \f$. deba@2440: /// deba@2440: /// \pre \ref run() must be called before using this function. deba@2440: Cost totalCost() const { deba@2440: Cost c = 0; kpeter@2574: for (EdgeIt e(_graph); e != INVALID; ++e) kpeter@2581: c += (*_flow)[e] * _cost[e]; deba@2440: return c; deba@2440: } deba@2440: kpeter@2581: /// @} kpeter@2581: kpeter@2574: private: deba@2440: kpeter@2620: /// Initialize the algorithm. kpeter@2574: bool init(bool scaling) { kpeter@2574: if (!_valid_supply) return false; kpeter@2581: kpeter@2581: // Initializing 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: for (EdgeIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0; kpeter@2581: for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0; deba@2440: kpeter@2581: _dijkstra = new ResidualDijkstra( _graph, *_flow, _res_cap, _cost, kpeter@2581: _excess, *_potential, _pred ); kpeter@2581: kpeter@2581: // Initializing delta value kpeter@2574: if (scaling) { kpeter@2535: // With scaling kpeter@2535: Supply max_sup = 0, max_dem = 0; kpeter@2574: for (NodeIt n(_graph); n != INVALID; ++n) { kpeter@2574: if ( _supply[n] > max_sup) max_sup = _supply[n]; kpeter@2574: if (-_supply[n] > max_dem) max_dem = -_supply[n]; kpeter@2535: } kpeter@2588: Capacity max_cap = 0; kpeter@2588: for (EdgeIt e(_graph); e != INVALID; ++e) { kpeter@2588: if (_capacity[e] > max_cap) max_cap = _capacity[e]; kpeter@2588: } kpeter@2588: max_sup = std::min(std::min(max_sup, max_dem), max_cap); kpeter@2574: _phase_num = 0; kpeter@2574: for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) kpeter@2574: ++_phase_num; kpeter@2535: } else { kpeter@2535: // Without scaling kpeter@2574: _delta = 1; deba@2440: } kpeter@2581: deba@2440: return true; deba@2440: } deba@2440: kpeter@2535: bool start() { kpeter@2574: if (_delta > 1) kpeter@2535: return startWithScaling(); kpeter@2535: else kpeter@2535: return startWithoutScaling(); kpeter@2535: } kpeter@2535: kpeter@2620: /// Execute the capacity scaling algorithm. kpeter@2535: bool startWithScaling() { kpeter@2535: // Processing capacity scaling phases kpeter@2535: Node s, t; kpeter@2535: int phase_cnt = 0; kpeter@2535: int factor = 4; kpeter@2535: while (true) { kpeter@2535: // Saturating all edges not satisfying the optimality condition kpeter@2574: for (EdgeIt e(_graph); e != INVALID; ++e) { kpeter@2574: Node u = _graph.source(e), v = _graph.target(e); kpeter@2581: Cost c = _cost[e] + (*_potential)[u] - (*_potential)[v]; kpeter@2574: if (c < 0 && _res_cap[e] >= _delta) { kpeter@2574: _excess[u] -= _res_cap[e]; kpeter@2574: _excess[v] += _res_cap[e]; kpeter@2581: (*_flow)[e] = _capacity[e]; kpeter@2574: _res_cap[e] = 0; kpeter@2535: } kpeter@2581: else if (c > 0 && (*_flow)[e] >= _delta) { kpeter@2581: _excess[u] += (*_flow)[e]; kpeter@2581: _excess[v] -= (*_flow)[e]; kpeter@2581: (*_flow)[e] = 0; kpeter@2574: _res_cap[e] = _capacity[e]; kpeter@2535: } kpeter@2535: } kpeter@2535: kpeter@2535: // Finding excess nodes and deficit nodes kpeter@2574: _excess_nodes.clear(); kpeter@2574: _deficit_nodes.clear(); kpeter@2574: for (NodeIt n(_graph); n != INVALID; ++n) { kpeter@2574: if (_excess[n] >= _delta) _excess_nodes.push_back(n); kpeter@2574: if (_excess[n] <= -_delta) _deficit_nodes.push_back(n); kpeter@2535: } kpeter@2620: int next_node = 0, next_def_node = 0; kpeter@2535: kpeter@2535: // Finding augmenting shortest paths kpeter@2574: while (next_node < int(_excess_nodes.size())) { kpeter@2535: // Checking deficit nodes kpeter@2574: if (_delta > 1) { kpeter@2535: bool delta_deficit = false; kpeter@2620: for ( ; next_def_node < int(_deficit_nodes.size()); kpeter@2620: ++next_def_node ) { kpeter@2620: if (_excess[_deficit_nodes[next_def_node]] <= -_delta) { kpeter@2535: delta_deficit = true; kpeter@2535: break; kpeter@2535: } kpeter@2535: } kpeter@2535: if (!delta_deficit) break; kpeter@2535: } kpeter@2535: kpeter@2535: // Running Dijkstra kpeter@2574: s = _excess_nodes[next_node]; kpeter@2581: if ((t = _dijkstra->run(s, _delta)) == INVALID) { kpeter@2574: if (_delta > 1) { kpeter@2535: ++next_node; kpeter@2535: continue; kpeter@2535: } kpeter@2535: return false; kpeter@2535: } kpeter@2535: kpeter@2535: // Augmenting along a shortest path from s to t. kpeter@2588: Capacity d = std::min(_excess[s], -_excess[t]); kpeter@2535: Node u = t; kpeter@2535: Edge e; kpeter@2574: if (d > _delta) { kpeter@2574: while ((e = _pred[u]) != INVALID) { kpeter@2535: Capacity rc; kpeter@2574: if (u == _graph.target(e)) { kpeter@2574: rc = _res_cap[e]; kpeter@2574: u = _graph.source(e); kpeter@2535: } else { kpeter@2581: rc = (*_flow)[e]; kpeter@2574: u = _graph.target(e); kpeter@2535: } kpeter@2535: if (rc < d) d = rc; kpeter@2535: } kpeter@2535: } kpeter@2535: u = t; kpeter@2574: while ((e = _pred[u]) != INVALID) { kpeter@2574: if (u == _graph.target(e)) { kpeter@2581: (*_flow)[e] += d; kpeter@2574: _res_cap[e] -= d; kpeter@2574: u = _graph.source(e); kpeter@2535: } else { kpeter@2581: (*_flow)[e] -= d; kpeter@2574: _res_cap[e] += d; kpeter@2574: u = _graph.target(e); kpeter@2535: } kpeter@2535: } kpeter@2574: _excess[s] -= d; kpeter@2574: _excess[t] += d; kpeter@2535: kpeter@2574: if (_excess[s] < _delta) ++next_node; kpeter@2535: } kpeter@2535: kpeter@2574: if (_delta == 1) break; kpeter@2574: if (++phase_cnt > _phase_num / 4) factor = 2; kpeter@2574: _delta = _delta <= factor ? 1 : _delta / factor; kpeter@2535: } kpeter@2535: kpeter@2556: // Handling non-zero lower bounds kpeter@2574: if (_lower) { kpeter@2574: for (EdgeIt e(_graph); e != INVALID; ++e) kpeter@2581: (*_flow)[e] += (*_lower)[e]; kpeter@2535: } kpeter@2535: return true; kpeter@2535: } kpeter@2535: kpeter@2620: /// Execute the successive shortest path algorithm. kpeter@2535: bool startWithoutScaling() { deba@2440: // Finding excess nodes kpeter@2574: for (NodeIt n(_graph); n != INVALID; ++n) kpeter@2574: if (_excess[n] > 0) _excess_nodes.push_back(n); kpeter@2574: if (_excess_nodes.size() == 0) return true; kpeter@2556: int next_node = 0; deba@2440: deba@2457: // Finding shortest paths kpeter@2535: Node s, t; kpeter@2574: while ( _excess[_excess_nodes[next_node]] > 0 || kpeter@2574: ++next_node < int(_excess_nodes.size()) ) deba@2440: { kpeter@2535: // Running Dijkstra kpeter@2574: s = _excess_nodes[next_node]; kpeter@2589: if ((t = _dijkstra->run(s)) == INVALID) return false; deba@2440: kpeter@2535: // Augmenting along a shortest path from s to t kpeter@2588: Capacity d = std::min(_excess[s], -_excess[t]); kpeter@2535: Node u = t; kpeter@2535: Edge e; kpeter@2588: if (d > 1) { kpeter@2588: while ((e = _pred[u]) != INVALID) { kpeter@2588: Capacity rc; kpeter@2588: if (u == _graph.target(e)) { kpeter@2588: rc = _res_cap[e]; kpeter@2588: u = _graph.source(e); kpeter@2588: } else { kpeter@2588: rc = (*_flow)[e]; kpeter@2588: u = _graph.target(e); kpeter@2588: } kpeter@2588: if (rc < d) d = rc; kpeter@2535: } kpeter@2535: } kpeter@2535: u = t; kpeter@2574: while ((e = _pred[u]) != INVALID) { kpeter@2574: if (u == _graph.target(e)) { kpeter@2581: (*_flow)[e] += d; kpeter@2574: _res_cap[e] -= d; kpeter@2574: u = _graph.source(e); kpeter@2535: } else { kpeter@2581: (*_flow)[e] -= d; kpeter@2574: _res_cap[e] += d; kpeter@2574: u = _graph.target(e); kpeter@2535: } kpeter@2535: } kpeter@2574: _excess[s] -= d; kpeter@2574: _excess[t] += d; deba@2440: } deba@2440: kpeter@2556: // Handling non-zero lower bounds kpeter@2574: if (_lower) { kpeter@2574: for (EdgeIt e(_graph); e != INVALID; ++e) kpeter@2581: (*_flow)[e] += (*_lower)[e]; deba@2440: } deba@2440: return true; deba@2440: } deba@2440: deba@2440: }; //class CapacityScaling deba@2440: deba@2440: ///@} deba@2440: deba@2440: } //namespace lemon deba@2440: deba@2440: #endif //LEMON_CAPACITY_SCALING_H