/* -*- C++ -*-

  *

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

  *

  * Copyright (C) 2003-2007

  * 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 <vector>

 #include <lemon/graph_adaptor.h>

 #include <lemon/dijkstra.h>

 #define WITH_SCALING

 #ifdef WITH_SCALING

 #define SCALING_FACTOR	  2

 #endif

 //#define _DEBUG_ITER_

 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 lemon::CapacityScaling "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 nonnegative 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

   {

     typedef typename Graph::Node Node;

     typedef typename Graph::NodeIt NodeIt;

     typedef typename Graph::Edge Edge;

     typedef typename Graph::EdgeIt EdgeIt;

     typedef typename Graph::InEdgeIt InEdgeIt;

     typedef typename Graph::OutEdgeIt OutEdgeIt;

     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> CapacityRefMap;

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

     typedef ResGraphAdaptor< const Graph, Capacity,

 			     CapacityRefMap, CapacityRefMap > ResGraph;

     typedef typename ResGraph::Node ResNode;

     typedef typename ResGraph::NodeIt ResNodeIt;

     typedef typename ResGraph::Edge ResEdge;

     typedef typename ResGraph::EdgeIt ResEdgeIt;

   public:

     /// \brief The type of the flow map.

     typedef CapacityRefMap FlowMap;

     /// \brief The type of the potential map.

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

   protected:

     /// \brief Map adaptor class for handling reduced edge costs.

     class ReducedCostMap : public MapBase<ResEdge, Cost>

     {

     private:

       const ResGraph &gr;

       const CostMap &cost_map;

       const PotentialMap &pot_map;

     public:

       ReducedCostMap( const ResGraph &_gr,

 		      const CostMap &_cost,

 		      const PotentialMap &_pot ) :

 	gr(_gr), cost_map(_cost), pot_map(_pot) {}

       Cost operator[](const ResEdge &e) const {

 	return ResGraph::forward(e) ?

 	   cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)] :

 	  -cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];

       }

     }; //class ReducedCostMap

     /// \brief Map class for the \ref lemon::Dijkstra "Dijkstra"

     /// algorithm to update node potentials.

     class PotentialUpdateMap : public MapBase<ResNode, Cost>

     {

     private:

       PotentialMap *pot;

       typedef std::pair<ResNode, Cost> Pair;

       std::vector<Pair> data;

     public:

       void potentialMap(PotentialMap &_pot) {

 	pot = &_pot;

       }

       void init() {

 	data.clear();

       }

       void set(const ResNode &n, const Cost &v) {

 	data.push_back(Pair(n, v));

       }

       void update() {

 	Cost val = data[data.size()-1].second;

 	for (int i = 0; i < data.size()-1; ++i)

 	  (*pot)[data[i].first] += val - data[i].second;

       }

     }; //class PotentialUpdateMap

 #ifdef WITH_SCALING

     /// \brief Map adaptor class for identifing deficit nodes.

     class DeficitBoolMap : public MapBase<ResNode, bool>

     {

     private:

       const SupplyRefMap &imb;

       const Capacity δ

     public:

       DeficitBoolMap(const SupplyRefMap &_imb, const Capacity &_delta) :

 	imb(_imb), delta(_delta) {}

       bool operator[](const ResNode &n) const {

 	return imb[n] <= -delta;

       }

     }; //class DeficitBoolMap

     /// \brief Map adaptor class for filtering edges with at least

     /// \c delta residual capacity.

     class DeltaFilterMap : public MapBase<ResEdge, bool>

     {

     private:

       const ResGraph &gr;

       const Capacity δ

     public:

       DeltaFilterMap(const ResGraph &_gr, const Capacity &_delta) :

 	gr(_gr), delta(_delta) {}

       bool operator[](const ResEdge &e) const {

 	return gr.rescap(e) >= delta;

       }

     }; //class DeltaFilterMap

     typedef EdgeSubGraphAdaptor<const ResGraph, const DeltaFilterMap>

       DeltaResGraph;

     /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class.

     class ResDijkstraTraits :

       public DijkstraDefaultTraits<DeltaResGraph, ReducedCostMap>

     {

     public:

       typedef PotentialUpdateMap DistMap;

       static DistMap *createDistMap(const DeltaResGraph&) {

 	return new DistMap();

       }

     }; //class ResDijkstraTraits

 #else //WITHOUT_CAPACITY_SCALING

     /// \brief Map adaptor class for identifing deficit nodes.

     class DeficitBoolMap : public MapBase<ResNode, bool>

     {

     private:

       const SupplyRefMap &imb;

     public:

       DeficitBoolMap(const SupplyRefMap &_imb) : imb(_imb) {}

       bool operator[](const ResNode &n) const {

 	return imb[n] < 0;

       }

     }; //class DeficitBoolMap

     /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class.

     class ResDijkstraTraits :

       public DijkstraDefaultTraits<ResGraph, ReducedCostMap>

     {

     public:

       typedef PotentialUpdateMap DistMap;

       static DistMap *createDistMap(const ResGraph&) {

 	return new DistMap();

       }

     }; //class ResDijkstraTraits

 #endif

   protected:

     /// \brief The directed graph the algorithm runs on.

     const Graph &graph;

     /// \brief The original lower bound map.

     const LowerMap *lower;

     /// \brief The modified capacity map.

     CapacityRefMap capacity;

     /// \brief The cost map.

     const CostMap &cost;

     /// \brief The modified supply map.

     SupplyRefMap supply;

     /// \brief The sum of supply values equals zero.

     bool valid_supply;

     /// \brief The edge map of the current flow.

     FlowMap flow;

     /// \brief The potential node map.

     PotentialMap potential;

     /// \brief The residual graph.

     ResGraph res_graph;

     /// \brief The reduced cost map.

     ReducedCostMap red_cost;

     /// \brief The imbalance map.

     SupplyRefMap imbalance;

     /// \brief The excess nodes.

     std::vector<ResNode> excess_nodes;

     /// \brief The index of the next excess node.

     int next_node;

 #ifdef WITH_SCALING

     typedef Dijkstra<DeltaResGraph, ReducedCostMap, ResDijkstraTraits>

       ResDijkstra;

     /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding

     /// shortest paths in the residual graph with respect to the

     /// reduced edge costs.

     ResDijkstra dijkstra;

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

     Capacity delta;

     /// \brief Edge filter map for filtering edges with at least

     /// \c delta residual capacity.

     DeltaFilterMap delta_filter;

     /// \brief The delta residual graph (i.e. the subgraph of the

     /// residual graph consisting of edges with at least \c delta

     /// residual capacity).

     DeltaResGraph dres_graph;

     /// \brief Map for identifing deficit nodes.

     DeficitBoolMap delta_deficit;

     /// \brief The deficit nodes.

     std::vector<ResNode> deficit_nodes;

 #else //WITHOUT_CAPACITY_SCALING

     typedef Dijkstra<ResGraph, ReducedCostMap, ResDijkstraTraits>

       ResDijkstra;

     /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding

     /// shortest paths in the residual graph with respect to the

     /// reduced edge costs.

     ResDijkstra dijkstra;

     /// \brief Map for identifing deficit nodes.

     DeficitBoolMap has_deficit;

 #endif

     /// \brief Pred map for the \ref lemon::Dijkstra "Dijkstra" class.

     typename ResDijkstra::PredMap pred;

     /// \brief Dist map for the \ref lemon::Dijkstra "Dijkstra" class to

     /// update node potentials.

     PotentialUpdateMap updater;

   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_graph(_graph, capacity, flow),

       red_cost(res_graph, cost, potential), imbalance(_graph),

 #ifdef WITH_SCALING

       delta(0), delta_filter(res_graph, delta),

       dres_graph(res_graph, delta_filter),

       dijkstra(dres_graph, red_cost), pred(dres_graph),

       delta_deficit(imbalance, delta)

 #else //WITHOUT_CAPACITY_SCALING

       dijkstra(res_graph, red_cost), pred(res_graph),

       has_deficit(imbalance)

 #endif

     {

       // Removing nonzero lower bounds

       capacity = subMap(_capacity, _lower);

       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_graph(_graph, capacity, flow),

       red_cost(res_graph, cost, potential), imbalance(_graph),

 #ifdef WITH_SCALING

       delta(0), delta_filter(res_graph, delta),

       dres_graph(res_graph, delta_filter),

       dijkstra(dres_graph, red_cost), pred(dres_graph),

       delta_deficit(imbalance, delta)

 #else //WITHOUT_CAPACITY_SCALING

       dijkstra(res_graph, red_cost), pred(res_graph),

       has_deficit(imbalance)

 #endif

     {

       // 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_graph(_graph, capacity, flow),

       red_cost(res_graph, cost, potential), imbalance(_graph),

 #ifdef WITH_SCALING

       delta(0), delta_filter(res_graph, delta),

       dres_graph(res_graph, delta_filter),

       dijkstra(dres_graph, red_cost), pred(dres_graph),

       delta_deficit(imbalance, delta)

 #else //WITHOUT_CAPACITY_SCALING

       dijkstra(res_graph, red_cost), pred(res_graph),

       has_deficit(imbalance)

 #endif

     {

       // Removing nonzero lower bounds

       capacity = subMap(_capacity, _lower);

       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_graph(_graph, capacity, flow),

       red_cost(res_graph, cost, potential), imbalance(_graph),

 #ifdef WITH_SCALING

       delta(0), delta_filter(res_graph, delta),

       dres_graph(res_graph, delta_filter),

       dijkstra(dres_graph, red_cost), pred(dres_graph),

       delta_deficit(imbalance, delta)

 #else //WITHOUT_CAPACITY_SCALING

       dijkstra(res_graph, red_cost), pred(res_graph),

       has_deficit(imbalance)

 #endif

     {

       supply[_s] =  _flow_value;

       supply[_t] = -_flow_value;

       valid_supply = true;

     }

     /// \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;

     }

     /// \brief Runs the algorithm.

     ///

     /// Runs the algorithm.

     ///

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

     bool run() {

       return init() && start();

     }

   protected:

     /// \brief Initializes the algorithm.

     bool init() {

       if (!valid_supply) return false;

       imbalance = supply;

       // Initalizing Dijkstra class

       updater.potentialMap(potential);

       dijkstra.distMap(updater).predMap(pred);

 #ifdef WITH_SCALING

       // Initilaizing delta value

       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;

       for ( delta = 1; SCALING_FACTOR * delta < max_sup;

 	    delta *= SCALING_FACTOR ) ;

 #endif

       return true;

     }

 #ifdef WITH_SCALING

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

     /// shortest path algorithm.

     bool start() {

       typedef typename DeltaResGraph::EdgeIt DeltaResEdgeIt;

       typedef typename DeltaResGraph::Edge DeltaResEdge;

 #ifdef _DEBUG_ITER_

       int dijk_num = 0;

 #endif

       // Processing capacity scaling phases

       ResNode s, t;

       for ( ; delta >= 1; delta = delta < SCALING_FACTOR && delta > 1 ?

 				  1 : delta / SCALING_FACTOR )

       {

 	// Saturating edges not satisfying the optimality condition

 	Capacity r;

 	for (DeltaResEdgeIt e(dres_graph); e != INVALID; ++e) {

 	  if (red_cost[e] < 0) {

 	    res_graph.augment(e, r = res_graph.rescap(e));

 	    imbalance[dres_graph.source(e)] -= r;

 	    imbalance[dres_graph.target(e)] += r;

 	  }

 	}

 	// Finding excess nodes

 	excess_nodes.clear();

 	deficit_nodes.clear();

 	for (ResNodeIt n(res_graph); n != INVALID; ++n) {

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

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

 	}

 	next_node = 0;

 	// Finding shortest paths

 	while (next_node < excess_nodes.size()) {

 	  // Checking deficit nodes

 	  if (delta > 1) {

 	    bool delta_def = false;

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

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

 		delta_def = true;

 		break;

 	      }

 	    }

 	    if (!delta_def) break;

 	  }

 	  // Running Dijkstra

 	  s = excess_nodes[next_node];

 	  updater.init();

 	  dijkstra.init();

 	  dijkstra.addSource(s);

 #ifdef _DEBUG_ITER_

 	  ++dijk_num;

 #endif

 	  if ((t = dijkstra.start(delta_deficit)) == INVALID) {

 	    if (delta > 1) {

 	      ++next_node;

 	      continue;

 	    }

 	    return false;

 	  }

 	  // Updating node potentials

 	  updater.update();

 	  // Augment along a shortest path from s to t

 	  Capacity d = imbalance[s] < -imbalance[t] ?

 	    imbalance[s] : -imbalance[t];

 	  ResNode u = t;

 	  ResEdge e;

 	  if (d > delta) {

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

 	      if (res_graph.rescap(e) < d) d = res_graph.rescap(e);

 	      u = dres_graph.source(e);

 	    }

 	  }

 	  u = t;

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

 	    res_graph.augment(e, d);

 	    u = dres_graph.source(e);

 	  }

 	  imbalance[s] -= d;

 	  imbalance[t] += d;

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

 	}

       }

 #ifdef _DEBUG_ITER_

       std::cout << "Capacity Scaling algorithm finished with running Dijkstra algorithm "

 	<< dijk_num << " times." << std::endl;

 #endif

       // Handling nonzero lower bounds

       if (lower) {

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

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

       }

       return true;

     }

 #else //WITHOUT_CAPACITY_SCALING

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

     /// capacity scaling.

     bool start() {

       // Finding excess nodes

       for (ResNodeIt n(res_graph); n != INVALID; ++n) {

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

       }

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

       next_node = 0;

       // Finding shortest paths

       ResNode s, t;

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

 	      ++next_node < excess_nodes.size() )

       {

 	// Running Dijkstra

 	s = excess_nodes[next_node];

 	updater.init();

 	dijkstra.init();

 	dijkstra.addSource(s);

 	if ((t = dijkstra.start(has_deficit)) == INVALID)

 	  return false;

 	// Updating node potentials

 	updater.update();

 	// Augmenting along a shortest path from s to t

 	Capacity delta = imbalance[s] < -imbalance[t] ?

 	  imbalance[s] : -imbalance[t];

 	ResNode u = t;

 	ResEdge e;

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

 	  if (res_graph.rescap(e) < delta) delta = res_graph.rescap(e);

 	  u = res_graph.source(e);

 	}

 	u = t;

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

 	  res_graph.augment(e, delta);

 	  u = res_graph.source(e);

 	}

 	imbalance[s] -= delta;

 	imbalance[t] += delta;

       }

       // Handling nonzero lower bounds

       if (lower) {

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

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

       }

       return true;

     }

 #endif

   }; //class CapacityScaling

   ///@}

 } //namespace lemon

 #endif //LEMON_CAPACITY_SCALING_H