deba@2211: /* -*- C++ -*-
deba@2211:  *
deba@2225:  * This file is a part of LEMON, a generic C++ optimization library
deba@2225:  *
deba@2225:  * Copyright (C) 2003-2006
deba@2225:  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
deba@2211:  * (Egervary Research Group on Combinatorial Optimization, EGRES).
deba@2211:  *
deba@2211:  * Permission to use, modify and distribute this software is granted
deba@2211:  * provided that this copyright notice appears in all copies. For
deba@2211:  * precise terms see the accompanying LICENSE file.
deba@2211:  *
deba@2211:  * This software is provided "AS IS" with no warranty of any kind,
deba@2211:  * express or implied, and with no claim as to its suitability for any
deba@2211:  * purpose.
deba@2211:  *
deba@2211:  */
deba@2211: 
deba@2211: #ifndef LEMON_HAO_ORLIN_H
deba@2211: #define LEMON_HAO_ORLIN_H
deba@2211: 
deba@2211: #include <vector>
deba@2211: #include <queue>
deba@2211: #include <limits>
deba@2211: 
deba@2211: #include <lemon/maps.h>
deba@2211: #include <lemon/graph_utils.h>
deba@2211: #include <lemon/graph_adaptor.h>
deba@2211: #include <lemon/iterable_maps.h>
deba@2211:  
deba@2211: 
deba@2211: /// \file
deba@2211: /// \ingroup flowalgs
deba@2225: /// \brief Implementation of the Hao-Orlin algorithm.
deba@2225: ///
deba@2225: /// Implementation of the HaoOrlin algorithms class for testing network 
deba@2211: /// reliability.
deba@2211: 
deba@2211: namespace lemon {
deba@2211: 
deba@2225:   /// \ingroup flowalgs
deba@2225:   ///
athos@2228:   /// \brief %Hao-Orlin algorithm to find a minimum cut in directed graphs.
deba@2211:   ///
athos@2273:   /// Hao-Orlin calculates a minimum cut in a directed graph 
athos@2273:   /// \f$ D=(V,A) \f$. It takes a fixed node \f$ source \in V \f$ and consists
athos@2228:   /// of two phases: in the first phase it determines a minimum cut
athos@2273:   /// with \f$ source \f$ on the source-side (i.e. a set \f$ X\subsetneq V \f$
athos@2273:   /// with \f$ source \in X \f$ and minimal out-degree) and in the
athos@2228:   /// second phase it determines a minimum cut with \f$ source \f$ on the
athos@2228:   /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \notin X \f$
athos@2228:   /// and minimal out-degree). Obviously, the smaller of these two
athos@2228:   /// cuts will be a minimum cut of \f$ D \f$. The algorithm is a
athos@2228:   /// modified push-relabel preflow algorithm and our implementation
athos@2228:   /// calculates the minimum cut in \f$ O(n^3) \f$ time (we use the
athos@2228:   /// highest-label rule). The purpose of such an algorithm is testing
athos@2228:   /// network reliability. For an undirected graph with \f$ n \f$
athos@2228:   /// nodes and \f$ e \f$ edges you can use the algorithm of Nagamochi
athos@2273:   /// and Ibaraki which solves the undirected problem in 
athos@2273:   /// \f$ O(ne + n^2 \log(n)) \f$ time: it is implemented in the MinCut 
athos@2273:   /// algorithm
athos@2228:   /// class.
deba@2225:   ///
deba@2225:   /// \param _Graph is the graph type of the algorithm.
deba@2225:   /// \param _CapacityMap is an edge map of capacities which should
deba@2225:   /// be any numreric type. The default type is _Graph::EdgeMap<int>.
deba@2225:   /// \param _Tolerance is the handler of the inexact computation. The
athos@2228:   /// default type for this is Tolerance<typename CapacityMap::Value>.
deba@2211:   ///
deba@2211:   /// \author Attila Bernath and Balazs Dezso
deba@2225: #ifdef DOXYGEN
deba@2225:   template <typename _Graph, typename _CapacityMap, typename _Tolerance>
deba@2225: #else
deba@2211:   template <typename _Graph,
deba@2211: 	    typename _CapacityMap = typename _Graph::template EdgeMap<int>,
deba@2211:             typename _Tolerance = Tolerance<typename _CapacityMap::Value> >
deba@2225: #endif
deba@2211:   class HaoOrlin {
deba@2211:   protected:
deba@2211: 
deba@2211:     typedef _Graph Graph;
deba@2211:     typedef _CapacityMap CapacityMap;
deba@2211:     typedef _Tolerance Tolerance;
deba@2211: 
deba@2211:     typedef typename CapacityMap::Value Value;
deba@2211: 
deba@2211:     
deba@2211:     typedef typename Graph::Node Node;
deba@2211:     typedef typename Graph::NodeIt NodeIt;
deba@2211:     typedef typename Graph::EdgeIt EdgeIt;
deba@2211:     typedef typename Graph::OutEdgeIt OutEdgeIt;
deba@2211:     typedef typename Graph::InEdgeIt InEdgeIt;
deba@2211: 
deba@2211:     const Graph* _graph;
deba@2225: 
deba@2211:     const CapacityMap* _capacity;
deba@2211: 
deba@2211:     typedef typename Graph::template EdgeMap<Value> FlowMap;
deba@2211: 
deba@2211:     FlowMap* _preflow;
deba@2211: 
deba@2211:     Node _source, _target;
deba@2211:     int _node_num;
deba@2211: 
deba@2211:     typedef ResGraphAdaptor<const Graph, Value, CapacityMap, 
deba@2225:                             FlowMap, Tolerance> OutResGraph;
deba@2225:     typedef typename OutResGraph::Edge OutResEdge;
deba@2225:     
deba@2225:     OutResGraph* _out_res_graph;
deba@2211: 
deba@2225:     typedef typename Graph::template NodeMap<OutResEdge> OutCurrentEdgeMap;
deba@2225:     OutCurrentEdgeMap* _out_current_edge;  
deba@2211: 
deba@2225:     typedef RevGraphAdaptor<const Graph> RevGraph;
deba@2225:     RevGraph* _rev_graph;
deba@2225: 
deba@2225:     typedef ResGraphAdaptor<const RevGraph, Value, CapacityMap, 
deba@2225:                             FlowMap, Tolerance> InResGraph;
deba@2225:     typedef typename InResGraph::Edge InResEdge;
deba@2225:     
deba@2225:     InResGraph* _in_res_graph;
deba@2225: 
deba@2225:     typedef typename Graph::template NodeMap<InResEdge> InCurrentEdgeMap;
deba@2225:     InCurrentEdgeMap* _in_current_edge;  
deba@2211: 
deba@2211: 
deba@2211:     typedef IterableBoolMap<Graph, Node> WakeMap;
deba@2211:     WakeMap* _wake;
deba@2211: 
deba@2211:     typedef typename Graph::template NodeMap<int> DistMap;
deba@2211:     DistMap* _dist;  
deba@2211:     
deba@2211:     typedef typename Graph::template NodeMap<Value> ExcessMap;
deba@2211:     ExcessMap* _excess;
deba@2211: 
deba@2211:     typedef typename Graph::template NodeMap<bool> SourceSetMap;
deba@2211:     SourceSetMap* _source_set;
deba@2211: 
deba@2211:     std::vector<int> _level_size;
deba@2211: 
deba@2211:     int _highest_active;
deba@2211:     std::vector<std::list<Node> > _active_nodes;
deba@2211: 
deba@2211:     int _dormant_max;
deba@2211:     std::vector<std::list<Node> > _dormant;
deba@2211: 
deba@2211: 
deba@2211:     Value _min_cut;
deba@2211: 
deba@2211:     typedef typename Graph::template NodeMap<bool> MinCutMap;
deba@2211:     MinCutMap* _min_cut_map;
deba@2211: 
deba@2211:     Tolerance _tolerance;
deba@2211: 
deba@2211:   public: 
deba@2211: 
deba@2225:     /// \brief Constructor
deba@2225:     ///
deba@2225:     /// Constructor of the algorithm class. 
deba@2211:     HaoOrlin(const Graph& graph, const CapacityMap& capacity, 
deba@2211:              const Tolerance& tolerance = Tolerance()) :
deba@2211:       _graph(&graph), _capacity(&capacity), 
deba@2225:       _preflow(0), _source(), _target(), 
deba@2225:       _out_res_graph(0), _out_current_edge(0),
deba@2225:       _rev_graph(0), _in_res_graph(0), _in_current_edge(0),
deba@2211:       _wake(0),_dist(0), _excess(0), _source_set(0), 
deba@2211:       _highest_active(), _active_nodes(), _dormant_max(), _dormant(), 
deba@2211:       _min_cut(), _min_cut_map(0), _tolerance(tolerance) {}
deba@2211: 
deba@2211:     ~HaoOrlin() {
deba@2211:       if (_min_cut_map) {
deba@2211:         delete _min_cut_map;
deba@2211:       } 
deba@2225:       if (_in_current_edge) {
deba@2225:         delete _in_current_edge;
deba@2225:       }
deba@2225:       if (_in_res_graph) {
deba@2225:         delete _in_res_graph;
deba@2225:       }
deba@2225:       if (_rev_graph) {
deba@2225:         delete _rev_graph;
deba@2225:       }
deba@2225:       if (_out_current_edge) {
deba@2225:         delete _out_current_edge;
deba@2225:       }
deba@2225:       if (_out_res_graph) {
deba@2225:         delete _out_res_graph;
deba@2211:       }
deba@2211:       if (_source_set) {
deba@2211:         delete _source_set;
deba@2211:       }
deba@2211:       if (_excess) {
deba@2211:         delete _excess;
deba@2211:       }
deba@2211:       if (_dist) {
deba@2211:         delete _dist;
deba@2211:       }
deba@2211:       if (_wake) {
deba@2211:         delete _wake;
deba@2211:       }
deba@2211:       if (_preflow) {
deba@2211:         delete _preflow;
deba@2211:       }
deba@2211:     }
deba@2211:     
deba@2211:   private:
deba@2211:     
deba@2225:     template <typename ResGraph, typename EdgeMap>
deba@2225:     void findMinCut(const Node& target, bool out, 
deba@2225:                     ResGraph& res_graph, EdgeMap& current_edge) {
deba@2225:       typedef typename ResGraph::Edge ResEdge;
deba@2225:       typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
deba@2225: 
deba@2225:       for (typename Graph::EdgeIt it(*_graph); it != INVALID; ++it) {
deba@2225:         (*_preflow)[it] = 0;      
deba@2225:       }
deba@2225:       for (NodeIt it(*_graph); it != INVALID; ++it) {
deba@2225:         (*_wake)[it] = true;
deba@2225:         (*_dist)[it] = 1;
deba@2225:         (*_excess)[it] = 0;
deba@2225:         (*_source_set)[it] = false;
deba@2225: 
deba@2225:         res_graph.firstOut(current_edge[it], it);
deba@2225:       }
deba@2225: 
deba@2225:       _target = target;
deba@2225:       (*_dist)[target] = 0;
deba@2225: 
deba@2225:       for (ResOutEdgeIt it(res_graph, _source); it != INVALID; ++it) {
deba@2225:         Value delta = res_graph.rescap(it);
deba@2225:         if (!_tolerance.positive(delta)) continue;
deba@2225:         
deba@2225:         (*_excess)[res_graph.source(it)] -= delta;
deba@2225:         res_graph.augment(it, delta);
deba@2225:         Node a = res_graph.target(it);
deba@2225:         if (!_tolerance.positive((*_excess)[a]) && 
deba@2225:             (*_wake)[a] && a != _target) {
deba@2225:           _active_nodes[(*_dist)[a]].push_front(a);
deba@2225:           if (_highest_active < (*_dist)[a]) {
deba@2225:             _highest_active = (*_dist)[a];
deba@2225:           }
deba@2225:         }
deba@2225:         (*_excess)[a] += delta;
deba@2225:       }
deba@2225: 
deba@2225:       _dormant[0].push_front(_source);
deba@2225:       (*_source_set)[_source] = true;
deba@2225:       _dormant_max = 0;
deba@2225:       (*_wake)[_source] = false;
deba@2225: 
deba@2225:       _level_size[0] = 1;
deba@2225:       _level_size[1] = _node_num - 1;
deba@2225: 
deba@2225:       do {
deba@2225: 	Node n;
deba@2225: 	while ((n = findActiveNode()) != INVALID) {
deba@2225: 	  ResEdge e;
deba@2225: 	  while (_tolerance.positive((*_excess)[n]) && 
deba@2225:                  (e = findAdmissibleEdge(n, res_graph, current_edge)) 
deba@2225:                  != INVALID){
deba@2225: 	    Value delta;
deba@2225: 	    if ((*_excess)[n] < res_graph.rescap(e)) {
deba@2225: 	      delta = (*_excess)[n];
deba@2225: 	    } else {
deba@2225: 	      delta = res_graph.rescap(e);
deba@2225: 	      res_graph.nextOut(current_edge[n]);
deba@2225: 	    }
deba@2225:             if (!_tolerance.positive(delta)) continue;
deba@2225: 	    res_graph.augment(e, delta);
deba@2225: 	    (*_excess)[res_graph.source(e)] -= delta;
deba@2225: 	    Node a = res_graph.target(e);
deba@2225: 	    if (!_tolerance.positive((*_excess)[a]) && a != _target) {
deba@2225: 	      _active_nodes[(*_dist)[a]].push_front(a);
deba@2225: 	    }
deba@2225: 	    (*_excess)[a] += delta;
deba@2225: 	  }
deba@2225: 	  if (_tolerance.positive((*_excess)[n])) {
deba@2225: 	    relabel(n, res_graph, current_edge);
deba@2225:           }
deba@2225: 	}
deba@2225: 
deba@2225: 	Value current_value = cutValue(out);
deba@2225:  	if (_min_cut > current_value){
deba@2225:           if (out) {
deba@2225:             for (NodeIt it(*_graph); it != INVALID; ++it) {
deba@2225:               _min_cut_map->set(it, !(*_wake)[it]);
deba@2225:             } 
deba@2225:           } else {
deba@2225:             for (NodeIt it(*_graph); it != INVALID; ++it) {
deba@2225:               _min_cut_map->set(it, (*_wake)[it]);
deba@2225:             } 
deba@2225:           }
deba@2225: 
deba@2225: 	  _min_cut = current_value;
deba@2225:  	}
deba@2225: 
deba@2225:       } while (selectNewSink(res_graph));
deba@2225:     }
deba@2225: 
deba@2225:     template <typename ResGraph, typename EdgeMap>
deba@2225:     void relabel(const Node& n, ResGraph& res_graph, EdgeMap& current_edge) {
deba@2225:       typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
deba@2225: 
deba@2225:       int k = (*_dist)[n];
deba@2211:       if (_level_size[k] == 1) {
deba@2211: 	++_dormant_max;
deba@2211: 	for (NodeIt it(*_graph); it != INVALID; ++it) {
deba@2211: 	  if ((*_wake)[it] && (*_dist)[it] >= k) {
deba@2211: 	    (*_wake)[it] = false;
deba@2211: 	    _dormant[_dormant_max].push_front(it);
deba@2211: 	    --_level_size[(*_dist)[it]];
deba@2211: 	  }
deba@2211: 	}
deba@2211: 	--_highest_active;
deba@2225:       } else {	
deba@2225:         int new_dist = _node_num;
deba@2225:         for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
deba@2225:           Node t = res_graph.target(e);
deba@2225:           if ((*_wake)[t] && new_dist > (*_dist)[t]) {
deba@2225:             new_dist = (*_dist)[t];
deba@2225:           }
deba@2225:         }
deba@2225:         if (new_dist == _node_num) {
deba@2211: 	  ++_dormant_max;
deba@2225: 	  (*_wake)[n] = false;
deba@2225: 	  _dormant[_dormant_max].push_front(n);
deba@2225: 	  --_level_size[(*_dist)[n]];
deba@2225: 	} else {	    
deba@2225: 	  --_level_size[(*_dist)[n]];
deba@2225: 	  (*_dist)[n] = new_dist + 1;
deba@2225: 	  _highest_active = (*_dist)[n];
deba@2225: 	  _active_nodes[_highest_active].push_front(n);
deba@2225: 	  ++_level_size[(*_dist)[n]];
deba@2225: 	  res_graph.firstOut(current_edge[n], n);
deba@2211: 	}
deba@2211:       }
deba@2211:     }
deba@2211: 
deba@2225:     template <typename ResGraph>
deba@2225:     bool selectNewSink(ResGraph& res_graph) {
deba@2225:       typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
deba@2225: 
deba@2211:       Node old_target = _target;
deba@2211:       (*_wake)[_target] = false;
deba@2211:       --_level_size[(*_dist)[_target]];
deba@2211:       _dormant[0].push_front(_target);
deba@2211:       (*_source_set)[_target] = true;
deba@2211:       if ((int)_dormant[0].size() == _node_num){
deba@2211:         _dormant[0].clear();
deba@2211: 	return false;
deba@2211:       }
deba@2211: 
deba@2211:       bool wake_was_empty = false;
deba@2211: 
deba@2211:       if(_wake->trueNum() == 0) {
deba@2211: 	while (!_dormant[_dormant_max].empty()){
deba@2211: 	  (*_wake)[_dormant[_dormant_max].front()] = true;
deba@2211: 	  ++_level_size[(*_dist)[_dormant[_dormant_max].front()]];
deba@2211: 	  _dormant[_dormant_max].pop_front();
deba@2211: 	}
deba@2211: 	--_dormant_max;
deba@2211: 	wake_was_empty = true;
deba@2211:       }
deba@2211: 
deba@2211:       int min_dist = std::numeric_limits<int>::max();
deba@2211:       for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
deba@2211: 	if (min_dist > (*_dist)[it]){
deba@2211: 	  _target = it;
deba@2211: 	  min_dist = (*_dist)[it];
deba@2211: 	}
deba@2211:       }
deba@2211: 
deba@2211:       if (wake_was_empty){
deba@2211: 	for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
deba@2211:           if (_tolerance.positive((*_excess)[it])) {
deba@2211: 	    if ((*_wake)[it] && it != _target) {
deba@2211: 	      _active_nodes[(*_dist)[it]].push_front(it);
deba@2211:             }
deba@2211: 	    if (_highest_active < (*_dist)[it]) {
deba@2211: 	      _highest_active = (*_dist)[it];		    
deba@2211:             }
deba@2211: 	  }
deba@2211: 	}
deba@2211:       }
deba@2211: 
deba@2225:       for (ResOutEdgeIt e(res_graph, old_target); e!=INVALID; ++e){
deba@2225: 	if (!(*_source_set)[res_graph.target(e)]) {
deba@2225:           Value delta = res_graph.rescap(e);
deba@2225:           if (!_tolerance.positive(delta)) continue;
deba@2225:           res_graph.augment(e, delta);
deba@2225:           (*_excess)[res_graph.source(e)] -= delta;
deba@2225:           Node a = res_graph.target(e);
deba@2225:           if (!_tolerance.positive((*_excess)[a]) && 
deba@2225:               (*_wake)[a] && a != _target) {
deba@2225:             _active_nodes[(*_dist)[a]].push_front(a);
deba@2225:             if (_highest_active < (*_dist)[a]) {
deba@2225:               _highest_active = (*_dist)[a];
deba@2225:             }
deba@2225:           }
deba@2225:           (*_excess)[a] += delta;
deba@2211: 	}
deba@2211:       }
deba@2211:       
deba@2211:       return true;
deba@2211:     }
deba@2211:     
deba@2211:     Node findActiveNode() {
deba@2211:       while (_highest_active > 0 && _active_nodes[_highest_active].empty()){ 
deba@2211: 	--_highest_active;
deba@2211:       }
deba@2211:       if( _highest_active > 0) {
deba@2211:        	Node n = _active_nodes[_highest_active].front();
deba@2211: 	_active_nodes[_highest_active].pop_front();
deba@2211: 	return n;
deba@2211:       } else {
deba@2211: 	return INVALID;
deba@2211:       }
deba@2211:     }
deba@2211: 
deba@2225:     template <typename ResGraph, typename EdgeMap>
deba@2225:     typename ResGraph::Edge findAdmissibleEdge(const Node& n, 
deba@2225:                                                ResGraph& res_graph, 
deba@2225:                                                EdgeMap& current_edge) {
deba@2225:       typedef typename ResGraph::Edge ResEdge;
deba@2225:       ResEdge e = current_edge[n];
deba@2211:       while (e != INVALID && 
deba@2225:              ((*_dist)[n] <= (*_dist)[res_graph.target(e)] || 
deba@2225:               !(*_wake)[res_graph.target(e)])) {
deba@2225: 	res_graph.nextOut(e);
deba@2211:       }
deba@2211:       if (e != INVALID) {
deba@2225: 	current_edge[n] = e;	
deba@2211: 	return e;
deba@2211:       } else {
deba@2211: 	return INVALID;
deba@2211:       }
deba@2211:     }
deba@2211: 
deba@2225:     Value cutValue(bool out) {
deba@2225:       Value value = 0;
deba@2225:       if (out) {
deba@2225:         for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
deba@2225:           for (InEdgeIt e(*_graph, it); e != INVALID; ++e) {
deba@2225:             if (!(*_wake)[_graph->source(e)]){
deba@2225:               value += (*_capacity)[e];
deba@2225:             }	
deba@2225:           }
deba@2225:         }
deba@2225:       } else {
deba@2225:         for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
deba@2225:           for (OutEdgeIt e(*_graph, it); e != INVALID; ++e) {
deba@2225:             if (!(*_wake)[_graph->target(e)]){
deba@2225:               value += (*_capacity)[e];
deba@2225:             }	
deba@2225:           }
deba@2211:         }
deba@2211:       }
deba@2225:       return value;
deba@2211:     }
deba@2225: 
deba@2211: 
deba@2211:   public:
deba@2211: 
deba@2225:     /// \name Execution control
deba@2225:     /// The simplest way to execute the algorithm is to use
deba@2225:     /// one of the member functions called \c run(...).
deba@2225:     /// \n
deba@2225:     /// If you need more control on the execution,
deba@2225:     /// first you must call \ref init(), then the \ref calculateIn() or
deba@2225:     /// \ref calculateIn() functions.
deba@2225: 
deba@2225:     /// @{
deba@2225: 
deba@2211:     /// \brief Initializes the internal data structures.
deba@2211:     ///
deba@2211:     /// Initializes the internal data structures. It creates
deba@2225:     /// the maps, residual graph adaptors and some bucket structures
deba@2211:     /// for the algorithm. 
deba@2211:     void init() {
deba@2211:       init(NodeIt(*_graph));
deba@2211:     }
deba@2211: 
deba@2211:     /// \brief Initializes the internal data structures.
deba@2211:     ///
deba@2211:     /// Initializes the internal data structures. It creates
deba@2211:     /// the maps, residual graph adaptor and some bucket structures
athos@2228:     /// for the algorithm. Node \c source  is used as the push-relabel
deba@2211:     /// algorithm's source.
deba@2211:     void init(const Node& source) {
deba@2211:       _source = source;
deba@2211:       _node_num = countNodes(*_graph);
deba@2211: 
deba@2211:       _dormant.resize(_node_num);
deba@2211:       _level_size.resize(_node_num, 0);
deba@2211:       _active_nodes.resize(_node_num);
deba@2211: 
deba@2211:       if (!_preflow) {
deba@2211:         _preflow = new FlowMap(*_graph);
deba@2211:       }
deba@2211:       if (!_wake) {
deba@2211:         _wake = new WakeMap(*_graph);
deba@2211:       }
deba@2211:       if (!_dist) {
deba@2211:         _dist = new DistMap(*_graph);
deba@2211:       }
deba@2211:       if (!_excess) {
deba@2211:         _excess = new ExcessMap(*_graph);
deba@2211:       }
deba@2211:       if (!_source_set) {
deba@2211:         _source_set = new SourceSetMap(*_graph);
deba@2211:       }
deba@2225:       if (!_out_res_graph) {
deba@2225:         _out_res_graph = new OutResGraph(*_graph, *_capacity, *_preflow);
deba@2225:       }
deba@2225:       if (!_out_current_edge) {
deba@2225:         _out_current_edge = new OutCurrentEdgeMap(*_graph);
deba@2225:       }
deba@2225:       if (!_rev_graph) {
deba@2225:         _rev_graph = new RevGraph(*_graph);
deba@2225:       }
deba@2225:       if (!_in_res_graph) {
deba@2225:         _in_res_graph = new InResGraph(*_rev_graph, *_capacity, *_preflow);
deba@2225:       }
deba@2225:       if (!_in_current_edge) {
deba@2225:         _in_current_edge = new InCurrentEdgeMap(*_graph);
deba@2211:       }
deba@2211:       if (!_min_cut_map) {
deba@2211:         _min_cut_map = new MinCutMap(*_graph);
deba@2211:       }
deba@2211: 
deba@2211:       _min_cut = std::numeric_limits<Value>::max();
deba@2211:     }
deba@2211: 
deba@2211: 
athos@2228:     /// \brief Calculates a minimum cut with \f$ source \f$ on the
athos@2228:     /// source-side.
deba@2211:     ///
athos@2228:     /// \brief Calculates a minimum cut with \f$ source \f$ on the
athos@2273:     /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
athos@2273:     ///  and minimal out-degree).
deba@2211:     void calculateOut() {
deba@2211:       for (NodeIt it(*_graph); it != INVALID; ++it) {
deba@2211:         if (it != _source) {
deba@2211:           calculateOut(it);
deba@2211:           return;
deba@2211:         }
deba@2211:       }
deba@2211:     }
deba@2211: 
athos@2228:     /// \brief Calculates a minimum cut with \f$ source \f$ on the
athos@2228:     /// source-side.
deba@2211:     ///
athos@2228:     /// \brief Calculates a minimum cut with \f$ source \f$ on the
athos@2273:     /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
athos@2273:     ///  and minimal out-degree). The \c target is the initial target
deba@2211:     /// for the push-relabel algorithm.
deba@2211:     void calculateOut(const Node& target) {
deba@2225:       findMinCut(target, true, *_out_res_graph, *_out_current_edge);
deba@2211:     }
deba@2211: 
athos@2228:     /// \brief Calculates a minimum cut with \f$ source \f$ on the
athos@2228:     /// sink-side.
deba@2225:     ///
athos@2228:     /// \brief Calculates a minimum cut with \f$ source \f$ on the
athos@2273:     /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with 
athos@2273:     /// \f$ source \notin X \f$
athos@2273:     /// and minimal out-degree).
deba@2211:     void calculateIn() {
deba@2211:       for (NodeIt it(*_graph); it != INVALID; ++it) {
deba@2211:         if (it != _source) {
deba@2211:           calculateIn(it);
deba@2211:           return;
deba@2211:         }
deba@2211:       }
deba@2211:     }
deba@2211: 
athos@2228:     /// \brief Calculates a minimum cut with \f$ source \f$ on the
athos@2228:     /// sink-side.
deba@2225:     ///
athos@2228:     /// \brief Calculates a minimum cut with \f$ source \f$ on the
athos@2273:     /// sink-side (i.e. a set \f$ X\subsetneq V 
athos@2273:     /// \f$ with \f$ source \notin X \f$ and minimal out-degree).  
athos@2273:     /// The \c target is the initial
athos@2228:     /// target for the push-relabel algorithm.
deba@2225:     void calculateIn(const Node& target) {
deba@2225:       findMinCut(target, false, *_in_res_graph, *_in_current_edge);
deba@2225:     }
deba@2225: 
deba@2225:     /// \brief Runs the algorithm.
deba@2225:     ///
athos@2228:     /// Runs the algorithm. It finds nodes \c source and \c target
athos@2228:     /// arbitrarily and then calls \ref init(), \ref calculateOut()
athos@2228:     /// and \ref calculateIn().
deba@2211:     void run() {
deba@2211:       init();
deba@2211:       for (NodeIt it(*_graph); it != INVALID; ++it) {
deba@2211:         if (it != _source) {
deba@2225:           calculateOut(it);
deba@2225:           calculateIn(it);
deba@2211:           return;
deba@2211:         }
deba@2211:       }
deba@2211:     }
deba@2211: 
deba@2225:     /// \brief Runs the algorithm.
deba@2225:     ///
athos@2228:     /// Runs the algorithm. It uses the given \c source node, finds a
athos@2228:     /// proper \c target and then calls the \ref init(), \ref
athos@2228:     /// calculateOut() and \ref calculateIn().
deba@2211:     void run(const Node& s) {
deba@2211:       init(s);
deba@2211:       for (NodeIt it(*_graph); it != INVALID; ++it) {
deba@2211:         if (it != _source) {
deba@2225:           calculateOut(it);
deba@2225:           calculateIn(it);
deba@2211:           return;
deba@2211:         }
deba@2211:       }
deba@2211:     }
deba@2211: 
deba@2225:     /// \brief Runs the algorithm.
deba@2225:     ///
deba@2225:     /// Runs the algorithm. It just calls the \ref init() and then
deba@2225:     /// \ref calculateOut() and \ref calculateIn().
deba@2211:     void run(const Node& s, const Node& t) {
deba@2225:       init(s); 
deba@2225:       calculateOut(t);
deba@2225:       calculateIn(t);
deba@2211:     }
deba@2225: 
deba@2225:     /// @}
deba@2211:     
athos@2275:     /// \name Query Functions 
athos@2275:     /// The result of the %HaoOrlin algorithm
deba@2225:     /// can be obtained using these functions.
deba@2225:     /// \n
athos@2275:     /// Before using these functions, either \ref run(), \ref
deba@2225:     /// calculateOut() or \ref calculateIn() must be called.
deba@2225:     
deba@2225:     /// @{
deba@2225: 
deba@2225:     /// \brief Returns the value of the minimum value cut.
deba@2211:     /// 
deba@2225:     /// Returns the value of the minimum value cut.
deba@2211:     Value minCut() const {
deba@2211:       return _min_cut;
deba@2211:     }
deba@2211: 
deba@2211: 
athos@2228:     /// \brief Returns a minimum cut.
deba@2211:     ///
deba@2211:     /// Sets \c nodeMap to the characteristic vector of a minimum
athos@2228:     /// value cut: it will give a nonempty set \f$ X\subsetneq V \f$
athos@2228:     /// with minimal out-degree (i.e. \c nodeMap will be true exactly
athos@2275:     /// for the nodes of \f$ X \f$).  \pre nodeMap should be a
athos@2228:     /// bool-valued node-map.
deba@2211:     template <typename NodeMap>
deba@2211:     Value minCut(NodeMap& nodeMap) const {
deba@2211:       for (NodeIt it(*_graph); it != INVALID; ++it) {
deba@2211: 	nodeMap.set(it, (*_min_cut_map)[it]);
deba@2211:       }
deba@2211:       return minCut();
deba@2211:     }
deba@2225: 
deba@2225:     /// @}
deba@2211:     
deba@2211:   }; //class HaoOrlin 
deba@2211: 
deba@2211: 
deba@2211: } //namespace lemon
deba@2211: 
deba@2211: #endif //LEMON_HAO_ORLIN_H