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@2556: /// \brief The capacity scaling algorithm for finding a minimum cost flow.
deba@2440: 
kpeter@2535: #include <lemon/graph_adaptor.h>
kpeter@2535: #include <lemon/bin_heap.h>
deba@2440: #include <vector>
deba@2457: 
deba@2440: namespace lemon {
deba@2440: 
deba@2440:   /// \addtogroup min_cost_flow
deba@2440:   /// @{
deba@2440: 
deba@2440:   /// \brief Implementation of the capacity scaling version of the
kpeter@2509:   /// successive shortest path algorithm for finding a minimum cost
kpeter@2509:   /// 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:   ///
deba@2440:   /// \param Graph The directed graph type the algorithm runs on.
deba@2440:   /// \param LowerMap The type of the lower bound map.
deba@2440:   /// \param CapacityMap The type of the capacity (upper bound) map.
deba@2440:   /// \param CostMap The type of the cost (length) map.
deba@2440:   /// \param SupplyMap The type of the supply map.
deba@2440:   ///
deba@2440:   /// \warning
kpeter@2556:   /// - Edge capacities and costs should be non-negative integers.
kpeter@2535:   ///   However \c CostMap::Value should be signed type.
kpeter@2509:   /// - Supply values should be signed integers.
deba@2440:   /// - \c LowerMap::Value must be convertible to
kpeter@2535:   ///   \c CapacityMap::Value and \c CapacityMap::Value must be
kpeter@2535:   ///   convertible to \c SupplyMap::Value.
deba@2440:   ///
deba@2440:   /// \author Peter Kovacs
deba@2440: 
kpeter@2533:   template < typename Graph,
kpeter@2535:              typename LowerMap = typename Graph::template EdgeMap<int>,
kpeter@2535:              typename CapacityMap = LowerMap,
kpeter@2535:              typename CostMap = typename Graph::template EdgeMap<int>,
kpeter@2535:              typename SupplyMap = typename Graph::template NodeMap
kpeter@2535:                                   <typename CapacityMap::Value> >
deba@2440:   class CapacityScaling
deba@2440:   {
kpeter@2556:     GRAPH_TYPEDEFS(typename Graph);
deba@2440: 
deba@2440:     typedef typename LowerMap::Value Lower;
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<Capacity> CapacityEdgeMap;
kpeter@2556:     typedef typename Graph::template NodeMap<Supply> SupplyNodeMap;
kpeter@2535:     typedef typename Graph::template NodeMap<Edge> PredMap;
deba@2440: 
deba@2440:   public:
deba@2440: 
kpeter@2556:     /// Type to enable or disable capacity scaling.
kpeter@2535:     enum ScalingEnum {
kpeter@2535:       WITH_SCALING = 0,
kpeter@2535:       WITHOUT_SCALING = -1
kpeter@2535:     };
kpeter@2535: 
kpeter@2556:     /// The type of the flow map.
kpeter@2556:     typedef typename Graph::template EdgeMap<Capacity> FlowMap;
kpeter@2556:     /// The type of the potential map.
deba@2440:     typedef typename Graph::template NodeMap<Cost> PotentialMap;
deba@2440: 
deba@2440:   protected:
deba@2440: 
kpeter@2535:     /// \brief Special implementation of the \ref Dijkstra algorithm
kpeter@2535:     /// for finding shortest paths in the residual network of the graph
kpeter@2535:     /// with respect to the reduced edge costs and modifying the
kpeter@2535:     /// node potentials according to the distance of the nodes.
kpeter@2535:     class ResidualDijkstra
deba@2440:     {
kpeter@2535:       typedef typename Graph::template NodeMap<Cost> CostNodeMap;
kpeter@2535:       typedef typename Graph::template NodeMap<Edge> PredMap;
deba@2440: 
kpeter@2535:       typedef typename Graph::template NodeMap<int> HeapCrossRef;
kpeter@2535:       typedef BinHeap<Cost, HeapCrossRef> Heap;
kpeter@2535: 
kpeter@2535:     protected:
kpeter@2535: 
kpeter@2556:       /// The directed graph the algorithm runs on.
kpeter@2535:       const Graph &graph;
kpeter@2535: 
kpeter@2556:       /// The flow map.
kpeter@2535:       const FlowMap &flow;
kpeter@2556:       /// The residual capacity map.
kpeter@2556:       const CapacityEdgeMap &res_cap;
kpeter@2556:       /// The cost map.
kpeter@2535:       const CostMap &cost;
kpeter@2556:       /// The excess map.
kpeter@2556:       const SupplyNodeMap &excess;
kpeter@2535: 
kpeter@2556:       /// The potential map.
kpeter@2535:       PotentialMap &potential;
kpeter@2535: 
kpeter@2556:       /// The distance map.
kpeter@2535:       CostNodeMap dist;
kpeter@2556:       /// The map of predecessors edges.
kpeter@2535:       PredMap &pred;
kpeter@2556:       /// The processed (i.e. permanently labeled) nodes.
kpeter@2535:       std::vector<Node> proc_nodes;
deba@2440: 
deba@2440:     public:
deba@2440: 
kpeter@2556:       /// The constructor of the class.
kpeter@2535:       ResidualDijkstra( const Graph &_graph,
kpeter@2535:                         const FlowMap &_flow,
kpeter@2556:                         const CapacityEdgeMap &_res_cap,
kpeter@2535:                         const CostMap &_cost,
kpeter@2535:                         const SupplyMap &_excess,
kpeter@2535:                         PotentialMap &_potential,
kpeter@2535:                         PredMap &_pred ) :
kpeter@2535:         graph(_graph), flow(_flow), res_cap(_res_cap), cost(_cost),
kpeter@2535:         excess(_excess), potential(_potential), dist(_graph),
kpeter@2535:         pred(_pred)
kpeter@2535:       {}
deba@2440: 
kpeter@2556:       /// Runs the algorithm from the given source node.
kpeter@2535:       Node run(Node s, Capacity delta) {
kpeter@2535:         HeapCrossRef heap_cross_ref(graph, Heap::PRE_HEAP);
kpeter@2535:         Heap heap(heap_cross_ref);
kpeter@2535:         heap.push(s, 0);
kpeter@2535:         pred[s] = INVALID;
kpeter@2535:         proc_nodes.clear();
kpeter@2535: 
kpeter@2535:         // Processing nodes
kpeter@2535:         while (!heap.empty() && excess[heap.top()] > -delta) {
kpeter@2535:           Node u = heap.top(), v;
kpeter@2535:           Cost d = heap.prio() - potential[u], nd;
kpeter@2535:           dist[u] = heap.prio();
kpeter@2535:           heap.pop();
kpeter@2535:           proc_nodes.push_back(u);
kpeter@2535: 
kpeter@2535:           // Traversing outgoing edges
kpeter@2535:           for (OutEdgeIt e(graph, u); e != INVALID; ++e) {
kpeter@2535:             if (res_cap[e] >= delta) {
kpeter@2535:               v = graph.target(e);
kpeter@2535:               switch(heap.state(v)) {
kpeter@2535:               case Heap::PRE_HEAP:
kpeter@2535:                 heap.push(v, d + cost[e] + potential[v]);
kpeter@2535:                 pred[v] = e;
kpeter@2535:                 break;
kpeter@2535:               case Heap::IN_HEAP:
kpeter@2535:                 nd = d + cost[e] + potential[v];
kpeter@2535:                 if (nd < heap[v]) {
kpeter@2535:                   heap.decrease(v, nd);
kpeter@2535:                   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@2535:           for (InEdgeIt e(graph, u); e != INVALID; ++e) {
kpeter@2535:             if (flow[e] >= delta) {
kpeter@2535:               v = graph.source(e);
kpeter@2535:               switch(heap.state(v)) {
kpeter@2535:               case Heap::PRE_HEAP:
kpeter@2535:                 heap.push(v, d - cost[e] + potential[v]);
kpeter@2535:                 pred[v] = e;
kpeter@2535:                 break;
kpeter@2535:               case Heap::IN_HEAP:
kpeter@2535:                 nd = d - cost[e] + potential[v];
kpeter@2535:                 if (nd < heap[v]) {
kpeter@2535:                   heap.decrease(v, nd);
kpeter@2535:                   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@2535:         Cost dt = heap.prio();
kpeter@2535:         for (int i = 0; i < proc_nodes.size(); ++i)
kpeter@2535:           potential[proc_nodes[i]] -= dist[proc_nodes[i]] - dt;
kpeter@2535: 
kpeter@2535:         return t;
deba@2440:       }
deba@2440: 
kpeter@2535:     }; //class ResidualDijkstra
deba@2440: 
deba@2440:   protected:
deba@2440: 
kpeter@2556:     /// The directed graph the algorithm runs on.
deba@2440:     const Graph &graph;
kpeter@2556:     /// The original lower bound map.
deba@2440:     const LowerMap *lower;
kpeter@2556:     /// The modified capacity map.
kpeter@2556:     CapacityEdgeMap capacity;
kpeter@2556:     /// The cost map.
deba@2440:     const CostMap &cost;
kpeter@2556:     /// The modified supply map.
kpeter@2556:     SupplyNodeMap supply;
deba@2440:     bool valid_supply;
deba@2440: 
kpeter@2556:     /// The edge map of the current flow.
deba@2440:     FlowMap flow;
kpeter@2556:     /// The potential node map.
deba@2440:     PotentialMap potential;
deba@2440: 
kpeter@2556:     /// The residual capacity map.
kpeter@2556:     CapacityEdgeMap res_cap;
kpeter@2556:     /// The excess map.
kpeter@2556:     SupplyNodeMap excess;
kpeter@2556:     /// The excess nodes (i.e. the nodes with positive excess).
kpeter@2535:     std::vector<Node> excess_nodes;
kpeter@2556:     /// The deficit nodes (i.e. the nodes with negative excess).
kpeter@2556:     std::vector<Node> deficit_nodes;
deba@2440: 
kpeter@2556:     /// The scaling status (enabled or disabled).
kpeter@2535:     ScalingEnum scaling;
kpeter@2556:     /// The \c delta parameter used for capacity scaling.
deba@2440:     Capacity delta;
kpeter@2556:     /// The maximum number of phases.
kpeter@2556:     int phase_num;
deba@2440: 
kpeter@2535:     /// \brief Implementation of the \ref Dijkstra algorithm for
kpeter@2535:     /// finding augmenting shortest paths in the residual network.
kpeter@2535:     ResidualDijkstra dijkstra;
kpeter@2556:     /// The map of predecessors edges.
kpeter@2535:     PredMap pred;
deba@2440: 
deba@2440:   public :
deba@2440: 
deba@2440:     /// \brief General constructor of the class (with lower bounds).
deba@2440:     ///
deba@2440:     /// General constructor of the class (with lower bounds).
deba@2440:     ///
deba@2440:     /// \param _graph The directed graph the algorithm runs on.
deba@2440:     /// \param _lower The lower bounds of the edges.
deba@2440:     /// \param _capacity The capacities (upper bounds) of the edges.
deba@2440:     /// \param _cost The cost (length) values of the edges.
deba@2440:     /// \param _supply The supply values of the nodes (signed).
deba@2440:     CapacityScaling( const Graph &_graph,
kpeter@2535:                      const LowerMap &_lower,
kpeter@2535:                      const CapacityMap &_capacity,
kpeter@2535:                      const CostMap &_cost,
kpeter@2535:                      const SupplyMap &_supply ) :
deba@2440:       graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
deba@2440:       supply(_graph), flow(_graph, 0), potential(_graph, 0),
kpeter@2535:       res_cap(_graph), excess(_graph), pred(_graph),
kpeter@2535:       dijkstra(graph, flow, res_cap, cost, excess, potential, pred)
deba@2440:     {
kpeter@2556:       // Removing non-zero lower bounds
deba@2440:       capacity = subMap(_capacity, _lower);
kpeter@2535:       res_cap = capacity;
deba@2440:       Supply sum = 0;
deba@2440:       for (NodeIt n(graph); n != INVALID; ++n) {
kpeter@2535:         Supply s = _supply[n];
kpeter@2535:         for (InEdgeIt e(graph, n); e != INVALID; ++e)
kpeter@2535:           s += _lower[e];
kpeter@2535:         for (OutEdgeIt e(graph, n); e != INVALID; ++e)
kpeter@2535:           s -= _lower[e];
kpeter@2535:         supply[n] = s;
kpeter@2535:         sum += s;
deba@2440:       }
deba@2440:       valid_supply = sum == 0;
deba@2440:     }
deba@2440: 
deba@2440:     /// \brief General constructor of the class (without lower bounds).
deba@2440:     ///
deba@2440:     /// General constructor of the class (without lower bounds).
deba@2440:     ///
deba@2440:     /// \param _graph The directed graph the algorithm runs on.
deba@2440:     /// \param _capacity The capacities (upper bounds) of the edges.
deba@2440:     /// \param _cost The cost (length) values of the edges.
deba@2440:     /// \param _supply The supply values of the nodes (signed).
deba@2440:     CapacityScaling( const Graph &_graph,
kpeter@2535:                      const CapacityMap &_capacity,
kpeter@2535:                      const CostMap &_cost,
kpeter@2535:                      const SupplyMap &_supply ) :
deba@2440:       graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
deba@2440:       supply(_supply), flow(_graph, 0), potential(_graph, 0),
kpeter@2535:       res_cap(_capacity), excess(_graph), pred(_graph),
kpeter@2535:       dijkstra(graph, flow, res_cap, cost, excess, potential)
deba@2440:     {
deba@2440:       // Checking the sum of supply values
deba@2440:       Supply sum = 0;
deba@2440:       for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n];
deba@2440:       valid_supply = sum == 0;
deba@2440:     }
deba@2440: 
deba@2440:     /// \brief Simple constructor of the class (with lower bounds).
deba@2440:     ///
deba@2440:     /// Simple constructor of the class (with lower bounds).
deba@2440:     ///
deba@2440:     /// \param _graph The directed graph the algorithm runs on.
deba@2440:     /// \param _lower The lower bounds of the edges.
deba@2440:     /// \param _capacity The capacities (upper bounds) of the edges.
deba@2440:     /// \param _cost The cost (length) values of the edges.
deba@2440:     /// \param _s The source node.
deba@2440:     /// \param _t The target node.
deba@2440:     /// \param _flow_value The required amount of flow from node \c _s
kpeter@2556:     /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t).
deba@2440:     CapacityScaling( const Graph &_graph,
kpeter@2535:                      const LowerMap &_lower,
kpeter@2535:                      const CapacityMap &_capacity,
kpeter@2535:                      const CostMap &_cost,
kpeter@2535:                      Node _s, Node _t,
kpeter@2535:                      Supply _flow_value ) :
deba@2440:       graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
deba@2440:       supply(_graph), flow(_graph, 0), potential(_graph, 0),
kpeter@2535:       res_cap(_graph), excess(_graph), pred(_graph),
kpeter@2535:       dijkstra(graph, flow, res_cap, cost, excess, potential)
deba@2440:     {
kpeter@2556:       // Removing non-zero lower bounds
deba@2440:       capacity = subMap(_capacity, _lower);
kpeter@2535:       res_cap = capacity;
deba@2440:       for (NodeIt n(graph); n != INVALID; ++n) {
kpeter@2535:         Supply s = 0;
kpeter@2535:         if (n == _s) s =  _flow_value;
kpeter@2535:         if (n == _t) s = -_flow_value;
kpeter@2535:         for (InEdgeIt e(graph, n); e != INVALID; ++e)
kpeter@2535:           s += _lower[e];
kpeter@2535:         for (OutEdgeIt e(graph, n); e != INVALID; ++e)
kpeter@2535:           s -= _lower[e];
kpeter@2535:         supply[n] = s;
deba@2440:       }
deba@2440:       valid_supply = true;
deba@2440:     }
deba@2440: 
deba@2440:     /// \brief Simple constructor of the class (without lower bounds).
deba@2440:     ///
deba@2440:     /// Simple constructor of the class (without lower bounds).
deba@2440:     ///
deba@2440:     /// \param _graph The directed graph the algorithm runs on.
deba@2440:     /// \param _capacity The capacities (upper bounds) of the edges.
deba@2440:     /// \param _cost The cost (length) values of the edges.
deba@2440:     /// \param _s The source node.
deba@2440:     /// \param _t The target node.
deba@2440:     /// \param _flow_value The required amount of flow from node \c _s
kpeter@2556:     /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t).
deba@2440:     CapacityScaling( const Graph &_graph,
kpeter@2535:                      const CapacityMap &_capacity,
kpeter@2535:                      const CostMap &_cost,
kpeter@2535:                      Node _s, Node _t,
kpeter@2535:                      Supply _flow_value ) :
deba@2440:       graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
deba@2440:       supply(_graph, 0), flow(_graph, 0), potential(_graph, 0),
kpeter@2535:       res_cap(_capacity), excess(_graph), pred(_graph),
kpeter@2535:       dijkstra(graph, flow, res_cap, cost, excess, potential)
deba@2440:     {
deba@2440:       supply[_s] =  _flow_value;
deba@2440:       supply[_t] = -_flow_value;
deba@2440:       valid_supply = true;
deba@2440:     }
deba@2440: 
kpeter@2556:     /// \brief Runs the algorithm.
kpeter@2556:     ///
kpeter@2556:     /// Runs the algorithm.
kpeter@2556:     ///
kpeter@2556:     /// \param scaling_mode The scaling mode. In case of WITH_SCALING
kpeter@2556:     /// capacity scaling is enabled in the algorithm (this is the
kpeter@2556:     /// default) otherwise it is disabled.
kpeter@2556:     /// If the maximum edge capacity and/or the amount of total supply
kpeter@2556:     /// is 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@2556:     bool run(int scaling_mode = WITH_SCALING) {
kpeter@2556:       return init(scaling_mode) && start();
kpeter@2556:     }
kpeter@2556: 
deba@2440:     /// \brief Returns a const reference to the flow map.
deba@2440:     ///
deba@2440:     /// Returns a const reference to the flow map.
deba@2440:     ///
deba@2440:     /// \pre \ref run() must be called before using this function.
deba@2440:     const FlowMap& flowMap() const {
deba@2440:       return flow;
deba@2440:     }
deba@2440: 
deba@2440:     /// \brief Returns a const reference to the potential map (the dual
deba@2440:     /// solution).
deba@2440:     ///
deba@2440:     /// Returns a const reference to the potential map (the dual
deba@2440:     /// solution).
deba@2440:     ///
deba@2440:     /// \pre \ref run() must be called before using this function.
deba@2440:     const PotentialMap& potentialMap() const {
deba@2440:       return potential;
deba@2440:     }
deba@2440: 
deba@2440:     /// \brief Returns the total cost of the found flow.
deba@2440:     ///
deba@2440:     /// Returns 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;
deba@2440:       for (EdgeIt e(graph); e != INVALID; ++e)
kpeter@2535:         c += flow[e] * cost[e];
deba@2440:       return c;
deba@2440:     }
deba@2440: 
deba@2440:   protected:
deba@2440: 
kpeter@2556:     /// Initializes the algorithm.
kpeter@2535:     bool init(int scaling_mode) {
deba@2440:       if (!valid_supply) return false;
kpeter@2535:       excess = supply;
deba@2440: 
deba@2440:       // Initilaizing delta value
kpeter@2535:       if (scaling_mode == WITH_SCALING) {
kpeter@2535:         // With scaling
kpeter@2535:         Supply max_sup = 0, max_dem = 0;
kpeter@2535:         for (NodeIt n(graph); n != INVALID; ++n) {
kpeter@2535:           if ( supply[n] > max_sup) max_sup =  supply[n];
kpeter@2535:           if (-supply[n] > max_dem) max_dem = -supply[n];
kpeter@2535:         }
kpeter@2535:         if (max_dem < max_sup) max_sup = max_dem;
kpeter@2535:         phase_num = 0;
kpeter@2535:         for (delta = 1; 2 * delta <= max_sup; delta *= 2)
kpeter@2535:           ++phase_num;
kpeter@2535:       } else {
kpeter@2535:         // Without scaling
kpeter@2535:         delta = 1;
deba@2440:       }
deba@2440:       return true;
deba@2440:     }
deba@2440: 
kpeter@2556:     /// Executes the algorithm.
kpeter@2535:     bool start() {
kpeter@2535:       if (delta > 1)
kpeter@2535:         return startWithScaling();
kpeter@2535:       else
kpeter@2535:         return startWithoutScaling();
kpeter@2535:     }
kpeter@2535: 
kpeter@2535:     /// \brief Executes the capacity scaling version of the successive
kpeter@2535:     /// shortest path 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@2535:         for (EdgeIt e(graph); e != INVALID; ++e) {
kpeter@2535:           Node u = graph.source(e), v = graph.target(e);
kpeter@2535:           Cost c = cost[e] - potential[u] + potential[v];
kpeter@2535:           if (c < 0 && res_cap[e] >= delta) {
kpeter@2535:             excess[u] -= res_cap[e];
kpeter@2535:             excess[v] += res_cap[e];
kpeter@2535:             flow[e] = capacity[e];
kpeter@2535:             res_cap[e] = 0;
kpeter@2535:           }
kpeter@2535:           else if (c > 0 && flow[e] >= delta) {
kpeter@2535:             excess[u] += flow[e];
kpeter@2535:             excess[v] -= flow[e];
kpeter@2535:             flow[e] = 0;
kpeter@2535:             res_cap[e] = capacity[e];
kpeter@2535:           }
kpeter@2535:         }
kpeter@2535: 
kpeter@2535:         // Finding excess nodes and deficit nodes
kpeter@2535:         excess_nodes.clear();
kpeter@2535:         deficit_nodes.clear();
kpeter@2535:         for (NodeIt n(graph); n != INVALID; ++n) {
kpeter@2535:           if (excess[n] >=  delta) excess_nodes.push_back(n);
kpeter@2535:           if (excess[n] <= -delta) deficit_nodes.push_back(n);
kpeter@2535:         }
kpeter@2556:         int next_node = 0;
kpeter@2535: 
kpeter@2535:         // Finding augmenting shortest paths
kpeter@2535:         while (next_node < excess_nodes.size()) {
kpeter@2535:           // Checking deficit nodes
kpeter@2535:           if (delta > 1) {
kpeter@2535:             bool delta_deficit = false;
kpeter@2535:             for (int i = 0; i < deficit_nodes.size(); ++i) {
kpeter@2535:               if (excess[deficit_nodes[i]] <= -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@2535:           s = excess_nodes[next_node];
kpeter@2535:           if ((t = dijkstra.run(s, delta)) == INVALID) {
kpeter@2535:             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@2535:           Capacity d = excess[s] < -excess[t] ? excess[s] : -excess[t];
kpeter@2535:           Node u = t;
kpeter@2535:           Edge e;
kpeter@2535:           if (d > delta) {
kpeter@2535:             while ((e = pred[u]) != INVALID) {
kpeter@2535:               Capacity rc;
kpeter@2535:               if (u == graph.target(e)) {
kpeter@2535:                 rc = res_cap[e];
kpeter@2535:                 u = graph.source(e);
kpeter@2535:               } else {
kpeter@2535:                 rc = flow[e];
kpeter@2535:                 u = graph.target(e);
kpeter@2535:               }
kpeter@2535:               if (rc < d) d = rc;
kpeter@2535:             }
kpeter@2535:           }
kpeter@2535:           u = t;
kpeter@2535:           while ((e = pred[u]) != INVALID) {
kpeter@2535:             if (u == graph.target(e)) {
kpeter@2535:               flow[e] += d;
kpeter@2535:               res_cap[e] -= d;
kpeter@2535:               u = graph.source(e);
kpeter@2535:             } else {
kpeter@2535:               flow[e] -= d;
kpeter@2535:               res_cap[e] += d;
kpeter@2535:               u = graph.target(e);
kpeter@2535:             }
kpeter@2535:           }
kpeter@2535:           excess[s] -= d;
kpeter@2535:           excess[t] += d;
kpeter@2535: 
kpeter@2535:           if (excess[s] < delta) ++next_node;
kpeter@2535:         }
kpeter@2535: 
kpeter@2535:         if (delta == 1) break;
kpeter@2535:         if (++phase_cnt > phase_num / 4) factor = 2;
kpeter@2535:         delta = delta <= factor ? 1 : delta / factor;
kpeter@2535:       }
kpeter@2535: 
kpeter@2556:       // Handling non-zero lower bounds
kpeter@2535:       if (lower) {
kpeter@2535:         for (EdgeIt e(graph); e != INVALID; ++e)
kpeter@2535:           flow[e] += (*lower)[e];
kpeter@2535:       }
kpeter@2535:       return true;
kpeter@2535:     }
kpeter@2535: 
deba@2440:     /// \brief Executes the successive shortest path algorithm without
deba@2440:     /// capacity scaling.
kpeter@2535:     bool startWithoutScaling() {
deba@2440:       // Finding excess nodes
kpeter@2556:       for (NodeIt n(graph); n != INVALID; ++n)
kpeter@2535:         if (excess[n] > 0) excess_nodes.push_back(n);
deba@2440:       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@2535:       while ( excess[excess_nodes[next_node]] > 0 ||
kpeter@2535:               ++next_node < excess_nodes.size() )
deba@2440:       {
kpeter@2535:         // Running Dijkstra
kpeter@2535:         s = excess_nodes[next_node];
kpeter@2535:         if ((t = dijkstra.run(s, 1)) == INVALID)
kpeter@2535:           return false;
deba@2440: 
kpeter@2535:         // Augmenting along a shortest path from s to t
kpeter@2535:         Capacity d = excess[s] < -excess[t] ? excess[s] : -excess[t];
kpeter@2535:         Node u = t;
kpeter@2535:         Edge e;
kpeter@2535:         while ((e = pred[u]) != INVALID) {
kpeter@2535:           Capacity rc;
kpeter@2535:           if (u == graph.target(e)) {
kpeter@2535:             rc = res_cap[e];
kpeter@2535:             u = graph.source(e);
kpeter@2535:           } else {
kpeter@2535:             rc = flow[e];
kpeter@2535:             u = graph.target(e);
kpeter@2535:           }
kpeter@2535:           if (rc < d) d = rc;
kpeter@2535:         }
kpeter@2535:         u = t;
kpeter@2535:         while ((e = pred[u]) != INVALID) {
kpeter@2535:           if (u == graph.target(e)) {
kpeter@2535:             flow[e] += d;
kpeter@2535:             res_cap[e] -= d;
kpeter@2535:             u = graph.source(e);
kpeter@2535:           } else {
kpeter@2535:             flow[e] -= d;
kpeter@2535:             res_cap[e] += d;
kpeter@2535:             u = graph.target(e);
kpeter@2535:           }
kpeter@2535:         }
kpeter@2535:         excess[s] -= d;
kpeter@2535:         excess[t] += d;
deba@2440:       }
deba@2440: 
kpeter@2556:       // Handling non-zero lower bounds
deba@2440:       if (lower) {
kpeter@2535:         for (EdgeIt e(graph); e != INVALID; ++e)
kpeter@2535:           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