/* -*- 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 <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 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 typename Graph::template NodeMap<Edge> PredMap;

  public:

    /// \brief Type to enable or disable capacity scaling.

    enum ScalingEnum {

      WITH_SCALING = 0,

      WITHOUT_SCALING = -1

    };

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

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

      const Graph &graph;

      /// \brief The flow map.

      const FlowMap &flow;

      /// \brief The residual capacity map.

      const CapacityRefMap &res_cap;

      /// \brief The cost map.

      const CostMap &cost;

      /// \brief The excess map.

      const SupplyRefMap &excess;

      /// \brief The potential map.

      PotentialMap &potential;

      /// \brief The distance map.

      CostNodeMap dist;

      /// \brief The map of predecessors edges.

      PredMap &pred;

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

      std::vector<Node> proc_nodes;

    public:

      /// \brief The constructor of the class.

      ResidualDijkstra( const Graph &_graph,

                        const FlowMap &_flow,

                        const CapacityRefMap &_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)

      {}

      /// \brief 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:

    /// \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 capacity map.

    CapacityRefMap res_cap;

    /// \brief The excess map.

    SupplyRefMap excess;

    /// \brief The excess nodes (i.e. nodes with positive excess).

    std::vector<Node> excess_nodes;

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

    int next_node;

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

    ScalingEnum scaling;

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

    Capacity delta;

    /// \brief The maximum number of phases.

    Capacity phase_num;

    /// \brief The deficit nodes.

    std::vector<Node> deficit_nodes;

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

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

    ResidualDijkstra dijkstra;

    /// \brief 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 nonzero 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 nonzero 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 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.

    ///

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

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

    /// default value) otherwise it is disabled.

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

    /// is small, the algorithm could be 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();

    }

  protected:

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

    }

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

        }

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

      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 nonzero 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