* Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

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

  *

  * Copyright (C) 2003-2006

  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

 /// Implementation of the Hao-Orlin algorithms class for testing network 

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

 ///

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

   /// \addtogroup flowalgs

   /// \ingroup flowalgs

   /// @{                                                   

   ///

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

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

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

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

   /// summary of the edge capacities go from the first set to the

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

   /// second set is the minimum.  The algorithm is a modified

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

   /// push-relabel preflow algorithm and it calculates the minimum cat

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

   /// in \f$ O(n^3) \f$ time. The purpose of such algorithm is testing

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

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

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

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

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

   /// problem in \f$ O(n^3) \f$ time. 

   /// 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>.

                             FlowMap, Tolerance> ResGraph;

                             FlowMap, Tolerance> OutResGraph;

     typedef typename ResGraph::Node ResNode;

     typedef typename OutResGraph::Edge OutResEdge;

     typedef typename ResGraph::NodeIt ResNodeIt;

     typedef typename ResGraph::EdgeIt ResEdgeIt;

     OutResGraph* _out_res_graph;

     typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;

     typedef typename ResGraph::Edge ResEdge;

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

     typedef typename ResGraph::InEdgeIt ResInEdgeIt;

     OutCurrentEdgeMap* _out_current_edge;  

     ResGraph* _res_graph;

     typedef RevGraphAdaptor<const Graph> RevGraph;

     RevGraph* _rev_graph;

     typedef typename Graph::template NodeMap<ResEdge> CurrentArcMap;

     CurrentArcMap* _current_arc;  

     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;  

       _preflow(0), _source(), _target(), _res_graph(0), _current_arc(0),

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

       _out_res_graph(0), _out_current_edge(0),

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

       if (_current_arc) {

       if (_in_current_edge) {

         delete _current_arc;

         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;

     void relabel(Node i) {

     template <typename ResGraph, typename EdgeMap>

       int k = (*_dist)[i];

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

       } else {

       } else {	

 	ResOutEdgeIt e(*_res_graph, i);

         int new_dist = _node_num;

 	while (e != INVALID && !(*_wake)[_res_graph->target(e)]) {

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

 	  ++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);

       }

 	if (e == INVALID){

     }

 	  ++_dormant_max;

 	  (*_wake)[i] = false;

     template <typename ResGraph>

 	  _dormant[_dormant_max].push_front(i);

     bool selectNewSink(ResGraph& res_graph) {

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

       typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;

 	} else{	    

 	  Node j = _res_graph->target(e);

 	  int new_dist = (*_dist)[j];

 	  ++e;

 	  while (e != INVALID){

 	    Node j = _res_graph->target(e);

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

 	      new_dist = (*_dist)[j];

             }

 	    ++e;

 	  }

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

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

 	  _highest_active = (*_dist)[i];

 	  _active_nodes[_highest_active].push_front(i);

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

 	  _res_graph->firstOut((*_current_arc)[i], i);

 	}

       }

     }

     bool selectNewSink(){

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

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

 	if (!(*_source_set)[_res_graph->target(e)]){

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

 	  push(e, _res_graph->rescap(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;

     ResEdge findAdmissibleEdge(const Node& n){

     template <typename ResGraph, typename EdgeMap>

       ResEdge e = (*_current_arc)[n];

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

                                                ResGraph& res_graph, 

                                                EdgeMap& current_edge) {

       typedef typename ResGraph::Edge ResEdge;

       ResEdge e = current_edge[n];

              ((*_dist)[n] <= (*_dist)[_res_graph->target(e)] || 

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

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

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

 	_res_graph->nextOut(e);

 	res_graph.nextOut(e);

 	(*_current_arc)[n] = e;	

 	current_edge[n] = e;	

     void push(ResEdge& e,const Value& delta){

     Value cutValue(bool out) {

       _res_graph->augment(e, delta);

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

       (*_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;

     }

     Value cutValue() {

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

       if (out) {

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

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

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

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

 	    value += (*_capacity)[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];

             }	

           }

         }

     }    

     }

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

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

       if (!_current_arc) {

       if (!_out_res_graph) {

         _current_arc = new CurrentArcMap(*_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);

     /// in the first partition.

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

     /// in the first partition.

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

     /// in the first partition.

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

     /// in the first partition. The \c target is the initial target

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

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

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

         (*_wake)[it] = true;

     }

         (*_dist)[it] = 1;

         (*_excess)[it] = 0;

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

         (*_source_set)[it] = false;

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

     ///

         _res_graph->firstOut((*_current_arc)[it], it);

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

       }

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

       _target = target;

       (*_dist)[target] = 0;

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

 	push(it, _res_graph->rescap(it));

       }

       _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)) != INVALID){

 	    Value delta;

 	    if ((*_excess)[n] < _res_graph->rescap(e)) {

 	      delta = (*_excess)[n];

 	    } else {

 	      delta = _res_graph->rescap(e);

 	      _res_graph->nextOut((*_current_arc)[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);

           }

 	}

 	Value current_value = cutValue();

  	if (_min_cut > current_value){

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

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

 	  }

 	  _min_cut = current_value;

  	}

       } while (selectNewSink());

     }

           startOut(it);

           calculateOut(it);

           //          startIn(it);

           calculateIn(it);

           startOut(it);

           calculateOut(it);

           //          startIn(it);

           calculateIn(it);

       init(s);

       init(s); 

       startOut(t);

       calculateOut(t);

       startIn(t);

       calculateIn(t);

     }

     }

     /// \brief Returns the value of the minimum value cut with node \c

     /// @}

     /// source on the source side.

     /// \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 with node \c source

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

     /// on the source side. This function can be called after running

     /// \ref findMinCut.

     /// value cut with node \c source on the source side. This

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

     /// function can be called after running \ref findMinCut.  

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