/* -*- C++ -*-

 *

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

 *

 * Copyright (C) 2003-2006

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

#define LEMON_HAO_ORLIN_H

#include <vector>

#include <queue>

#include <limits>

#include <lemon/maps.h>

#include <lemon/graph_utils.h>

#include <lemon/graph_adaptor.h>

#include <lemon/iterable_maps.h>

/// \file

/// \ingroup flowalgs

/// \brief Implementation of the Hao-Orlin algorithm.

///

/// Implementation of the HaoOrlin algorithms class for testing network 

/// reliability.

namespace lemon {

  /// \ingroup flowalgs

  ///

  /// \brief %Hao-Orlin algorithm for calculate minimum cut in directed graphs.

  ///

  /// Hao-Orlin calculates the minimum cut in directed graphs. It

  /// separates the nodes of the graph into two disjoint sets, 

  /// \f$ V_{out} \f$ and \f$ V_{in} \f$. This separation is the minimum

  /// cut if the summary of the edge capacities which source is in

  /// \f$ V_{out} \f$ and the target is in \f$ V_{in} \f$ is the

  /// minimum.  The algorithm is a modified push-relabel preflow

  /// algorithm and it calculates the minimum cut in \f$ O(n^3) \f$

  /// time. The purpose of such algorithm is testing network

  /// reliability. For sparse undirected graph you can use the

  /// algorithm of Nagamochi and Ibaraki which solves the undirected

  /// problem in \f$ O(ne + n^2 \log(n)) \f$ time and it is implemented in the

  /// MinCut algorithm class.

  ///

  /// \param _Graph is the graph type of the algorithm.

  /// \param _CapacityMap is an edge map of capacities which should

  /// be any numreric type. The default type is _Graph::EdgeMap<int>.

  /// \param _Tolerance is the handler of the inexact computation. The

  /// default type for it is Tolerance<typename CapacityMap::Value>.

  ///

  /// \author Attila Bernath and Balazs Dezso

#ifdef DOXYGEN

  template <typename _Graph, typename _CapacityMap, typename _Tolerance>

#else

  template <typename _Graph,

	    typename _CapacityMap = typename _Graph::template EdgeMap<int>,

            typename _Tolerance = Tolerance<typename _CapacityMap::Value> >

#endif

  class HaoOrlin {

  protected:

    typedef _Graph Graph;

    typedef _CapacityMap CapacityMap;

    typedef _Tolerance Tolerance;

    typedef typename CapacityMap::Value Value;

    typedef typename Graph::Node Node;

    typedef typename Graph::NodeIt NodeIt;

    typedef typename Graph::EdgeIt EdgeIt;

    typedef typename Graph::OutEdgeIt OutEdgeIt;

    typedef typename Graph::InEdgeIt InEdgeIt;

    const Graph* _graph;

    const CapacityMap* _capacity;

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

    FlowMap* _preflow;

    Node _source, _target;

    int _node_num;

    typedef ResGraphAdaptor<const Graph, Value, CapacityMap, 

                            FlowMap, Tolerance> OutResGraph;

    typedef typename OutResGraph::Edge OutResEdge;

    OutResGraph* _out_res_graph;

    typedef typename Graph::template NodeMap<OutResEdge> OutCurrentEdgeMap;

    OutCurrentEdgeMap* _out_current_edge;  

    typedef RevGraphAdaptor<const Graph> RevGraph;

    RevGraph* _rev_graph;

    typedef ResGraphAdaptor<const RevGraph, Value, CapacityMap, 

                            FlowMap, Tolerance> InResGraph;

    typedef typename InResGraph::Edge InResEdge;

    InResGraph* _in_res_graph;

    typedef typename Graph::template NodeMap<InResEdge> InCurrentEdgeMap;

    InCurrentEdgeMap* _in_current_edge;  

    typedef IterableBoolMap<Graph, Node> WakeMap;

    WakeMap* _wake;

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

    DistMap* _dist;  

    typedef typename Graph::template NodeMap<Value> ExcessMap;

    ExcessMap* _excess;

    typedef typename Graph::template NodeMap<bool> SourceSetMap;

    SourceSetMap* _source_set;

    std::vector<int> _level_size;

    int _highest_active;

    std::vector<std::list<Node> > _active_nodes;

    int _dormant_max;

    std::vector<std::list<Node> > _dormant;

    Value _min_cut;

    typedef typename Graph::template NodeMap<bool> MinCutMap;

    MinCutMap* _min_cut_map;

    Tolerance _tolerance;

  public: 

    /// \brief Constructor

    ///

    /// Constructor of the algorithm class. 

    HaoOrlin(const Graph& graph, const CapacityMap& capacity, 

             const Tolerance& tolerance = Tolerance()) :

      _graph(&graph), _capacity(&capacity), 

      _preflow(0), _source(), _target(), 

      _out_res_graph(0), _out_current_edge(0),

      _rev_graph(0), _in_res_graph(0), _in_current_edge(0),

      _wake(0),_dist(0), _excess(0), _source_set(0), 

      _highest_active(), _active_nodes(), _dormant_max(), _dormant(), 

      _min_cut(), _min_cut_map(0), _tolerance(tolerance) {}

    ~HaoOrlin() {

      if (_min_cut_map) {

        delete _min_cut_map;

      } 

      if (_in_current_edge) {

        delete _in_current_edge;

      }

      if (_in_res_graph) {

        delete _in_res_graph;

      }

      if (_rev_graph) {

        delete _rev_graph;

      }

      if (_out_current_edge) {

        delete _out_current_edge;

      }

      if (_out_res_graph) {

        delete _out_res_graph;

      }

      if (_source_set) {

        delete _source_set;

      }

      if (_excess) {

        delete _excess;

      }

      if (_dist) {

        delete _dist;

      }

      if (_wake) {

        delete _wake;

      }

      if (_preflow) {

        delete _preflow;

      }

    }

  private:

    template <typename ResGraph, typename EdgeMap>

    void findMinCut(const Node& target, bool out, 

                    ResGraph& res_graph, EdgeMap& current_edge) {

      typedef typename ResGraph::Edge ResEdge;

      typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;

      for (typename Graph::EdgeIt it(*_graph); it != INVALID; ++it) {

        (*_preflow)[it] = 0;      

      }

      for (NodeIt it(*_graph); it != INVALID; ++it) {

        (*_wake)[it] = true;

        (*_dist)[it] = 1;

        (*_excess)[it] = 0;

        (*_source_set)[it] = false;

        res_graph.firstOut(current_edge[it], it);

      }

      _target = target;

      (*_dist)[target] = 0;

      for (ResOutEdgeIt it(res_graph, _source); it != INVALID; ++it) {

        Value delta = res_graph.rescap(it);

        if (!_tolerance.positive(delta)) continue;

        (*_excess)[res_graph.source(it)] -= delta;

        res_graph.augment(it, delta);

        Node a = res_graph.target(it);

        if (!_tolerance.positive((*_excess)[a]) && 

            (*_wake)[a] && a != _target) {

          _active_nodes[(*_dist)[a]].push_front(a);

          if (_highest_active < (*_dist)[a]) {

            _highest_active = (*_dist)[a];

          }

        }

        (*_excess)[a] += delta;

      }

      _dormant[0].push_front(_source);

      (*_source_set)[_source] = true;

      _dormant_max = 0;

      (*_wake)[_source] = false;

      _level_size[0] = 1;

      _level_size[1] = _node_num - 1;

      do {

	Node n;

	while ((n = findActiveNode()) != INVALID) {

	  ResEdge e;

	  while (_tolerance.positive((*_excess)[n]) && 

                 (e = findAdmissibleEdge(n, res_graph, current_edge)) 

                 != INVALID){

	    Value delta;

	    if ((*_excess)[n] < res_graph.rescap(e)) {

	      delta = (*_excess)[n];

	    } else {

	      delta = res_graph.rescap(e);

	      res_graph.nextOut(current_edge[n]);

	    }

            if (!_tolerance.positive(delta)) continue;

	    res_graph.augment(e, delta);

	    (*_excess)[res_graph.source(e)] -= delta;

	    Node a = res_graph.target(e);

	    if (!_tolerance.positive((*_excess)[a]) && a != _target) {

	      _active_nodes[(*_dist)[a]].push_front(a);

	    }

	    (*_excess)[a] += delta;

	  }

	  if (_tolerance.positive((*_excess)[n])) {

	    relabel(n, res_graph, current_edge);

          }

	}

	Value current_value = cutValue(out);

 	if (_min_cut > current_value){

          if (out) {

            for (NodeIt it(*_graph); it != INVALID; ++it) {

              _min_cut_map->set(it, !(*_wake)[it]);

            } 

          } else {

            for (NodeIt it(*_graph); it != INVALID; ++it) {

              _min_cut_map->set(it, (*_wake)[it]);

            } 

          }

	  _min_cut = current_value;

 	}

      } while (selectNewSink(res_graph));

    }

    template <typename ResGraph, typename EdgeMap>

    void relabel(const Node& n, ResGraph& res_graph, EdgeMap& current_edge) {

      typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;

      int k = (*_dist)[n];

      if (_level_size[k] == 1) {

	++_dormant_max;

	for (NodeIt it(*_graph); it != INVALID; ++it) {

	  if ((*_wake)[it] && (*_dist)[it] >= k) {

	    (*_wake)[it] = false;

	    _dormant[_dormant_max].push_front(it);

	    --_level_size[(*_dist)[it]];

	  }

	}

	--_highest_active;

      } else {	

        int new_dist = _node_num;

        for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {

          Node t = res_graph.target(e);

          if ((*_wake)[t] && new_dist > (*_dist)[t]) {

            new_dist = (*_dist)[t];

          }

        }

        if (new_dist == _node_num) {

	  ++_dormant_max;

	  (*_wake)[n] = false;

	  _dormant[_dormant_max].push_front(n);

	  --_level_size[(*_dist)[n]];

	} else {	    

	  --_level_size[(*_dist)[n]];

	  (*_dist)[n] = new_dist + 1;

	  _highest_active = (*_dist)[n];

	  _active_nodes[_highest_active].push_front(n);

	  ++_level_size[(*_dist)[n]];

	  res_graph.firstOut(current_edge[n], n);

	}

      }

    }

    template <typename ResGraph>

    bool selectNewSink(ResGraph& res_graph) {

      typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;

      Node old_target = _target;

      (*_wake)[_target] = false;

      --_level_size[(*_dist)[_target]];

      _dormant[0].push_front(_target);

      (*_source_set)[_target] = true;

      if ((int)_dormant[0].size() == _node_num){

        _dormant[0].clear();

	return false;

      }

      bool wake_was_empty = false;

      if(_wake->trueNum() == 0) {

	while (!_dormant[_dormant_max].empty()){

	  (*_wake)[_dormant[_dormant_max].front()] = true;

	  ++_level_size[(*_dist)[_dormant[_dormant_max].front()]];

	  _dormant[_dormant_max].pop_front();

	}

	--_dormant_max;

	wake_was_empty = true;

      }

      int min_dist = std::numeric_limits<int>::max();

      for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {

	if (min_dist > (*_dist)[it]){

	  _target = it;

	  min_dist = (*_dist)[it];

	}

      }

      if (wake_was_empty){

	for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {

          if (_tolerance.positive((*_excess)[it])) {

	    if ((*_wake)[it] && it != _target) {

	      _active_nodes[(*_dist)[it]].push_front(it);

            }

	    if (_highest_active < (*_dist)[it]) {

	      _highest_active = (*_dist)[it];		    

            }

	  }

	}

      }

      for (ResOutEdgeIt e(res_graph, old_target); e!=INVALID; ++e){

	if (!(*_source_set)[res_graph.target(e)]) {

          Value delta = res_graph.rescap(e);

          if (!_tolerance.positive(delta)) continue;

          res_graph.augment(e, delta);

          (*_excess)[res_graph.source(e)] -= delta;

          Node a = res_graph.target(e);

          if (!_tolerance.positive((*_excess)[a]) && 

              (*_wake)[a] && a != _target) {

            _active_nodes[(*_dist)[a]].push_front(a);

            if (_highest_active < (*_dist)[a]) {

              _highest_active = (*_dist)[a];

            }

          }

          (*_excess)[a] += delta;

	}

      }

      return true;

    }

    Node findActiveNode() {

      while (_highest_active > 0 && _active_nodes[_highest_active].empty()){ 

	--_highest_active;

      }

      if( _highest_active > 0) {

       	Node n = _active_nodes[_highest_active].front();

	_active_nodes[_highest_active].pop_front();

	return n;

      } else {

	return INVALID;

      }

    }

    template <typename ResGraph, typename EdgeMap>

    typename ResGraph::Edge findAdmissibleEdge(const Node& n, 

                                               ResGraph& res_graph, 

                                               EdgeMap& current_edge) {

      typedef typename ResGraph::Edge ResEdge;

      ResEdge e = current_edge[n];

      while (e != INVALID && 

             ((*_dist)[n] <= (*_dist)[res_graph.target(e)] || 

              !(*_wake)[res_graph.target(e)])) {

	res_graph.nextOut(e);

      }

      if (e != INVALID) {

	current_edge[n] = e;	

	return e;

      } else {

	return INVALID;

      }

    }

    Value cutValue(bool out) {

      Value value = 0;

      if (out) {

        for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {

          for (InEdgeIt e(*_graph, it); e != INVALID; ++e) {

            if (!(*_wake)[_graph->source(e)]){

              value += (*_capacity)[e];

            }	

          }

        }

      } else {

        for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {

          for (OutEdgeIt e(*_graph, it); e != INVALID; ++e) {

            if (!(*_wake)[_graph->target(e)]){

              value += (*_capacity)[e];

            }	

          }

        }

      }

      return value;

    }

  public:

    /// \name Execution control

    /// The simplest way to execute the algorithm is to use

    /// one of the member functions called \c run(...).

    /// \n

    /// If you need more control on the execution,

    /// first you must call \ref init(), then the \ref calculateIn() or

    /// \ref calculateIn() functions.

    /// @{

    /// \brief Initializes the internal data structures.

    ///

    /// Initializes the internal data structures. It creates

    /// the maps, residual graph adaptors and some bucket structures

    /// for the algorithm. 

    void init() {

      init(NodeIt(*_graph));

    }

    /// \brief Initializes the internal data structures.

    ///

    /// Initializes the internal data structures. It creates

    /// the maps, residual graph adaptor and some bucket structures

    /// for the algorithm. The \c source node is used as the push-relabel

    /// algorithm's source.

    void init(const Node& source) {

      _source = source;

      _node_num = countNodes(*_graph);

      _dormant.resize(_node_num);

      _level_size.resize(_node_num, 0);

      _active_nodes.resize(_node_num);

      if (!_preflow) {

        _preflow = new FlowMap(*_graph);

      }

      if (!_wake) {

        _wake = new WakeMap(*_graph);

      }

      if (!_dist) {

        _dist = new DistMap(*_graph);

      }

      if (!_excess) {

        _excess = new ExcessMap(*_graph);

      }

      if (!_source_set) {

        _source_set = new SourceSetMap(*_graph);

      }

      if (!_out_res_graph) {

        _out_res_graph = new OutResGraph(*_graph, *_capacity, *_preflow);

      }

      if (!_out_current_edge) {

        _out_current_edge = new OutCurrentEdgeMap(*_graph);

      }

      if (!_rev_graph) {

        _rev_graph = new RevGraph(*_graph);

      }

      if (!_in_res_graph) {

        _in_res_graph = new InResGraph(*_rev_graph, *_capacity, *_preflow);

      }

      if (!_in_current_edge) {

        _in_current_edge = new InCurrentEdgeMap(*_graph);

      }

      if (!_min_cut_map) {

        _min_cut_map = new MinCutMap(*_graph);

      }

      _min_cut = std::numeric_limits<Value>::max();

    }

    /// \brief Calculates the minimum cut with the \c source node

    /// in the \f$ V_{out} \f$.

    ///

    /// Calculates the minimum cut with the \c source node

    /// in the \f$ V_{out} \f$.

    void calculateOut() {

      for (NodeIt it(*_graph); it != INVALID; ++it) {

        if (it != _source) {

          calculateOut(it);

          return;

        }

      }

    }

    /// \brief Calculates the minimum cut with the \c source node

    /// in the \f$ V_{out} \f$.

    ///

    /// Calculates the minimum cut with the \c source node

    /// in the \f$ V_{out} \f$. The \c target is the initial target

    /// for the push-relabel algorithm.

    void calculateOut(const Node& target) {

      findMinCut(target, true, *_out_res_graph, *_out_current_edge);

    }

    /// \brief Calculates the minimum cut with the \c source node

    /// in the \f$ V_{in} \f$.

    ///

    /// Calculates the minimum cut with the \c source node

    /// in the \f$ V_{in} \f$.

    void calculateIn() {

      for (NodeIt it(*_graph); it != INVALID; ++it) {

        if (it != _source) {

          calculateIn(it);

          return;

        }

      }

    }

    /// \brief Calculates the minimum cut with the \c source node

    /// in the \f$ V_{in} \f$.

    ///

    /// Calculates the minimum cut with the \c source node

    /// in the \f$ V_{in} \f$. The \c target is the initial target

    /// for the push-relabel algorithm.

    void calculateIn(const Node& target) {

      findMinCut(target, false, *_in_res_graph, *_in_current_edge);

    }

    /// \brief Runs the algorithm.

    ///

    /// Runs the algorithm. It finds a proper \c source and \c target

    /// and then calls the \ref init(), \ref calculateOut() and \ref

    /// calculateIn().

    void run() {

      init();

      for (NodeIt it(*_graph); it != INVALID; ++it) {

        if (it != _source) {

          calculateOut(it);

          calculateIn(it);

          return;

        }

      }

    }

    /// \brief Runs the algorithm.

    ///

    /// Runs the algorithm. It finds a proper \c target and then calls

    /// the \ref init(), \ref calculateOut() and \ref calculateIn().

    void run(const Node& s) {

      init(s);

      for (NodeIt it(*_graph); it != INVALID; ++it) {

        if (it != _source) {

          calculateOut(it);

          calculateIn(it);

          return;

        }

      }

    }

    /// \brief Runs the algorithm.

    ///

    /// Runs the algorithm. It just calls the \ref init() and then

    /// \ref calculateOut() and \ref calculateIn().

    void run(const Node& s, const Node& t) {

      init(s); 

      calculateOut(t);

      calculateIn(t);

    }

    /// @}

    /// \name Query Functions The result of the %HaoOrlin algorithm

    /// can be obtained using these functions.

    /// \n

    /// Before the use of these functions, either \ref run(), \ref

    /// calculateOut() or \ref calculateIn() must be called.

    /// @{

    /// \brief Returns the value of the minimum value cut.

    /// 

    /// Returns the value of the minimum value cut.

    Value minCut() const {

      return _min_cut;

    }

    /// \brief Returns a minimum value cut.

    ///

    /// Sets \c nodeMap to the characteristic vector of a minimum

    /// value cut. The nodes in \f$ V_{out} \f$ will be set true and

    /// the nodes in \f$ V_{in} \f$ will be set false. 

    /// \pre nodeMap should be a bool-valued node-map.

    template <typename NodeMap>

    Value minCut(NodeMap& nodeMap) const {

      for (NodeIt it(*_graph); it != INVALID; ++it) {

	nodeMap.set(it, (*_min_cut_map)[it]);

      }

      return minCut();

    }

    /// @}

  }; //class HaoOrlin 

} //namespace lemon

#endif //LEMON_HAO_ORLIN_H