lemon/hao_orlin.h
changeset 2375 e30a0fdad0d7
parent 2275 ff46747676ed
child 2376 0ed45a6c74b1
equal deleted inserted replaced
4:05e17bc2435d 5:26bac530db45
    17  */
    17  */
    18 
    18 
    19 #ifndef LEMON_HAO_ORLIN_H
    19 #ifndef LEMON_HAO_ORLIN_H
    20 #define LEMON_HAO_ORLIN_H
    20 #define LEMON_HAO_ORLIN_H
    21 
    21 
       
    22 #include <cassert>
       
    23  
       
    24 
       
    25 
    22 #include <vector>
    26 #include <vector>
    23 #include <queue>
    27 #include <queue>
       
    28 #include <list>
    24 #include <limits>
    29 #include <limits>
    25 
    30 
    26 #include <lemon/maps.h>
    31 #include <lemon/maps.h>
    27 #include <lemon/graph_utils.h>
    32 #include <lemon/graph_utils.h>
    28 #include <lemon/graph_adaptor.h>
    33 #include <lemon/graph_adaptor.h>
    29 #include <lemon/iterable_maps.h>
    34 #include <lemon/iterable_maps.h>
    30  
       
    31 
    35 
    32 /// \file
    36 /// \file
    33 /// \ingroup flowalgs
    37 /// \ingroup flowalgs
    34 /// \brief Implementation of the Hao-Orlin algorithm.
    38 /// \brief Implementation of the Hao-Orlin algorithm.
    35 ///
    39 ///
   106                             FlowMap, Tolerance> OutResGraph;
   110                             FlowMap, Tolerance> OutResGraph;
   107     typedef typename OutResGraph::Edge OutResEdge;
   111     typedef typename OutResGraph::Edge OutResEdge;
   108     
   112     
   109     OutResGraph* _out_res_graph;
   113     OutResGraph* _out_res_graph;
   110 
   114 
   111     typedef typename Graph::template NodeMap<OutResEdge> OutCurrentEdgeMap;
       
   112     OutCurrentEdgeMap* _out_current_edge;  
       
   113 
       
   114     typedef RevGraphAdaptor<const Graph> RevGraph;
   115     typedef RevGraphAdaptor<const Graph> RevGraph;
   115     RevGraph* _rev_graph;
   116     RevGraph* _rev_graph;
   116 
   117 
   117     typedef ResGraphAdaptor<const RevGraph, Value, CapacityMap, 
   118     typedef ResGraphAdaptor<const RevGraph, Value, CapacityMap, 
   118                             FlowMap, Tolerance> InResGraph;
   119                             FlowMap, Tolerance> InResGraph;
   119     typedef typename InResGraph::Edge InResEdge;
   120     typedef typename InResGraph::Edge InResEdge;
   120     
   121     
   121     InResGraph* _in_res_graph;
   122     InResGraph* _in_res_graph;
   122 
       
   123     typedef typename Graph::template NodeMap<InResEdge> InCurrentEdgeMap;
       
   124     InCurrentEdgeMap* _in_current_edge;  
       
   125 
       
   126 
   123 
   127     typedef IterableBoolMap<Graph, Node> WakeMap;
   124     typedef IterableBoolMap<Graph, Node> WakeMap;
   128     WakeMap* _wake;
   125     WakeMap* _wake;
   129 
   126 
   130     typedef typename Graph::template NodeMap<int> DistMap;
   127     typedef typename Graph::template NodeMap<int> DistMap;
   159     /// Constructor of the algorithm class. 
   156     /// Constructor of the algorithm class. 
   160     HaoOrlin(const Graph& graph, const CapacityMap& capacity, 
   157     HaoOrlin(const Graph& graph, const CapacityMap& capacity, 
   161              const Tolerance& tolerance = Tolerance()) :
   158              const Tolerance& tolerance = Tolerance()) :
   162       _graph(&graph), _capacity(&capacity), 
   159       _graph(&graph), _capacity(&capacity), 
   163       _preflow(0), _source(), _target(), 
   160       _preflow(0), _source(), _target(), 
   164       _out_res_graph(0), _out_current_edge(0),
   161       _out_res_graph(0), _rev_graph(0), _in_res_graph(0),
   165       _rev_graph(0), _in_res_graph(0), _in_current_edge(0),
       
   166       _wake(0),_dist(0), _excess(0), _source_set(0), 
   162       _wake(0),_dist(0), _excess(0), _source_set(0), 
   167       _highest_active(), _active_nodes(), _dormant_max(), _dormant(), 
   163       _highest_active(), _active_nodes(), _dormant_max(), _dormant(), 
   168       _min_cut(), _min_cut_map(0), _tolerance(tolerance) {}
   164       _min_cut(), _min_cut_map(0), _tolerance(tolerance) {}
   169 
   165 
   170     ~HaoOrlin() {
   166     ~HaoOrlin() {
   171       if (_min_cut_map) {
   167       if (_min_cut_map) {
   172         delete _min_cut_map;
   168         delete _min_cut_map;
   173       } 
   169       } 
   174       if (_in_current_edge) {
       
   175         delete _in_current_edge;
       
   176       }
       
   177       if (_in_res_graph) {
   170       if (_in_res_graph) {
   178         delete _in_res_graph;
   171         delete _in_res_graph;
   179       }
   172       }
   180       if (_rev_graph) {
   173       if (_rev_graph) {
   181         delete _rev_graph;
   174         delete _rev_graph;
   182       }
   175       }
   183       if (_out_current_edge) {
       
   184         delete _out_current_edge;
       
   185       }
       
   186       if (_out_res_graph) {
   176       if (_out_res_graph) {
   187         delete _out_res_graph;
   177         delete _out_res_graph;
   188       }
   178       }
   189       if (_source_set) {
   179       if (_source_set) {
   190         delete _source_set;
   180         delete _source_set;
   203       }
   193       }
   204     }
   194     }
   205     
   195     
   206   private:
   196   private:
   207     
   197     
   208     template <typename ResGraph, typename EdgeMap>
   198     template <typename ResGraph>
   209     void findMinCut(const Node& target, bool out, 
   199     void findMinCut(const Node& target, bool out, ResGraph& res_graph) {
   210                     ResGraph& res_graph, EdgeMap& current_edge) {
   200       typedef typename Graph::Node Node;
   211       typedef typename ResGraph::Edge ResEdge;
   201       typedef typename ResGraph::Edge ResEdge;
   212       typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
   202       typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
   213 
   203 
   214       for (typename Graph::EdgeIt it(*_graph); it != INVALID; ++it) {
   204       for (typename Graph::EdgeIt it(*_graph); it != INVALID; ++it) {
   215         (*_preflow)[it] = 0;      
   205         (*_preflow)[it] = 0;      
   217       for (NodeIt it(*_graph); it != INVALID; ++it) {
   207       for (NodeIt it(*_graph); it != INVALID; ++it) {
   218         (*_wake)[it] = true;
   208         (*_wake)[it] = true;
   219         (*_dist)[it] = 1;
   209         (*_dist)[it] = 1;
   220         (*_excess)[it] = 0;
   210         (*_excess)[it] = 0;
   221         (*_source_set)[it] = false;
   211         (*_source_set)[it] = false;
   222 
   212       }
   223         res_graph.firstOut(current_edge[it], it);
   213 
   224       }
   214       _dormant[0].push_front(_source);
       
   215       (*_source_set)[_source] = true;
       
   216       _dormant_max = 0;
       
   217       (*_wake)[_source] = false;
       
   218 
       
   219       _level_size[0] = 1;
       
   220       _level_size[1] = _node_num - 1;
   225 
   221 
   226       _target = target;
   222       _target = target;
   227       (*_dist)[target] = 0;
   223       (*_dist)[target] = 0;
   228 
   224 
   229       for (ResOutEdgeIt it(res_graph, _source); it != INVALID; ++it) {
   225       for (ResOutEdgeIt it(res_graph, _source); it != INVALID; ++it) {
   230         Value delta = res_graph.rescap(it);
   226         Value delta = res_graph.rescap(it);
   231         if (!_tolerance.positive(delta)) continue;
   227         (*_excess)[_source] -= delta;
   232         
       
   233         (*_excess)[res_graph.source(it)] -= delta;
       
   234         res_graph.augment(it, delta);
   228         res_graph.augment(it, delta);
   235         Node a = res_graph.target(it);
   229         Node a = res_graph.target(it);
   236         if (!_tolerance.positive((*_excess)[a]) && 
   230         if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) {
   237             (*_wake)[a] && a != _target) {
       
   238           _active_nodes[(*_dist)[a]].push_front(a);
   231           _active_nodes[(*_dist)[a]].push_front(a);
   239           if (_highest_active < (*_dist)[a]) {
   232           if (_highest_active < (*_dist)[a]) {
   240             _highest_active = (*_dist)[a];
   233             _highest_active = (*_dist)[a];
   241           }
   234           }
   242         }
   235         }
   243         (*_excess)[a] += delta;
   236         (*_excess)[a] += delta;
   244       }
   237       }
   245 
   238 
   246       _dormant[0].push_front(_source);
       
   247       (*_source_set)[_source] = true;
       
   248       _dormant_max = 0;
       
   249       (*_wake)[_source] = false;
       
   250 
       
   251       _level_size[0] = 1;
       
   252       _level_size[1] = _node_num - 1;
       
   253 
   239 
   254       do {
   240       do {
   255 	Node n;
   241 	Node n;
   256 	while ((n = findActiveNode()) != INVALID) {
   242 	while ((n = findActiveNode()) != INVALID) {
   257 	  ResEdge e;
   243 	  for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
   258 	  while (_tolerance.positive((*_excess)[n]) && 
   244             Node a = res_graph.target(e);
   259                  (e = findAdmissibleEdge(n, res_graph, current_edge)) 
   245             if ((*_dist)[a] >= (*_dist)[n] || !(*_wake)[a]) continue;
   260                  != INVALID){
   246 	    Value delta = res_graph.rescap(e);
   261 	    Value delta;
   247 	    if (_tolerance.positive((*_excess)[n] - delta)) {
   262 	    if ((*_excess)[n] < res_graph.rescap(e)) {
   248               (*_excess)[n] -= delta;
       
   249 	    } else {
   263 	      delta = (*_excess)[n];
   250 	      delta = (*_excess)[n];
   264 	    } else {
   251               (*_excess)[n] = 0;
   265 	      delta = res_graph.rescap(e);
       
   266 	      res_graph.nextOut(current_edge[n]);
       
   267 	    }
   252 	    }
   268             if (!_tolerance.positive(delta)) continue;
       
   269 	    res_graph.augment(e, delta);
   253 	    res_graph.augment(e, delta);
   270 	    (*_excess)[res_graph.source(e)] -= delta;
   254 	    if ((*_excess)[a] == 0 && a != _target) {
   271 	    Node a = res_graph.target(e);
       
   272 	    if (!_tolerance.positive((*_excess)[a]) && a != _target) {
       
   273 	      _active_nodes[(*_dist)[a]].push_front(a);
   255 	      _active_nodes[(*_dist)[a]].push_front(a);
   274 	    }
   256 	    }
   275 	    (*_excess)[a] += delta;
   257 	    (*_excess)[a] += delta;
       
   258             if ((*_excess)[n] == 0) break;
   276 	  }
   259 	  }
   277 	  if (_tolerance.positive((*_excess)[n])) {
   260 	  if ((*_excess)[n] != 0) {
   278 	    relabel(n, res_graph, current_edge);
   261 	    relabel(n, res_graph);
   279           }
   262           }
   280 	}
   263 	}
   281 
   264 
   282 	Value current_value = cutValue(out);
   265 	Value current_value = cutValue(out);
   283  	if (_min_cut > current_value){
   266  	if (_min_cut > current_value){
   295  	}
   278  	}
   296 
   279 
   297       } while (selectNewSink(res_graph));
   280       } while (selectNewSink(res_graph));
   298     }
   281     }
   299 
   282 
   300     template <typename ResGraph, typename EdgeMap>
   283     template <typename ResGraph>
   301     void relabel(const Node& n, ResGraph& res_graph, EdgeMap& current_edge) {
   284     void relabel(const Node& n, ResGraph& res_graph) {
   302       typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
   285       typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
   303 
   286 
   304       int k = (*_dist)[n];
   287       int k = (*_dist)[n];
   305       if (_level_size[k] == 1) {
   288       if (_level_size[k] == 1) {
   306 	++_dormant_max;
   289 	++_dormant_max;
   309 	    (*_wake)[it] = false;
   292 	    (*_wake)[it] = false;
   310 	    _dormant[_dormant_max].push_front(it);
   293 	    _dormant[_dormant_max].push_front(it);
   311 	    --_level_size[(*_dist)[it]];
   294 	    --_level_size[(*_dist)[it]];
   312 	  }
   295 	  }
   313 	}
   296 	}
   314 	--_highest_active;
   297         --_highest_active;
   315       } else {	
   298       } else {	
   316         int new_dist = _node_num;
   299         int new_dist = _node_num;
   317         for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
   300         for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
   318           Node t = res_graph.target(e);
   301           Node t = res_graph.target(e);
   319           if ((*_wake)[t] && new_dist > (*_dist)[t]) {
   302           if ((*_wake)[t] && new_dist > (*_dist)[t]) {
   329 	  --_level_size[(*_dist)[n]];
   312 	  --_level_size[(*_dist)[n]];
   330 	  (*_dist)[n] = new_dist + 1;
   313 	  (*_dist)[n] = new_dist + 1;
   331 	  _highest_active = (*_dist)[n];
   314 	  _highest_active = (*_dist)[n];
   332 	  _active_nodes[_highest_active].push_front(n);
   315 	  _active_nodes[_highest_active].push_front(n);
   333 	  ++_level_size[(*_dist)[n]];
   316 	  ++_level_size[(*_dist)[n]];
   334 	  res_graph.firstOut(current_edge[n], n);
       
   335 	}
   317 	}
   336       }
   318       }
   337     }
   319     }
   338 
   320 
   339     template <typename ResGraph>
   321     template <typename ResGraph>
   370 	}
   352 	}
   371       }
   353       }
   372 
   354 
   373       if (wake_was_empty){
   355       if (wake_was_empty){
   374 	for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
   356 	for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
   375           if (_tolerance.positive((*_excess)[it])) {
   357           if ((*_excess)[it] != 0 && it != _target) {
   376 	    if ((*_wake)[it] && it != _target) {
   358             _active_nodes[(*_dist)[it]].push_front(it);
   377 	      _active_nodes[(*_dist)[it]].push_front(it);
   359             if (_highest_active < (*_dist)[it]) {
   378             }
   360               _highest_active = (*_dist)[it];		    
   379 	    if (_highest_active < (*_dist)[it]) {
       
   380 	      _highest_active = (*_dist)[it];		    
       
   381             }
   361             }
   382 	  }
   362 	  }
   383 	}
   363 	}
   384       }
   364       }
   385 
   365 
   386       for (ResOutEdgeIt e(res_graph, old_target); e!=INVALID; ++e){
   366       Node n = old_target;
   387 	if (!(*_source_set)[res_graph.target(e)]) {
   367       for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e){
   388           Value delta = res_graph.rescap(e);
   368         Node a = res_graph.target(e);
   389           if (!_tolerance.positive(delta)) continue;
   369 	if (!(*_wake)[a]) continue;
   390           res_graph.augment(e, delta);
   370         Value delta = res_graph.rescap(e);
   391           (*_excess)[res_graph.source(e)] -= delta;
   371         res_graph.augment(e, delta);
   392           Node a = res_graph.target(e);
   372         (*_excess)[n] -= delta;
   393           if (!_tolerance.positive((*_excess)[a]) && 
   373         if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) {
   394               (*_wake)[a] && a != _target) {
   374           _active_nodes[(*_dist)[a]].push_front(a);
   395             _active_nodes[(*_dist)[a]].push_front(a);
   375           if (_highest_active < (*_dist)[a]) {
   396             if (_highest_active < (*_dist)[a]) {
   376             _highest_active = (*_dist)[a];
   397               _highest_active = (*_dist)[a];
       
   398             }
       
   399           }
   377           }
   400           (*_excess)[a] += delta;
   378         }
   401 	}
   379         (*_excess)[a] += delta;
   402       }
   380       }
   403       
   381       
   404       return true;
   382       return true;
   405     }
   383     }
   406     
   384 
   407     Node findActiveNode() {
   385     Node findActiveNode() {
   408       while (_highest_active > 0 && _active_nodes[_highest_active].empty()){ 
   386       while (_highest_active > 0 && _active_nodes[_highest_active].empty()){ 
   409 	--_highest_active;
   387 	--_highest_active;
   410       }
   388       }
   411       if( _highest_active > 0) {
   389       if( _highest_active > 0) {
   412        	Node n = _active_nodes[_highest_active].front();
   390        	Node n = _active_nodes[_highest_active].front();
   413 	_active_nodes[_highest_active].pop_front();
   391 	_active_nodes[_highest_active].pop_front();
   414 	return n;
   392 	return n;
   415       } else {
       
   416 	return INVALID;
       
   417       }
       
   418     }
       
   419 
       
   420     template <typename ResGraph, typename EdgeMap>
       
   421     typename ResGraph::Edge findAdmissibleEdge(const Node& n, 
       
   422                                                ResGraph& res_graph, 
       
   423                                                EdgeMap& current_edge) {
       
   424       typedef typename ResGraph::Edge ResEdge;
       
   425       ResEdge e = current_edge[n];
       
   426       while (e != INVALID && 
       
   427              ((*_dist)[n] <= (*_dist)[res_graph.target(e)] || 
       
   428               !(*_wake)[res_graph.target(e)])) {
       
   429 	res_graph.nextOut(e);
       
   430       }
       
   431       if (e != INVALID) {
       
   432 	current_edge[n] = e;	
       
   433 	return e;
       
   434       } else {
   393       } else {
   435 	return INVALID;
   394 	return INVALID;
   436       }
   395       }
   437     }
   396     }
   438 
   397 
   510         _source_set = new SourceSetMap(*_graph);
   469         _source_set = new SourceSetMap(*_graph);
   511       }
   470       }
   512       if (!_out_res_graph) {
   471       if (!_out_res_graph) {
   513         _out_res_graph = new OutResGraph(*_graph, *_capacity, *_preflow);
   472         _out_res_graph = new OutResGraph(*_graph, *_capacity, *_preflow);
   514       }
   473       }
   515       if (!_out_current_edge) {
       
   516         _out_current_edge = new OutCurrentEdgeMap(*_graph);
       
   517       }
       
   518       if (!_rev_graph) {
   474       if (!_rev_graph) {
   519         _rev_graph = new RevGraph(*_graph);
   475         _rev_graph = new RevGraph(*_graph);
   520       }
   476       }
   521       if (!_in_res_graph) {
   477       if (!_in_res_graph) {
   522         _in_res_graph = new InResGraph(*_rev_graph, *_capacity, *_preflow);
   478         _in_res_graph = new InResGraph(*_rev_graph, *_capacity, *_preflow);
   523       }
       
   524       if (!_in_current_edge) {
       
   525         _in_current_edge = new InCurrentEdgeMap(*_graph);
       
   526       }
   479       }
   527       if (!_min_cut_map) {
   480       if (!_min_cut_map) {
   528         _min_cut_map = new MinCutMap(*_graph);
   481         _min_cut_map = new MinCutMap(*_graph);
   529       }
   482       }
   530 
   483 
   553     /// \brief Calculates a minimum cut with \f$ source \f$ on the
   506     /// \brief Calculates a minimum cut with \f$ source \f$ on the
   554     /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
   507     /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
   555     ///  and minimal out-degree). The \c target is the initial target
   508     ///  and minimal out-degree). The \c target is the initial target
   556     /// for the push-relabel algorithm.
   509     /// for the push-relabel algorithm.
   557     void calculateOut(const Node& target) {
   510     void calculateOut(const Node& target) {
   558       findMinCut(target, true, *_out_res_graph, *_out_current_edge);
   511       findMinCut(target, true, *_out_res_graph);
   559     }
   512     }
   560 
   513 
   561     /// \brief Calculates a minimum cut with \f$ source \f$ on the
   514     /// \brief Calculates a minimum cut with \f$ source \f$ on the
   562     /// sink-side.
   515     /// sink-side.
   563     ///
   516     ///
   581     /// sink-side (i.e. a set \f$ X\subsetneq V 
   534     /// sink-side (i.e. a set \f$ X\subsetneq V 
   582     /// \f$ with \f$ source \notin X \f$ and minimal out-degree).  
   535     /// \f$ with \f$ source \notin X \f$ and minimal out-degree).  
   583     /// The \c target is the initial
   536     /// The \c target is the initial
   584     /// target for the push-relabel algorithm.
   537     /// target for the push-relabel algorithm.
   585     void calculateIn(const Node& target) {
   538     void calculateIn(const Node& target) {
   586       findMinCut(target, false, *_in_res_graph, *_in_current_edge);
   539       findMinCut(target, false, *_in_res_graph);
   587     }
   540     }
   588 
   541 
   589     /// \brief Runs the algorithm.
   542     /// \brief Runs the algorithm.
   590     ///
   543     ///
   591     /// Runs the algorithm. It finds nodes \c source and \c target
   544     /// Runs the algorithm. It finds nodes \c source and \c target