/* -*- C++ -*-

  *

  * This file is a part of LEMON, a generic C++ optimization library

  *

  * Copyright (C) 2003-2008

  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

  * (Egervary Research Group on Combinatorial Optimization, EGRES).

  *

  * Permission to use, modify and distribute this software is granted

  * provided that this copyright notice appears in all copies. For

  * precise terms see the accompanying LICENSE file.

  *

  * This software is provided "AS IS" with no warranty of any kind,

  * express or implied, and with no claim as to its suitability for any

  * purpose.

  *

  */

 #ifndef LEMON_CAPACITY_SCALING_H

 #define LEMON_CAPACITY_SCALING_H

 /// \ingroup min_cost_flow

 ///

 /// \file

 /// \brief The capacity scaling algorithm for finding a minimum cost flow.

 #include <lemon/graph_adaptor.h>

 #include <lemon/bin_heap.h>

 #include <vector>

 namespace lemon {

   /// \addtogroup min_cost_flow

   /// @{

   /// \brief Implementation of the capacity scaling version of the

   /// successive shortest path algorithm for finding a minimum cost

   /// flow.

   ///

   /// \ref CapacityScaling implements the capacity scaling version

   /// of the successive shortest path algorithm for finding a minimum

   /// cost flow.

   ///

   /// \param Graph The directed graph type the algorithm runs on.

   /// \param LowerMap The type of the lower bound map.

   /// \param CapacityMap The type of the capacity (upper bound) map.

   /// \param CostMap The type of the cost (length) map.

   /// \param SupplyMap The type of the supply map.

   ///

   /// \warning

   /// - Edge capacities and costs should be non-negative integers.

   ///   However \c CostMap::Value should be signed type.

   /// - Supply values should be signed integers.

   /// - \c LowerMap::Value must be convertible to

   ///   \c CapacityMap::Value and \c CapacityMap::Value must be

   ///   convertible to \c SupplyMap::Value.

   ///

   /// \author Peter Kovacs

   template < typename Graph,

              typename LowerMap = typename Graph::template EdgeMap<int>,

              typename CapacityMap = LowerMap,

              typename CostMap = typename Graph::template EdgeMap<int>,

              typename SupplyMap = typename Graph::template NodeMap

                                   <typename CapacityMap::Value> >

   class CapacityScaling

   {

     GRAPH_TYPEDEFS(typename Graph);

     typedef typename LowerMap::Value Lower;

     typedef typename CapacityMap::Value Capacity;

     typedef typename CostMap::Value Cost;

     typedef typename SupplyMap::Value Supply;

     typedef typename Graph::template EdgeMap<Capacity> CapacityEdgeMap;

     typedef typename Graph::template NodeMap<Supply> SupplyNodeMap;

     typedef typename Graph::template NodeMap<Edge> PredMap;

   public:

     /// Type to enable or disable capacity scaling.

     enum ScalingEnum {

       WITH_SCALING = 0,

       WITHOUT_SCALING = -1

     };

     /// The type of the flow map.

     typedef typename Graph::template EdgeMap<Capacity> FlowMap;

     /// The type of the potential map.

     typedef typename Graph::template NodeMap<Cost> PotentialMap;

   protected:

     /// \brief Special implementation of the \ref Dijkstra algorithm

     /// for finding shortest paths in the residual network of the graph

     /// with respect to the reduced edge costs and modifying the

     /// node potentials according to the distance of the nodes.

     class ResidualDijkstra

     {

       typedef typename Graph::template NodeMap<Cost> CostNodeMap;

       typedef typename Graph::template NodeMap<Edge> PredMap;

       typedef typename Graph::template NodeMap<int> HeapCrossRef;

       typedef BinHeap<Cost, HeapCrossRef> Heap;

     protected:

       /// The directed graph the algorithm runs on.

       const Graph &graph;

       /// The flow map.

       const FlowMap &flow;

       /// The residual capacity map.

       const CapacityEdgeMap &res_cap;

       /// The cost map.

       const CostMap &cost;

       /// The excess map.

       const SupplyNodeMap &excess;

       /// The potential map.

       PotentialMap &potential;

       /// The distance map.

       CostNodeMap dist;

       /// The map of predecessors edges.

       PredMap &pred;

       /// The processed (i.e. permanently labeled) nodes.

       std::vector<Node> proc_nodes;

     public:

       /// The constructor of the class.

       ResidualDijkstra( const Graph &_graph,

                         const FlowMap &_flow,

                         const CapacityEdgeMap &_res_cap,

                         const CostMap &_cost,

                         const SupplyMap &_excess,

                         PotentialMap &_potential,

                         PredMap &_pred ) :

         graph(_graph), flow(_flow), res_cap(_res_cap), cost(_cost),

         excess(_excess), potential(_potential), dist(_graph),

         pred(_pred)

       {}

       /// Runs the algorithm from the given source node.

       Node run(Node s, Capacity delta) {

         HeapCrossRef heap_cross_ref(graph, Heap::PRE_HEAP);

         Heap heap(heap_cross_ref);

         heap.push(s, 0);

         pred[s] = INVALID;

         proc_nodes.clear();

         // Processing nodes

         while (!heap.empty() && excess[heap.top()] > -delta) {

           Node u = heap.top(), v;

           Cost d = heap.prio() - potential[u], nd;

           dist[u] = heap.prio();

           heap.pop();

           proc_nodes.push_back(u);

           // Traversing outgoing edges

           for (OutEdgeIt e(graph, u); e != INVALID; ++e) {

             if (res_cap[e] >= delta) {

               v = graph.target(e);

               switch(heap.state(v)) {

               case Heap::PRE_HEAP:

                 heap.push(v, d + cost[e] + potential[v]);

                 pred[v] = e;

                 break;

               case Heap::IN_HEAP:

                 nd = d + cost[e] + potential[v];

                 if (nd < heap[v]) {

                   heap.decrease(v, nd);

                   pred[v] = e;

                 }

                 break;

               case Heap::POST_HEAP:

                 break;

               }

             }

           }

           // Traversing incoming edges

           for (InEdgeIt e(graph, u); e != INVALID; ++e) {

             if (flow[e] >= delta) {

               v = graph.source(e);

               switch(heap.state(v)) {

               case Heap::PRE_HEAP:

                 heap.push(v, d - cost[e] + potential[v]);

                 pred[v] = e;

                 break;

               case Heap::IN_HEAP:

                 nd = d - cost[e] + potential[v];

                 if (nd < heap[v]) {

                   heap.decrease(v, nd);

                   pred[v] = e;

                 }

                 break;

               case Heap::POST_HEAP:

                 break;

               }

             }

           }

         }

         if (heap.empty()) return INVALID;

         // Updating potentials of processed nodes

         Node t = heap.top();

         Cost dt = heap.prio();

         for (int i = 0; i < proc_nodes.size(); ++i)

           potential[proc_nodes[i]] -= dist[proc_nodes[i]] - dt;

         return t;

       }

     }; //class ResidualDijkstra

   protected:

     /// The directed graph the algorithm runs on.

     const Graph &graph;

     /// The original lower bound map.

     const LowerMap *lower;

     /// The modified capacity map.

     CapacityEdgeMap capacity;

     /// The cost map.

     const CostMap &cost;

     /// The modified supply map.

     SupplyNodeMap supply;

     bool valid_supply;

     /// The edge map of the current flow.

     FlowMap flow;

     /// The potential node map.

     PotentialMap potential;

     /// The residual capacity map.

     CapacityEdgeMap res_cap;

     /// The excess map.

     SupplyNodeMap excess;

     /// The excess nodes (i.e. the nodes with positive excess).

     std::vector<Node> excess_nodes;

     /// The deficit nodes (i.e. the nodes with negative excess).

     std::vector<Node> deficit_nodes;

     /// The scaling status (enabled or disabled).

     ScalingEnum scaling;

     /// The \c delta parameter used for capacity scaling.

     Capacity delta;

     /// The maximum number of phases.

     int phase_num;

     /// \brief Implementation of the \ref Dijkstra algorithm for

     /// finding augmenting shortest paths in the residual network.

     ResidualDijkstra dijkstra;

     /// The map of predecessors edges.

     PredMap pred;

   public :

     /// \brief General constructor of the class (with lower bounds).

     ///

     /// General constructor of the class (with lower bounds).

     ///

     /// \param _graph The directed graph the algorithm runs on.

     /// \param _lower The lower bounds of the edges.

     /// \param _capacity The capacities (upper bounds) of the edges.

     /// \param _cost The cost (length) values of the edges.

     /// \param _supply The supply values of the nodes (signed).

     CapacityScaling( const Graph &_graph,

                      const LowerMap &_lower,

                      const CapacityMap &_capacity,

                      const CostMap &_cost,

                      const SupplyMap &_supply ) :

       graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),

       supply(_graph), flow(_graph, 0), potential(_graph, 0),

       res_cap(_graph), excess(_graph), pred(_graph),

       dijkstra(graph, flow, res_cap, cost, excess, potential, pred)

     {

       // Removing non-zero lower bounds

       capacity = subMap(_capacity, _lower);

       res_cap = capacity;

       Supply sum = 0;

       for (NodeIt n(graph); n != INVALID; ++n) {

         Supply s = _supply[n];

         for (InEdgeIt e(graph, n); e != INVALID; ++e)

           s += _lower[e];

         for (OutEdgeIt e(graph, n); e != INVALID; ++e)

           s -= _lower[e];

         supply[n] = s;

         sum += s;

       }

       valid_supply = sum == 0;

     }

     /// \brief General constructor of the class (without lower bounds).

     ///

     /// General constructor of the class (without lower bounds).

     ///

     /// \param _graph The directed graph the algorithm runs on.

     /// \param _capacity The capacities (upper bounds) of the edges.

     /// \param _cost The cost (length) values of the edges.

     /// \param _supply The supply values of the nodes (signed).

     CapacityScaling( const Graph &_graph,

                      const CapacityMap &_capacity,

                      const CostMap &_cost,

                      const SupplyMap &_supply ) :

       graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),

       supply(_supply), flow(_graph, 0), potential(_graph, 0),

       res_cap(_capacity), excess(_graph), pred(_graph),

       dijkstra(graph, flow, res_cap, cost, excess, potential)

     {

       // Checking the sum of supply values

       Supply sum = 0;

       for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n];

       valid_supply = sum == 0;

     }

     /// \brief Simple constructor of the class (with lower bounds).

     ///

     /// Simple constructor of the class (with lower bounds).

     ///

     /// \param _graph The directed graph the algorithm runs on.

     /// \param _lower The lower bounds of the edges.

     /// \param _capacity The capacities (upper bounds) of the edges.

     /// \param _cost The cost (length) values of the edges.

     /// \param _s The source node.

     /// \param _t The target node.

     /// \param _flow_value The required amount of flow from node \c _s

     /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t).

     CapacityScaling( const Graph &_graph,

                      const LowerMap &_lower,

                      const CapacityMap &_capacity,

                      const CostMap &_cost,

                      Node _s, Node _t,

                      Supply _flow_value ) :

       graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),

       supply(_graph), flow(_graph, 0), potential(_graph, 0),

       res_cap(_graph), excess(_graph), pred(_graph),

       dijkstra(graph, flow, res_cap, cost, excess, potential)

     {

       // Removing non-zero lower bounds

       capacity = subMap(_capacity, _lower);

       res_cap = capacity;

       for (NodeIt n(graph); n != INVALID; ++n) {

         Supply s = 0;

         if (n == _s) s =  _flow_value;

         if (n == _t) s = -_flow_value;

         for (InEdgeIt e(graph, n); e != INVALID; ++e)

           s += _lower[e];

         for (OutEdgeIt e(graph, n); e != INVALID; ++e)

           s -= _lower[e];

         supply[n] = s;

       }

       valid_supply = true;

     }

     /// \brief Simple constructor of the class (without lower bounds).

     ///

     /// Simple constructor of the class (without lower bounds).

     ///

     /// \param _graph The directed graph the algorithm runs on.

     /// \param _capacity The capacities (upper bounds) of the edges.

     /// \param _cost The cost (length) values of the edges.

     /// \param _s The source node.

     /// \param _t The target node.

     /// \param _flow_value The required amount of flow from node \c _s

     /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t).

     CapacityScaling( const Graph &_graph,

                      const CapacityMap &_capacity,

                      const CostMap &_cost,

                      Node _s, Node _t,

                      Supply _flow_value ) :

       graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),

       supply(_graph, 0), flow(_graph, 0), potential(_graph, 0),

       res_cap(_capacity), excess(_graph), pred(_graph),

       dijkstra(graph, flow, res_cap, cost, excess, potential)

     {

       supply[_s] =  _flow_value;

       supply[_t] = -_flow_value;

       valid_supply = true;

     }

     /// \brief Runs the algorithm.

     ///

     /// Runs the algorithm.

     ///

     /// \param scaling_mode The scaling mode. In case of WITH_SCALING

     /// capacity scaling is enabled in the algorithm (this is the

     /// default) otherwise it is disabled.

     /// If the maximum edge capacity and/or the amount of total supply

     /// is small, the algorithm could be slightly faster without

     /// scaling.

     ///

     /// \return \c true if a feasible flow can be found.

     bool run(int scaling_mode = WITH_SCALING) {

       return init(scaling_mode) && start();

     }

     /// \brief Returns a const reference to the flow map.

     ///

     /// Returns a const reference to the flow map.

     ///

     /// \pre \ref run() must be called before using this function.

     const FlowMap& flowMap() const {

       return flow;

     }

     /// \brief Returns a const reference to the potential map (the dual

     /// solution).

     ///

     /// Returns a const reference to the potential map (the dual

     /// solution).

     ///

     /// \pre \ref run() must be called before using this function.

     const PotentialMap& potentialMap() const {

       return potential;

     }

     /// \brief Returns the total cost of the found flow.

     ///

     /// Returns the total cost of the found flow. The complexity of the

     /// function is \f$ O(e) \f$.

     ///

     /// \pre \ref run() must be called before using this function.

     Cost totalCost() const {

       Cost c = 0;

       for (EdgeIt e(graph); e != INVALID; ++e)

         c += flow[e] * cost[e];

       return c;

     }

   protected:

     /// Initializes the algorithm.

     bool init(int scaling_mode) {

       if (!valid_supply) return false;

       excess = supply;

       // Initilaizing delta value

       if (scaling_mode == WITH_SCALING) {

         // With scaling

         Supply max_sup = 0, max_dem = 0;

         for (NodeIt n(graph); n != INVALID; ++n) {

           if ( supply[n] > max_sup) max_sup =  supply[n];

           if (-supply[n] > max_dem) max_dem = -supply[n];

         }

         if (max_dem < max_sup) max_sup = max_dem;

         phase_num = 0;

         for (delta = 1; 2 * delta <= max_sup; delta *= 2)

           ++phase_num;

       } else {

         // Without scaling

         delta = 1;

       }

       return true;

     }

     /// Executes the algorithm.

     bool start() {

       if (delta > 1)

         return startWithScaling();

       else

         return startWithoutScaling();

     }

     /// \brief Executes the capacity scaling version of the successive

     /// shortest path algorithm.

     bool startWithScaling() {

       // Processing capacity scaling phases

       Node s, t;

       int phase_cnt = 0;

       int factor = 4;

       while (true) {

         // Saturating all edges not satisfying the optimality condition

         for (EdgeIt e(graph); e != INVALID; ++e) {

           Node u = graph.source(e), v = graph.target(e);

           Cost c = cost[e] - potential[u] + potential[v];

           if (c < 0 && res_cap[e] >= delta) {

             excess[u] -= res_cap[e];

             excess[v] += res_cap[e];

             flow[e] = capacity[e];

             res_cap[e] = 0;

           }

           else if (c > 0 && flow[e] >= delta) {

             excess[u] += flow[e];

             excess[v] -= flow[e];

             flow[e] = 0;

             res_cap[e] = capacity[e];

           }

         }

         // Finding excess nodes and deficit nodes

         excess_nodes.clear();

         deficit_nodes.clear();

         for (NodeIt n(graph); n != INVALID; ++n) {

           if (excess[n] >=  delta) excess_nodes.push_back(n);

           if (excess[n] <= -delta) deficit_nodes.push_back(n);

         }

         int next_node = 0;

         // Finding augmenting shortest paths

         while (next_node < excess_nodes.size()) {

           // Checking deficit nodes

           if (delta > 1) {

             bool delta_deficit = false;

             for (int i = 0; i < deficit_nodes.size(); ++i) {

               if (excess[deficit_nodes[i]] <= -delta) {

                 delta_deficit = true;

                 break;

               }

             }

             if (!delta_deficit) break;

           }

           // Running Dijkstra

           s = excess_nodes[next_node];

           if ((t = dijkstra.run(s, delta)) == INVALID) {

             if (delta > 1) {

               ++next_node;

               continue;

             }

             return false;

           }

           // Augmenting along a shortest path from s to t.

           Capacity d = excess[s] < -excess[t] ? excess[s] : -excess[t];

           Node u = t;

           Edge e;

           if (d > delta) {

             while ((e = pred[u]) != INVALID) {

               Capacity rc;

               if (u == graph.target(e)) {

                 rc = res_cap[e];

                 u = graph.source(e);

               } else {

                 rc = flow[e];

                 u = graph.target(e);

               }

               if (rc < d) d = rc;

             }

           }

           u = t;

           while ((e = pred[u]) != INVALID) {

             if (u == graph.target(e)) {

               flow[e] += d;

               res_cap[e] -= d;

               u = graph.source(e);

             } else {

               flow[e] -= d;

               res_cap[e] += d;

               u = graph.target(e);

             }

           }

           excess[s] -= d;

           excess[t] += d;

           if (excess[s] < delta) ++next_node;

         }

         if (delta == 1) break;

         if (++phase_cnt > phase_num / 4) factor = 2;

         delta = delta <= factor ? 1 : delta / factor;

       }

       // Handling non-zero lower bounds

       if (lower) {

         for (EdgeIt e(graph); e != INVALID; ++e)

           flow[e] += (*lower)[e];

       }

       return true;

     }

     /// \brief Executes the successive shortest path algorithm without

     /// capacity scaling.

     bool startWithoutScaling() {

       // Finding excess nodes

       for (NodeIt n(graph); n != INVALID; ++n)

         if (excess[n] > 0) excess_nodes.push_back(n);

       if (excess_nodes.size() == 0) return true;

       int next_node = 0;

       // Finding shortest paths

       Node s, t;

       while ( excess[excess_nodes[next_node]] > 0 ||

               ++next_node < excess_nodes.size() )

       {

         // Running Dijkstra

         s = excess_nodes[next_node];

         if ((t = dijkstra.run(s, 1)) == INVALID)

           return false;

         // Augmenting along a shortest path from s to t

         Capacity d = excess[s] < -excess[t] ? excess[s] : -excess[t];

         Node u = t;

         Edge e;

         while ((e = pred[u]) != INVALID) {

           Capacity rc;

           if (u == graph.target(e)) {

             rc = res_cap[e];

             u = graph.source(e);

           } else {

             rc = flow[e];

             u = graph.target(e);

           }

           if (rc < d) d = rc;

         }

         u = t;

         while ((e = pred[u]) != INVALID) {

           if (u == graph.target(e)) {

             flow[e] += d;

             res_cap[e] -= d;

             u = graph.source(e);

           } else {

             flow[e] -= d;

             res_cap[e] += d;

             u = graph.target(e);

           }

         }

         excess[s] -= d;

         excess[t] += d;

       }

       // Handling non-zero lower bounds

       if (lower) {

         for (EdgeIt e(graph); e != INVALID; ++e)

           flow[e] += (*lower)[e];

       }

       return true;

     }

   }; //class CapacityScaling

   ///@}

 } //namespace lemon

 #endif //LEMON_CAPACITY_SCALING_H