/* -*- 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 to find a minimum cut in directed graphs.
  ///
  /// Hao-Orlin calculates a minimum cut in a directed graph \f$
  /// D=(V,A) \f$. It takes a fixed node \f$ source \in V \f$ and consists
  /// of two phases: in the first phase it determines a minimum cut
  /// with \f$ source \f$ on the source-side (i.e. a set \f$ X\subsetneq V
  /// \f$ with \f$ source \in X \f$ and minimal out-degree) and in the
  /// second phase it determines a minimum cut with \f$ source \f$ on the
  /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \notin X \f$
  /// and minimal out-degree). Obviously, the smaller of these two
  /// cuts will be a minimum cut of \f$ D \f$. The algorithm is a
  /// modified push-relabel preflow algorithm and our implementation
  /// calculates the minimum cut in \f$ O(n^3) \f$ time (we use the
  /// highest-label rule). The purpose of such an algorithm is testing
  /// network reliability. For an undirected graph with \f$ n \f$
  /// nodes and \f$ e \f$ edges you can use the algorithm of Nagamochi
  /// and Ibaraki which solves the undirected problem in \f$ O(ne +
  /// n^2 \log(n)) \f$ time: 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 this 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. Node \c source  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 a minimum cut with \f$ source \f$ on the
    /// source-side.
    ///
    /// \brief Calculates a minimum cut with \f$ source \f$ on the
    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X
    /// \f$ and minimal out-degree).
    void calculateOut() {
      for (NodeIt it(*_graph); it != INVALID; ++it) {
        if (it != _source) {
          calculateOut(it);
          return;
        }
      }
    }

    /// \brief Calculates a minimum cut with \f$ source \f$ on the
    /// source-side.
    ///
    /// \brief Calculates a minimum cut with \f$ source \f$ on the
    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X
    /// \f$ and minimal out-degree). 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 a minimum cut with \f$ source \f$ on the
    /// sink-side.
    ///
    /// \brief Calculates a minimum cut with \f$ source \f$ on the
    /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \notin X
    /// \f$ and minimal out-degree).
    void calculateIn() {
      for (NodeIt it(*_graph); it != INVALID; ++it) {
        if (it != _source) {
          calculateIn(it);
          return;
        }
      }
    }

    /// \brief Calculates a minimum cut with \f$ source \f$ on the
    /// sink-side.
    ///
    /// \brief Calculates a minimum cut with \f$ source \f$ on the
    /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \notin
    /// X \f$ and minimal out-degree).  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 nodes \c source and \c target
    /// arbitrarily and then calls \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 uses the given \c source node, 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 cut.
    ///
    /// Sets \c nodeMap to the characteristic vector of a minimum
    /// value cut: it will give a nonempty set \f$ X\subsetneq V \f$
    /// with minimal out-degree (i.e. \c nodeMap will be true exactly
    /// for the nodes of \f$ X \f$.  \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