Changeset 2225:bb3d5e6f9fcb in lemon0.x for lemon/hao_orlin.h
 Timestamp:
 09/29/06 13:36:30 (13 years ago)
 Branch:
 default
 Phase:
 public
 Convert:
 svn:c9d7d8f590d60310b91f818b3a526b0e/lemon/trunk@2965
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/hao_orlin.h
r2211 r2225 1 1 /* * C++ * 2 * lemon/hao_orlin.h  Part of LEMON, a generic C++ optimization library3 2 * 4 * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport 3 * This file is a part of LEMON, a generic C++ optimization library 4 * 5 * Copyright (C) 20032006 6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport 5 7 * (Egervary Research Group on Combinatorial Optimization, EGRES). 6 8 * … … 30 32 /// \file 31 33 /// \ingroup flowalgs 32 /// Implementation of the HaoOrlin algorithms class for testing network 34 /// \brief Implementation of the HaoOrlin algorithm. 35 /// 36 /// Implementation of the HaoOrlin algorithms class for testing network 33 37 /// reliability. 34 38 35 39 namespace lemon { 36 40 37 /// \addtogroup flowalgs 38 /// @{ 39 40 /// %HaoOrlin algorithm for calculate minimum cut in directed graphs. 41 /// \ingroup flowalgs 42 /// 43 /// \brief %HaoOrlin algorithm for calculate minimum cut in directed graphs. 41 44 /// 42 45 /// HaoOrlin calculates the minimum cut in directed graphs. It 43 /// separates the nodes of the graph into two disjoint sets and the 44 /// summary of the edge capacities go from the first set to the 45 /// second set is the minimum. The algorithm is a modified 46 /// pushrelabel preflow algorithm and it calculates the minimum cat 47 /// in \f$ O(n^3) \f$ time. The purpose of such algorithm is testing 48 /// network reliability. For sparse undirected graph you can use the 49 /// algorithm of Nagamochi and Ibraki which solves the undirected 50 /// problem in \f$ O(n^3) \f$ time. 46 /// separates the nodes of the graph into two disjoint sets, 47 /// \f$ V_{out} \f$ and \f$ V_{in} \f$. This separation is the minimum 48 /// cut if the summary of the edge capacities which source is in 49 /// \f$ V_{out} \f$ and the target is in \f$ V_{in} \f$ is the 50 /// minimum. The algorithm is a modified pushrelabel preflow 51 /// algorithm and it calculates the minimum cut in \f$ O(n^3) \f$ 52 /// time. The purpose of such algorithm is testing network 53 /// reliability. For sparse undirected graph you can use the 54 /// algorithm of Nagamochi and Ibaraki which solves the undirected 55 /// problem in \f$ O(ne + n^2 \log(n)) \f$ time and it is implemented in the 56 /// MinCut algorithm class. 57 /// 58 /// \param _Graph is the graph type of the algorithm. 59 /// \param _CapacityMap is an edge map of capacities which should 60 /// be any numreric type. The default type is _Graph::EdgeMap<int>. 61 /// \param _Tolerance is the handler of the inexact computation. The 62 /// default type for it is Tolerance<typename CapacityMap::Value>. 51 63 /// 52 64 /// \author Attila Bernath and Balazs Dezso 65 #ifdef DOXYGEN 66 template <typename _Graph, typename _CapacityMap, typename _Tolerance> 67 #else 53 68 template <typename _Graph, 54 69 typename _CapacityMap = typename _Graph::template EdgeMap<int>, 55 70 typename _Tolerance = Tolerance<typename _CapacityMap::Value> > 71 #endif 56 72 class HaoOrlin { 57 73 protected: … … 71 87 72 88 const Graph* _graph; 89 73 90 const CapacityMap* _capacity; 74 91 … … 81 98 82 99 typedef ResGraphAdaptor<const Graph, Value, CapacityMap, 83 FlowMap, Tolerance> ResGraph; 84 typedef typename ResGraph::Node ResNode; 85 typedef typename ResGraph::NodeIt ResNodeIt; 86 typedef typename ResGraph::EdgeIt ResEdgeIt; 87 typedef typename ResGraph::OutEdgeIt ResOutEdgeIt; 88 typedef typename ResGraph::Edge ResEdge; 89 typedef typename ResGraph::InEdgeIt ResInEdgeIt; 90 91 ResGraph* _res_graph; 92 93 typedef typename Graph::template NodeMap<ResEdge> CurrentArcMap; 94 CurrentArcMap* _current_arc; 100 FlowMap, Tolerance> OutResGraph; 101 typedef typename OutResGraph::Edge OutResEdge; 102 103 OutResGraph* _out_res_graph; 104 105 typedef typename Graph::template NodeMap<OutResEdge> OutCurrentEdgeMap; 106 OutCurrentEdgeMap* _out_current_edge; 107 108 typedef RevGraphAdaptor<const Graph> RevGraph; 109 RevGraph* _rev_graph; 110 111 typedef ResGraphAdaptor<const RevGraph, Value, CapacityMap, 112 FlowMap, Tolerance> InResGraph; 113 typedef typename InResGraph::Edge InResEdge; 114 115 InResGraph* _in_res_graph; 116 117 typedef typename Graph::template NodeMap<InResEdge> InCurrentEdgeMap; 118 InCurrentEdgeMap* _in_current_edge; 95 119 96 120 … … 125 149 public: 126 150 151 /// \brief Constructor 152 /// 153 /// Constructor of the algorithm class. 127 154 HaoOrlin(const Graph& graph, const CapacityMap& capacity, 128 155 const Tolerance& tolerance = Tolerance()) : 129 156 _graph(&graph), _capacity(&capacity), 130 _preflow(0), _source(), _target(), _res_graph(0), _current_arc(0), 157 _preflow(0), _source(), _target(), 158 _out_res_graph(0), _out_current_edge(0), 159 _rev_graph(0), _in_res_graph(0), _in_current_edge(0), 131 160 _wake(0),_dist(0), _excess(0), _source_set(0), 132 161 _highest_active(), _active_nodes(), _dormant_max(), _dormant(), … … 134 163 135 164 ~HaoOrlin() { 136 if (_res_graph) {137 delete _res_graph;138 }139 165 if (_min_cut_map) { 140 166 delete _min_cut_map; 141 167 } 142 if (_current_arc) { 143 delete _current_arc; 168 if (_in_current_edge) { 169 delete _in_current_edge; 170 } 171 if (_in_res_graph) { 172 delete _in_res_graph; 173 } 174 if (_rev_graph) { 175 delete _rev_graph; 176 } 177 if (_out_current_edge) { 178 delete _out_current_edge; 179 } 180 if (_out_res_graph) { 181 delete _out_res_graph; 144 182 } 145 183 if (_source_set) { … … 162 200 private: 163 201 164 void relabel(Node i) { 165 int k = (*_dist)[i]; 202 template <typename ResGraph, typename EdgeMap> 203 void findMinCut(const Node& target, bool out, 204 ResGraph& res_graph, EdgeMap& current_edge) { 205 typedef typename ResGraph::Edge ResEdge; 206 typedef typename ResGraph::OutEdgeIt ResOutEdgeIt; 207 208 for (typename Graph::EdgeIt it(*_graph); it != INVALID; ++it) { 209 (*_preflow)[it] = 0; 210 } 211 for (NodeIt it(*_graph); it != INVALID; ++it) { 212 (*_wake)[it] = true; 213 (*_dist)[it] = 1; 214 (*_excess)[it] = 0; 215 (*_source_set)[it] = false; 216 217 res_graph.firstOut(current_edge[it], it); 218 } 219 220 _target = target; 221 (*_dist)[target] = 0; 222 223 for (ResOutEdgeIt it(res_graph, _source); it != INVALID; ++it) { 224 Value delta = res_graph.rescap(it); 225 if (!_tolerance.positive(delta)) continue; 226 227 (*_excess)[res_graph.source(it)] = delta; 228 res_graph.augment(it, delta); 229 Node a = res_graph.target(it); 230 if (!_tolerance.positive((*_excess)[a]) && 231 (*_wake)[a] && a != _target) { 232 _active_nodes[(*_dist)[a]].push_front(a); 233 if (_highest_active < (*_dist)[a]) { 234 _highest_active = (*_dist)[a]; 235 } 236 } 237 (*_excess)[a] += delta; 238 } 239 240 _dormant[0].push_front(_source); 241 (*_source_set)[_source] = true; 242 _dormant_max = 0; 243 (*_wake)[_source] = false; 244 245 _level_size[0] = 1; 246 _level_size[1] = _node_num  1; 247 248 do { 249 Node n; 250 while ((n = findActiveNode()) != INVALID) { 251 ResEdge e; 252 while (_tolerance.positive((*_excess)[n]) && 253 (e = findAdmissibleEdge(n, res_graph, current_edge)) 254 != INVALID){ 255 Value delta; 256 if ((*_excess)[n] < res_graph.rescap(e)) { 257 delta = (*_excess)[n]; 258 } else { 259 delta = res_graph.rescap(e); 260 res_graph.nextOut(current_edge[n]); 261 } 262 if (!_tolerance.positive(delta)) continue; 263 res_graph.augment(e, delta); 264 (*_excess)[res_graph.source(e)] = delta; 265 Node a = res_graph.target(e); 266 if (!_tolerance.positive((*_excess)[a]) && a != _target) { 267 _active_nodes[(*_dist)[a]].push_front(a); 268 } 269 (*_excess)[a] += delta; 270 } 271 if (_tolerance.positive((*_excess)[n])) { 272 relabel(n, res_graph, current_edge); 273 } 274 } 275 276 Value current_value = cutValue(out); 277 if (_min_cut > current_value){ 278 if (out) { 279 for (NodeIt it(*_graph); it != INVALID; ++it) { 280 _min_cut_map>set(it, !(*_wake)[it]); 281 } 282 } else { 283 for (NodeIt it(*_graph); it != INVALID; ++it) { 284 _min_cut_map>set(it, (*_wake)[it]); 285 } 286 } 287 288 _min_cut = current_value; 289 } 290 291 } while (selectNewSink(res_graph)); 292 } 293 294 template <typename ResGraph, typename EdgeMap> 295 void relabel(const Node& n, ResGraph& res_graph, EdgeMap& current_edge) { 296 typedef typename ResGraph::OutEdgeIt ResOutEdgeIt; 297 298 int k = (*_dist)[n]; 166 299 if (_level_size[k] == 1) { 167 300 ++_dormant_max; … … 174 307 } 175 308 _highest_active; 176 } else { 177 ResOutEdgeIt e(*_res_graph, i); 178 while (e != INVALID && !(*_wake)[_res_graph>target(e)]) { 179 ++e; 309 } else { 310 int new_dist = _node_num; 311 for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) { 312 Node t = res_graph.target(e); 313 if ((*_wake)[t] && new_dist > (*_dist)[t]) { 314 new_dist = (*_dist)[t]; 315 } 316 } 317 if (new_dist == _node_num) { 318 ++_dormant_max; 319 (*_wake)[n] = false; 320 _dormant[_dormant_max].push_front(n); 321 _level_size[(*_dist)[n]]; 322 } else { 323 _level_size[(*_dist)[n]]; 324 (*_dist)[n] = new_dist + 1; 325 _highest_active = (*_dist)[n]; 326 _active_nodes[_highest_active].push_front(n); 327 ++_level_size[(*_dist)[n]]; 328 res_graph.firstOut(current_edge[n], n); 180 329 } 181 182 if (e == INVALID){ 183 ++_dormant_max; 184 (*_wake)[i] = false; 185 _dormant[_dormant_max].push_front(i); 186 _level_size[(*_dist)[i]]; 187 } else{ 188 Node j = _res_graph>target(e); 189 int new_dist = (*_dist)[j]; 190 ++e; 191 while (e != INVALID){ 192 Node j = _res_graph>target(e); 193 if ((*_wake)[j] && new_dist > (*_dist)[j]) { 194 new_dist = (*_dist)[j]; 195 } 196 ++e; 197 } 198 _level_size[(*_dist)[i]]; 199 (*_dist)[i] = new_dist + 1; 200 _highest_active = (*_dist)[i]; 201 _active_nodes[_highest_active].push_front(i); 202 ++_level_size[(*_dist)[i]]; 203 _res_graph>firstOut((*_current_arc)[i], i); 204 } 205 } 206 } 207 208 bool selectNewSink(){ 330 } 331 } 332 333 template <typename ResGraph> 334 bool selectNewSink(ResGraph& res_graph) { 335 typedef typename ResGraph::OutEdgeIt ResOutEdgeIt; 336 209 337 Node old_target = _target; 210 338 (*_wake)[_target] = false; … … 250 378 } 251 379 252 for (ResOutEdgeIt e(*_res_graph, old_target); e!=INVALID; ++e){ 253 if (!(*_source_set)[_res_graph>target(e)]){ 254 push(e, _res_graph>rescap(e)); 380 for (ResOutEdgeIt e(res_graph, old_target); e!=INVALID; ++e){ 381 if (!(*_source_set)[res_graph.target(e)]) { 382 Value delta = res_graph.rescap(e); 383 if (!_tolerance.positive(delta)) continue; 384 res_graph.augment(e, delta); 385 (*_excess)[res_graph.source(e)] = delta; 386 Node a = res_graph.target(e); 387 if (!_tolerance.positive((*_excess)[a]) && 388 (*_wake)[a] && a != _target) { 389 _active_nodes[(*_dist)[a]].push_front(a); 390 if (_highest_active < (*_dist)[a]) { 391 _highest_active = (*_dist)[a]; 392 } 393 } 394 (*_excess)[a] += delta; 255 395 } 256 396 } … … 272 412 } 273 413 274 ResEdge findAdmissibleEdge(const Node& n){ 275 ResEdge e = (*_current_arc)[n]; 414 template <typename ResGraph, typename EdgeMap> 415 typename ResGraph::Edge findAdmissibleEdge(const Node& n, 416 ResGraph& res_graph, 417 EdgeMap& current_edge) { 418 typedef typename ResGraph::Edge ResEdge; 419 ResEdge e = current_edge[n]; 276 420 while (e != INVALID && 277 ((*_dist)[n] <= (*_dist)[ _res_graph>target(e)] 278 !(*_wake)[ _res_graph>target(e)])) {279 _res_graph>nextOut(e);421 ((*_dist)[n] <= (*_dist)[res_graph.target(e)]  422 !(*_wake)[res_graph.target(e)])) { 423 res_graph.nextOut(e); 280 424 } 281 425 if (e != INVALID) { 282 (*_current_arc)[n] = e;426 current_edge[n] = e; 283 427 return e; 284 428 } else { … … 287 431 } 288 432 289 void push(ResEdge& e,const Value& delta){ 290 _res_graph>augment(e, delta); 291 if (!_tolerance.positive(delta)) return; 292 293 (*_excess)[_res_graph>source(e)] = delta; 294 Node a = _res_graph>target(e); 295 if (!_tolerance.positive((*_excess)[a]) && (*_wake)[a] && a != _target) { 296 _active_nodes[(*_dist)[a]].push_front(a); 297 if (_highest_active < (*_dist)[a]) { 298 _highest_active = (*_dist)[a]; 299 } 300 } 301 (*_excess)[a] += delta; 302 } 303 304 Value cutValue() { 433 Value cutValue(bool out) { 305 434 Value value = 0; 306 for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) { 307 for (InEdgeIt e(*_graph, it); e != INVALID; ++e) { 308 if (!(*_wake)[_graph>source(e)]){ 309 value += (*_capacity)[e]; 310 } 311 } 435 if (out) { 436 for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) { 437 for (InEdgeIt e(*_graph, it); e != INVALID; ++e) { 438 if (!(*_wake)[_graph>source(e)]){ 439 value += (*_capacity)[e]; 440 } 441 } 442 } 443 } else { 444 for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) { 445 for (OutEdgeIt e(*_graph, it); e != INVALID; ++e) { 446 if (!(*_wake)[_graph>target(e)]){ 447 value += (*_capacity)[e]; 448 } 449 } 450 } 312 451 } 313 452 return value; 314 } 453 } 454 315 455 316 456 public: 317 457 458 /// \name Execution control 459 /// The simplest way to execute the algorithm is to use 460 /// one of the member functions called \c run(...). 461 /// \n 462 /// If you need more control on the execution, 463 /// first you must call \ref init(), then the \ref calculateIn() or 464 /// \ref calculateIn() functions. 465 466 /// @{ 467 318 468 /// \brief Initializes the internal data structures. 319 469 /// 320 470 /// Initializes the internal data structures. It creates 321 /// the maps, residual graph adaptor and some bucket structures471 /// the maps, residual graph adaptors and some bucket structures 322 472 /// for the algorithm. 323 473 void init() { … … 354 504 _source_set = new SourceSetMap(*_graph); 355 505 } 356 if (!_current_arc) { 357 _current_arc = new CurrentArcMap(*_graph); 506 if (!_out_res_graph) { 507 _out_res_graph = new OutResGraph(*_graph, *_capacity, *_preflow); 508 } 509 if (!_out_current_edge) { 510 _out_current_edge = new OutCurrentEdgeMap(*_graph); 511 } 512 if (!_rev_graph) { 513 _rev_graph = new RevGraph(*_graph); 514 } 515 if (!_in_res_graph) { 516 _in_res_graph = new InResGraph(*_rev_graph, *_capacity, *_preflow); 517 } 518 if (!_in_current_edge) { 519 _in_current_edge = new InCurrentEdgeMap(*_graph); 358 520 } 359 521 if (!_min_cut_map) { 360 522 _min_cut_map = new MinCutMap(*_graph); 361 523 } 362 if (!_res_graph) {363 _res_graph = new ResGraph(*_graph, *_capacity, *_preflow);364 }365 524 366 525 _min_cut = std::numeric_limits<Value>::max(); … … 369 528 370 529 /// \brief Calculates the minimum cut with the \c source node 371 /// in the first partition.530 /// in the \f$ V_{out} \f$. 372 531 /// 373 532 /// Calculates the minimum cut with the \c source node 374 /// in the first partition.533 /// in the \f$ V_{out} \f$. 375 534 void calculateOut() { 376 535 for (NodeIt it(*_graph); it != INVALID; ++it) { … … 383 542 384 543 /// \brief Calculates the minimum cut with the \c source node 385 /// in the first partition.544 /// in the \f$ V_{out} \f$. 386 545 /// 387 546 /// Calculates the minimum cut with the \c source node 388 /// in the first partition. The \c target is the initial target547 /// in the \f$ V_{out} \f$. The \c target is the initial target 389 548 /// for the pushrelabel algorithm. 390 549 void calculateOut(const Node& target) { 391 for (NodeIt it(*_graph); it != INVALID; ++it) { 392 (*_wake)[it] = true; 393 (*_dist)[it] = 1; 394 (*_excess)[it] = 0; 395 (*_source_set)[it] = false; 396 397 _res_graph>firstOut((*_current_arc)[it], it); 398 } 399 400 _target = target; 401 (*_dist)[target] = 0; 402 403 for (ResOutEdgeIt it(*_res_graph, _source); it != INVALID; ++it) { 404 push(it, _res_graph>rescap(it)); 405 } 406 407 _dormant[0].push_front(_source); 408 (*_source_set)[_source] = true; 409 _dormant_max = 0; 410 (*_wake)[_source]=false; 411 412 _level_size[0] = 1; 413 _level_size[1] = _node_num  1; 414 415 do { 416 Node n; 417 while ((n = findActiveNode()) != INVALID) { 418 ResEdge e; 419 while (_tolerance.positive((*_excess)[n]) && 420 (e = findAdmissibleEdge(n)) != INVALID){ 421 Value delta; 422 if ((*_excess)[n] < _res_graph>rescap(e)) { 423 delta = (*_excess)[n]; 424 } else { 425 delta = _res_graph>rescap(e); 426 _res_graph>nextOut((*_current_arc)[n]); 427 } 428 if (!_tolerance.positive(delta)) continue; 429 _res_graph>augment(e, delta); 430 (*_excess)[_res_graph>source(e)] = delta; 431 Node a = _res_graph>target(e); 432 if (!_tolerance.positive((*_excess)[a]) && a != _target) { 433 _active_nodes[(*_dist)[a]].push_front(a); 434 } 435 (*_excess)[a] += delta; 436 } 437 if (_tolerance.positive((*_excess)[n])) { 438 relabel(n); 439 } 440 } 441 442 Value current_value = cutValue(); 443 if (_min_cut > current_value){ 444 for (NodeIt it(*_graph); it != INVALID; ++it) { 445 _min_cut_map>set(it, !(*_wake)[it]); 446 } 447 448 _min_cut = current_value; 449 } 450 451 } while (selectNewSink()); 452 } 453 550 findMinCut(target, true, *_out_res_graph, *_out_current_edge); 551 } 552 553 /// \brief Calculates the minimum cut with the \c source node 554 /// in the \f$ V_{in} \f$. 555 /// 556 /// Calculates the minimum cut with the \c source node 557 /// in the \f$ V_{in} \f$. 454 558 void calculateIn() { 455 559 for (NodeIt it(*_graph); it != INVALID; ++it) { … … 461 565 } 462 566 567 /// \brief Calculates the minimum cut with the \c source node 568 /// in the \f$ V_{in} \f$. 569 /// 570 /// Calculates the minimum cut with the \c source node 571 /// in the \f$ V_{in} \f$. The \c target is the initial target 572 /// for the pushrelabel algorithm. 573 void calculateIn(const Node& target) { 574 findMinCut(target, false, *_in_res_graph, *_in_current_edge); 575 } 576 577 /// \brief Runs the algorithm. 578 /// 579 /// Runs the algorithm. It finds a proper \c source and \c target 580 /// and then calls the \ref init(), \ref calculateOut() and \ref 581 /// calculateIn(). 463 582 void run() { 464 583 init(); 465 584 for (NodeIt it(*_graph); it != INVALID; ++it) { 466 585 if (it != _source) { 467 startOut(it);468 // startIn(it);586 calculateOut(it); 587 calculateIn(it); 469 588 return; 470 589 } … … 472 591 } 473 592 593 /// \brief Runs the algorithm. 594 /// 595 /// Runs the algorithm. It finds a proper \c target and then calls 596 /// the \ref init(), \ref calculateOut() and \ref calculateIn(). 474 597 void run(const Node& s) { 475 598 init(s); 476 599 for (NodeIt it(*_graph); it != INVALID; ++it) { 477 600 if (it != _source) { 478 startOut(it);479 // startIn(it);601 calculateOut(it); 602 calculateIn(it); 480 603 return; 481 604 } … … 483 606 } 484 607 608 /// \brief Runs the algorithm. 609 /// 610 /// Runs the algorithm. It just calls the \ref init() and then 611 /// \ref calculateOut() and \ref calculateIn(). 485 612 void run(const Node& s, const Node& t) { 486 init(s); 487 startOut(t); 488 startIn(t); 489 } 490 491 /// \brief Returns the value of the minimum value cut with node \c 492 /// source on the source side. 613 init(s); 614 calculateOut(t); 615 calculateIn(t); 616 } 617 618 /// @} 619 620 /// \name Query Functions The result of the %HaoOrlin algorithm 621 /// can be obtained using these functions. 622 /// \n 623 /// Before the use of these functions, either \ref run(), \ref 624 /// calculateOut() or \ref calculateIn() must be called. 625 626 /// @{ 627 628 /// \brief Returns the value of the minimum value cut. 493 629 /// 494 /// Returns the value of the minimum value cut with node \c source 495 /// on the source side. This function can be called after running 496 /// \ref findMinCut. 630 /// Returns the value of the minimum value cut. 497 631 Value minCut() const { 498 632 return _min_cut; … … 503 637 /// 504 638 /// Sets \c nodeMap to the characteristic vector of a minimum 505 /// value cut with node \c source on the source side. This506 /// function can be called after running \ref findMinCut.639 /// value cut. The nodes in \f$ V_{out} \f$ will be set true and 640 /// the nodes in \f$ V_{in} \f$ will be set false. 507 641 /// \pre nodeMap should be a boolvalued nodemap. 508 642 template <typename NodeMap> … … 513 647 return minCut(); 514 648 } 649 650 /// @} 515 651 516 652 }; //class HaoOrlin
Note: See TracChangeset
for help on using the changeset viewer.