Changeset 2514:57143c09dc20 in lemon-0.x for lemon/preflow.h
- Timestamp:
- 11/17/07 21:58:11 (17 years ago)
- Branch:
- default
- Phase:
- public
- Convert:
- svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3379
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
lemon/preflow.h
r2473 r2514 20 20 #define LEMON_PREFLOW_H 21 21 22 #include <vector>23 #include <queue>24 25 22 #include <lemon/error.h> 26 23 #include <lemon/bits/invalid.h> … … 28 25 #include <lemon/maps.h> 29 26 #include <lemon/graph_utils.h> 27 #include <lemon/elevator.h> 30 28 31 29 /// \file … … 33 31 /// \brief Implementation of the preflow algorithm. 34 32 35 namespace lemon { 36 37 ///\ingroup max_flow 38 ///\brief %Preflow algorithms class. 33 namespace lemon { 34 35 /// \brief Default traits class of Preflow class. 39 36 /// 40 ///This class provides an implementation of the \e preflow \e 41 ///algorithm producing a flow of maximum value in a directed 42 ///graph. The preflow algorithms are the fastest known max flow algorithms. 43 ///The \e source node, the \e target node, the \e 44 ///capacity of the edges and the \e starting \e flow value of the 45 ///edges should be passed to the algorithm through the 46 ///constructor. It is possible to change these quantities using the 47 ///functions \ref source, \ref target, \ref capacityMap and \ref 48 ///flowMap. 37 /// Default traits class of Preflow class. 38 /// \param _Graph Graph type. 39 /// \param _CapacityMap Type of capacity map. 40 template <typename _Graph, typename _CapacityMap> 41 struct PreflowDefaultTraits { 42 43 /// \brief The graph type the algorithm runs on. 44 typedef _Graph Graph; 45 46 /// \brief The type of the map that stores the edge capacities. 47 /// 48 /// The type of the map that stores the edge capacities. 49 /// It must meet the \ref concepts::ReadMap "ReadMap" concept. 50 typedef _CapacityMap CapacityMap; 51 52 /// \brief The type of the length of the edges. 53 typedef typename CapacityMap::Value Value; 54 55 /// \brief The map type that stores the flow values. 56 /// 57 /// The map type that stores the flow values. 58 /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept. 59 typedef typename Graph::template EdgeMap<Value> FlowMap; 60 61 /// \brief Instantiates a FlowMap. 62 /// 63 /// This function instantiates a \ref FlowMap. 64 /// \param graph The graph, to which we would like to define the flow map. 65 static FlowMap* createFlowMap(const Graph& graph) { 66 return new FlowMap(graph); 67 } 68 69 /// \brief The eleavator type used by Preflow algorithm. 70 /// 71 /// The elevator type used by Preflow algorithm. 72 /// 73 /// \sa Elevator 74 /// \sa LinkedElevator 75 typedef LinkedElevator<Graph, typename Graph::Node> Elevator; 76 77 /// \brief Instantiates an Elevator. 78 /// 79 /// This function instantiates a \ref Elevator. 80 /// \param graph The graph, to which we would like to define the elevator. 81 /// \param max_level The maximum level of the elevator. 82 static Elevator* createElevator(const Graph& graph, int max_level) { 83 return new Elevator(graph, max_level); 84 } 85 86 /// \brief The tolerance used by the algorithm 87 /// 88 /// The tolerance used by the algorithm to handle inexact computation. 89 typedef Tolerance<Value> Tolerance; 90 91 }; 92 93 94 /// \ingroup max_flow 49 95 /// 50 ///After running \ref lemon::Preflow::phase1() "phase1()" 51 ///or \ref lemon::Preflow::run() "run()", the maximal flow 52 ///value can be obtained by calling \ref flowValue(). The minimum 53 ///value cut can be written into a <tt>bool</tt> node map by 54 ///calling \ref minCut(). (\ref minMinCut() and \ref maxMinCut() writes 55 ///the inclusionwise minimum and maximum of the minimum value cuts, 56 ///resp.) 96 /// \brief %Preflow algorithms class. 57 97 /// 58 ///\param Graph The directed graph type the algorithm runs on. 59 ///\param Num The number type of the capacities and the flow values. 60 ///\param CapacityMap The capacity map type. 61 ///\param FlowMap The flow map type. 62 ///\param Tol The tolerance type. 98 /// This class provides an implementation of the Goldberg's \e 99 /// preflow \e algorithm producing a flow of maximum value in a 100 /// directed graph. The preflow algorithms are the fastest known max 101 /// flow algorithms. The current implementation use a mixture of the 102 /// \e "highest label" and the \e "bound decrease" heuristics. 103 /// The worst case time complexity of the algorithm is \f$O(n^3)\f$. 63 104 /// 64 ///\author Jacint Szabo 65 ///\todo Second template parameter is superfluous 66 template <typename Graph, typename Num, 67 typename CapacityMap=typename Graph::template EdgeMap<Num>, 68 typename FlowMap=typename Graph::template EdgeMap<Num>, 69 typename Tol=Tolerance<Num> > 105 /// The algorithm consists from two phases. After the first phase 106 /// the maximal flow value and the minimum cut can be obtained. The 107 /// second phase constructs the feasible maximum flow on each edge. 108 /// 109 /// \param _Graph The directed graph type the algorithm runs on. 110 /// \param _CapacityMap The flow map type. 111 /// \param _Traits Traits class to set various data types used by 112 /// the algorithm. The default traits class is \ref 113 /// PreflowDefaultTraits. See \ref PreflowDefaultTraits for the 114 /// documentation of a %Preflow traits class. 115 /// 116 ///\author Jacint Szabo and Balazs Dezso 117 #ifdef DOXYGEN 118 template <typename _Graph, typename _CapacityMap, typename _Traits> 119 #else 120 template <typename _Graph, 121 typename _CapacityMap = typename _Graph::template EdgeMap<int>, 122 typename _Traits = PreflowDefaultTraits<_Graph, _CapacityMap> > 123 #endif 70 124 class Preflow { 71 protected: 72 typedef typename Graph::Node Node; 73 typedef typename Graph::NodeIt NodeIt; 74 typedef typename Graph::EdgeIt EdgeIt; 75 typedef typename Graph::OutEdgeIt OutEdgeIt; 76 typedef typename Graph::InEdgeIt InEdgeIt; 77 78 typedef typename Graph::template NodeMap<Node> NNMap; 79 typedef typename std::vector<Node> VecNode; 80 81 const Graph* _g; 82 Node _source; 83 Node _target; 84 const CapacityMap* _capacity; 85 FlowMap* _flow; 86 87 Tol _surely; 125 public: 88 126 89 int _node_num; //the number of nodes of G 90 91 typename Graph::template NodeMap<int> level; 92 typename Graph::template NodeMap<Num> excess; 93 94 // constants used for heuristics 95 static const int H0=20; 96 static const int H1=1; 97 98 public: 99 100 ///\ref Exception for the case when s=t. 101 102 ///\ref Exception for the case when the source equals the target. 103 class InvalidArgument : public lemon::LogicError { 127 typedef _Traits Traits; 128 typedef typename Traits::Graph Graph; 129 typedef typename Traits::CapacityMap CapacityMap; 130 typedef typename Traits::Value Value; 131 132 typedef typename Traits::FlowMap FlowMap; 133 typedef typename Traits::Elevator Elevator; 134 typedef typename Traits::Tolerance Tolerance; 135 136 /// \brief \ref Exception for uninitialized parameters. 137 /// 138 /// This error represents problems in the initialization 139 /// of the parameters of the algorithms. 140 class UninitializedParameter : public lemon::UninitializedParameter { 104 141 public: 105 142 virtual const char* what() const throw() { 106 return "lemon::Preflow:: InvalidArgument";143 return "lemon::Preflow::UninitializedParameter"; 107 144 } 108 145 }; 146 147 private: 148 149 GRAPH_TYPEDEFS(typename Graph); 150 151 const Graph& _graph; 152 const CapacityMap* _capacity; 153 154 int _node_num; 155 156 Node _source, _target; 157 158 FlowMap* _flow; 159 bool _local_flow; 160 161 Elevator* _level; 162 bool _local_level; 163 164 typedef typename Graph::template NodeMap<Value> ExcessMap; 165 ExcessMap* _excess; 166 167 Tolerance _tolerance; 168 169 bool _phase; 170 171 void createStructures() { 172 _node_num = countNodes(_graph); 173 174 if (!_flow) { 175 _flow = Traits::createFlowMap(_graph); 176 _local_flow = true; 177 } 178 if (!_level) { 179 _level = Traits::createElevator(_graph, _node_num); 180 _local_level = true; 181 } 182 if (!_excess) { 183 _excess = new ExcessMap(_graph); 184 } 185 } 186 187 void destroyStructures() { 188 if (_local_flow) { 189 delete _flow; 190 } 191 if (_local_level) { 192 delete _level; 193 } 194 if (_excess) { 195 delete _excess; 196 } 197 } 198 199 public: 200 201 typedef Preflow Create; 202 203 ///\name Named template parameters 204 205 ///@{ 206 207 template <typename _FlowMap> 208 struct DefFlowMapTraits : public Traits { 209 typedef _FlowMap FlowMap; 210 static FlowMap *createFlowMap(const Graph&) { 211 throw UninitializedParameter(); 212 } 213 }; 214 215 /// \brief \ref named-templ-param "Named parameter" for setting 216 /// FlowMap type 217 /// 218 /// \ref named-templ-param "Named parameter" for setting FlowMap 219 /// type 220 template <typename _FlowMap> 221 struct DefFlowMap 222 : public Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > { 223 typedef Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > Create; 224 }; 225 226 template <typename _Elevator> 227 struct DefElevatorTraits : public Traits { 228 typedef _Elevator Elevator; 229 static Elevator *createElevator(const Graph&, int) { 230 throw UninitializedParameter(); 231 } 232 }; 233 234 /// \brief \ref named-templ-param "Named parameter" for setting 235 /// Elevator type 236 /// 237 /// \ref named-templ-param "Named parameter" for setting Elevator 238 /// type 239 template <typename _Elevator> 240 struct DefElevator 241 : public Preflow<Graph, CapacityMap, DefElevatorTraits<_Elevator> > { 242 typedef Preflow<Graph, CapacityMap, DefElevatorTraits<_Elevator> > Create; 243 }; 244 245 template <typename _Elevator> 246 struct DefStandardElevatorTraits : public Traits { 247 typedef _Elevator Elevator; 248 static Elevator *createElevator(const Graph& graph, int max_level) { 249 return new Elevator(graph, max_level); 250 } 251 }; 252 253 /// \brief \ref named-templ-param "Named parameter" for setting 254 /// Elevator type 255 /// 256 /// \ref named-templ-param "Named parameter" for setting Elevator 257 /// type. The Elevator should be standard constructor interface, ie. 258 /// the graph and the maximum level should be passed to it. 259 template <typename _Elevator> 260 struct DefStandardElevator 261 : public Preflow<Graph, CapacityMap, 262 DefStandardElevatorTraits<_Elevator> > { 263 typedef Preflow<Graph, CapacityMap, 264 DefStandardElevatorTraits<_Elevator> > Create; 265 }; 266 267 /// @} 268 269 /// \brief The constructor of the class. 270 /// 271 /// The constructor of the class. 272 /// \param graph The directed graph the algorithm runs on. 273 /// \param capacity The capacity of the edges. 274 /// \param source The source node. 275 /// \param target The target node. 276 Preflow(const Graph& graph, const CapacityMap& capacity, 277 Node source, Node target) 278 : _graph(graph), _capacity(&capacity), 279 _node_num(0), _source(source), _target(target), 280 _flow(0), _local_flow(false), 281 _level(0), _local_level(false), 282 _excess(0), _tolerance(), _phase() {} 283 284 /// \brief Destrcutor. 285 /// 286 /// Destructor. 287 ~Preflow() { 288 destroyStructures(); 289 } 290 291 /// \brief Sets the capacity map. 292 /// 293 /// Sets the capacity map. 294 /// \return \c (*this) 295 Preflow& capacityMap(const CapacityMap& map) { 296 _capacity = ↦ 297 return *this; 298 } 299 300 /// \brief Sets the flow map. 301 /// 302 /// Sets the flow map. 303 /// \return \c (*this) 304 Preflow& flowMap(FlowMap& map) { 305 if (_local_flow) { 306 delete _flow; 307 _local_flow = false; 308 } 309 _flow = ↦ 310 return *this; 311 } 312 313 /// \brief Returns the flow map. 314 /// 315 /// \return The flow map. 316 const FlowMap& flowMap() { 317 return *_flow; 318 } 319 320 /// \brief Sets the elevator. 321 /// 322 /// Sets the elevator. 323 /// \return \c (*this) 324 Preflow& elevator(Elevator& elevator) { 325 if (_local_level) { 326 delete _level; 327 _local_level = false; 328 } 329 _level = &elevator; 330 return *this; 331 } 332 333 /// \brief Returns the elevator. 334 /// 335 /// \return The elevator. 336 const Elevator& elevator() { 337 return *_level; 338 } 339 340 /// \brief Sets the source node. 341 /// 342 /// Sets the source node. 343 /// \return \c (*this) 344 Preflow& source(const Node& node) { 345 _source = node; 346 return *this; 347 } 348 349 /// \brief Sets the target node. 350 /// 351 /// Sets the target node. 352 /// \return \c (*this) 353 Preflow& target(const Node& node) { 354 _target = node; 355 return *this; 356 } 357 358 /// \brief Sets the tolerance used by algorithm. 359 /// 360 /// Sets the tolerance used by algorithm. 361 Preflow& tolerance(const Tolerance& tolerance) const { 362 _tolerance = tolerance; 363 return *this; 364 } 365 366 /// \brief Returns the tolerance used by algorithm. 367 /// 368 /// Returns the tolerance used by algorithm. 369 const Tolerance& tolerance() const { 370 return tolerance; 371 } 372 373 /// \name Execution control The simplest way to execute the 374 /// algorithm is to use one of the member functions called \c 375 /// run(). 376 /// \n 377 /// If you need more control on initial solution or 378 /// execution then you have to call one \ref init() function and then 379 /// the startFirstPhase() and if you need the startSecondPhase(). 109 380 381 ///@{ 382 383 /// \brief Initializes the internal data structures. 384 /// 385 /// Initializes the internal data structures. 386 /// 387 void init() { 388 createStructures(); 389 390 _phase = true; 391 for (NodeIt n(_graph); n != INVALID; ++n) { 392 _excess->set(n, 0); 393 } 394 395 for (EdgeIt e(_graph); e != INVALID; ++e) { 396 _flow->set(e, 0); 397 } 398 399 typename Graph::template NodeMap<bool> reached(_graph, false); 400 401 _level->initStart(); 402 _level->initAddItem(_target); 403 404 std::vector<Node> queue; 405 reached.set(_source, true); 406 407 queue.push_back(_target); 408 reached.set(_target, true); 409 while (!queue.empty()) { 410 std::vector<Node> nqueue; 411 for (int i = 0; i < int(queue.size()); ++i) { 412 Node n = queue[i]; 413 for (InEdgeIt e(_graph, n); e != INVALID; ++e) { 414 Node u = _graph.source(e); 415 if (!reached[u] && _tolerance.positive((*_capacity)[e])) { 416 reached.set(u, true); 417 _level->initAddItem(u); 418 nqueue.push_back(u); 419 } 420 } 421 } 422 queue.swap(nqueue); 423 } 424 _level->initFinish(); 425 426 for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) { 427 if (_tolerance.positive((*_capacity)[e])) { 428 Node u = _graph.target(e); 429 if ((*_level)[u] == _level->maxLevel()) continue; 430 _flow->set(e, (*_capacity)[e]); 431 _excess->set(u, (*_excess)[u] + (*_capacity)[e]); 432 if (u != _target && !_level->active(u)) { 433 _level->activate(u); 434 } 435 } 436 } 437 } 438 439 /// \brief Initializes the internal data structures. 440 /// 441 /// Initializes the internal data structures and sets the initial 442 /// flow to the given \c flowMap. The \c flowMap should contain a 443 /// flow or at least a preflow, ie. in each node excluding the 444 /// target the incoming flow should greater or equal to the 445 /// outgoing flow. 446 /// \return %False when the given \c flowMap is not a preflow. 447 template <typename FlowMap> 448 bool flowInit(const FlowMap& flowMap) { 449 createStructures(); 450 451 for (EdgeIt e(_graph); e != INVALID; ++e) { 452 _flow->set(e, flowMap[e]); 453 } 454 455 for (NodeIt n(_graph); n != INVALID; ++n) { 456 Value excess = 0; 457 for (InEdgeIt e(_graph, n); e != INVALID; ++e) { 458 excess += (*_flow)[e]; 459 } 460 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) { 461 excess -= (*_flow)[e]; 462 } 463 if (excess < 0 && n != _source) return false; 464 _excess->set(n, excess); 465 } 466 467 typename Graph::template NodeMap<bool> reached(_graph, false); 468 469 _level->initStart(); 470 _level->initAddItem(_target); 471 472 std::vector<Node> queue; 473 reached.set(_source, true); 474 475 queue.push_back(_target); 476 reached.set(_target, true); 477 while (!queue.empty()) { 478 _level->initNewLevel(); 479 std::vector<Node> nqueue; 480 for (int i = 0; i < int(queue.size()); ++i) { 481 Node n = queue[i]; 482 for (InEdgeIt e(_graph, n); e != INVALID; ++e) { 483 Node u = _graph.source(e); 484 if (!reached[u] && 485 _tolerance.positive((*_capacity)[e] - (*_flow)[e])) { 486 reached.set(u, true); 487 _level->initAddItem(u); 488 nqueue.push_back(u); 489 } 490 } 491 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) { 492 Node v = _graph.target(e); 493 if (!reached[v] && _tolerance.positive((*_flow)[e])) { 494 reached.set(v, true); 495 _level->initAddItem(v); 496 nqueue.push_back(v); 497 } 498 } 499 } 500 queue.swap(nqueue); 501 } 502 _level->initFinish(); 503 504 for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) { 505 Value rem = (*_capacity)[e] - (*_flow)[e]; 506 if (_tolerance.positive(rem)) { 507 Node u = _graph.target(e); 508 if ((*_level)[u] == _level->maxLevel()) continue; 509 _flow->set(e, (*_capacity)[e]); 510 _excess->set(u, (*_excess)[u] + rem); 511 if (u != _target && !_level->active(u)) { 512 _level->activate(u); 513 } 514 } 515 } 516 for (InEdgeIt e(_graph, _source); e != INVALID; ++e) { 517 Value rem = (*_flow)[e]; 518 if (_tolerance.positive(rem)) { 519 Node v = _graph.source(e); 520 if ((*_level)[v] == _level->maxLevel()) continue; 521 _flow->set(e, 0); 522 _excess->set(v, (*_excess)[v] + rem); 523 if (v != _target && !_level->active(v)) { 524 _level->activate(v); 525 } 526 } 527 } 528 return true; 529 } 530 531 /// \brief Starts the first phase of the preflow algorithm. 532 /// 533 /// The preflow algorithm consists of two phases, this method runs 534 /// the first phase. After the first phase the maximum flow value 535 /// and a minimum value cut can already be computed, although a 536 /// maximum flow is not yet obtained. So after calling this method 537 /// \ref flowValue() returns the value of a maximum flow and \ref 538 /// minCut() returns a minimum cut. 539 /// \pre One of the \ref init() functions should be called. 540 void startFirstPhase() { 541 _phase = true; 542 543 Node n = _level->highestActive(); 544 int level = _level->highestActiveLevel(); 545 while (n != INVALID) { 546 int num = _node_num; 547 548 while (num > 0 && n != INVALID) { 549 Value excess = (*_excess)[n]; 550 int new_level = _level->maxLevel(); 551 552 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) { 553 Value rem = (*_capacity)[e] - (*_flow)[e]; 554 if (!_tolerance.positive(rem)) continue; 555 Node v = _graph.target(e); 556 if ((*_level)[v] < level) { 557 if (!_level->active(v) && v != _target) { 558 _level->activate(v); 559 } 560 if (!_tolerance.less(rem, excess)) { 561 _flow->set(e, (*_flow)[e] + excess); 562 _excess->set(v, (*_excess)[v] + excess); 563 excess = 0; 564 goto no_more_push_1; 565 } else { 566 excess -= rem; 567 _excess->set(v, (*_excess)[v] + rem); 568 _flow->set(e, (*_capacity)[e]); 569 } 570 } else if (new_level > (*_level)[v]) { 571 new_level = (*_level)[v]; 572 } 573 } 574 575 for (InEdgeIt e(_graph, n); e != INVALID; ++e) { 576 Value rem = (*_flow)[e]; 577 if (!_tolerance.positive(rem)) continue; 578 Node v = _graph.source(e); 579 if ((*_level)[v] < level) { 580 if (!_level->active(v) && v != _target) { 581 _level->activate(v); 582 } 583 if (!_tolerance.less(rem, excess)) { 584 _flow->set(e, (*_flow)[e] - excess); 585 _excess->set(v, (*_excess)[v] + excess); 586 excess = 0; 587 goto no_more_push_1; 588 } else { 589 excess -= rem; 590 _excess->set(v, (*_excess)[v] + rem); 591 _flow->set(e, 0); 592 } 593 } else if (new_level > (*_level)[v]) { 594 new_level = (*_level)[v]; 595 } 596 } 597 598 no_more_push_1: 599 600 _excess->set(n, excess); 601 602 if (excess != 0) { 603 if (new_level + 1 < _level->maxLevel()) { 604 _level->liftHighestActive(new_level + 1); 605 } else { 606 _level->liftHighestActiveToTop(); 607 } 608 if (_level->emptyLevel(level)) { 609 _level->liftToTop(level); 610 } 611 } else { 612 _level->deactivate(n); 613 } 614 615 n = _level->highestActive(); 616 level = _level->highestActiveLevel(); 617 --num; 618 } 619 620 num = _node_num * 20; 621 while (num > 0 && n != INVALID) { 622 Value excess = (*_excess)[n]; 623 int new_level = _level->maxLevel(); 624 625 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) { 626 Value rem = (*_capacity)[e] - (*_flow)[e]; 627 if (!_tolerance.positive(rem)) continue; 628 Node v = _graph.target(e); 629 if ((*_level)[v] < level) { 630 if (!_level->active(v) && v != _target) { 631 _level->activate(v); 632 } 633 if (!_tolerance.less(rem, excess)) { 634 _flow->set(e, (*_flow)[e] + excess); 635 _excess->set(v, (*_excess)[v] + excess); 636 excess = 0; 637 goto no_more_push_2; 638 } else { 639 excess -= rem; 640 _excess->set(v, (*_excess)[v] + rem); 641 _flow->set(e, (*_capacity)[e]); 642 } 643 } else if (new_level > (*_level)[v]) { 644 new_level = (*_level)[v]; 645 } 646 } 647 648 for (InEdgeIt e(_graph, n); e != INVALID; ++e) { 649 Value rem = (*_flow)[e]; 650 if (!_tolerance.positive(rem)) continue; 651 Node v = _graph.source(e); 652 if ((*_level)[v] < level) { 653 if (!_level->active(v) && v != _target) { 654 _level->activate(v); 655 } 656 if (!_tolerance.less(rem, excess)) { 657 _flow->set(e, (*_flow)[e] - excess); 658 _excess->set(v, (*_excess)[v] + excess); 659 excess = 0; 660 goto no_more_push_2; 661 } else { 662 excess -= rem; 663 _excess->set(v, (*_excess)[v] + rem); 664 _flow->set(e, 0); 665 } 666 } else if (new_level > (*_level)[v]) { 667 new_level = (*_level)[v]; 668 } 669 } 670 671 no_more_push_2: 672 673 _excess->set(n, excess); 674 675 if (excess != 0) { 676 if (new_level + 1 < _level->maxLevel()) { 677 _level->liftActiveOn(level, new_level + 1); 678 } else { 679 _level->liftActiveToTop(level); 680 } 681 if (_level->emptyLevel(level)) { 682 _level->liftToTop(level); 683 } 684 } else { 685 _level->deactivate(n); 686 } 687 688 while (level >= 0 && _level->activeFree(level)) { 689 --level; 690 } 691 if (level == -1) { 692 n = _level->highestActive(); 693 level = _level->highestActiveLevel(); 694 } else { 695 n = _level->activeOn(level); 696 } 697 --num; 698 } 699 } 700 } 701 702 /// \brief Starts the second phase of the preflow algorithm. 703 /// 704 /// The preflow algorithm consists of two phases, this method runs 705 /// the second phase. After calling \ref init() and \ref 706 /// startFirstPhase() and then \ref startSecondPhase(), \ref 707 /// flowMap() return a maximum flow, \ref flowValue() returns the 708 /// value of a maximum flow, \ref minCut() returns a minimum cut 709 /// \pre The \ref init() and startFirstPhase() functions should be 710 /// called before. 711 void startSecondPhase() { 712 _phase = false; 713 714 typename Graph::template NodeMap<bool> reached(_graph, false); 715 for (NodeIt n(_graph); n != INVALID; ++n) { 716 reached.set(n, (*_level)[n] < _level->maxLevel()); 717 } 718 719 _level->initStart(); 720 _level->initAddItem(_source); 721 722 std::vector<Node> queue; 723 queue.push_back(_source); 724 reached.set(_source, true); 725 726 while (!queue.empty()) { 727 _level->initNewLevel(); 728 std::vector<Node> nqueue; 729 for (int i = 0; i < int(queue.size()); ++i) { 730 Node n = queue[i]; 731 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) { 732 Node v = _graph.target(e); 733 if (!reached[v] && _tolerance.positive((*_flow)[e])) { 734 reached.set(v, true); 735 _level->initAddItem(v); 736 nqueue.push_back(v); 737 } 738 } 739 for (InEdgeIt e(_graph, n); e != INVALID; ++e) { 740 Node u = _graph.source(e); 741 if (!reached[u] && 742 _tolerance.positive((*_capacity)[e] - (*_flow)[e])) { 743 reached.set(u, true); 744 _level->initAddItem(u); 745 nqueue.push_back(u); 746 } 747 } 748 } 749 queue.swap(nqueue); 750 } 751 _level->initFinish(); 752 753 for (NodeIt n(_graph); n != INVALID; ++n) { 754 if ((*_excess)[n] > 0 && _target != n) { 755 _level->activate(n); 756 } 757 } 758 759 Node n; 760 while ((n = _level->highestActive()) != INVALID) { 761 Value excess = (*_excess)[n]; 762 int level = _level->highestActiveLevel(); 763 int new_level = _level->maxLevel(); 764 765 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) { 766 Value rem = (*_capacity)[e] - (*_flow)[e]; 767 if (!_tolerance.positive(rem)) continue; 768 Node v = _graph.target(e); 769 if ((*_level)[v] < level) { 770 if (!_level->active(v) && v != _source) { 771 _level->activate(v); 772 } 773 if (!_tolerance.less(rem, excess)) { 774 _flow->set(e, (*_flow)[e] + excess); 775 _excess->set(v, (*_excess)[v] + excess); 776 excess = 0; 777 goto no_more_push; 778 } else { 779 excess -= rem; 780 _excess->set(v, (*_excess)[v] + rem); 781 _flow->set(e, (*_capacity)[e]); 782 } 783 } else if (new_level > (*_level)[v]) { 784 new_level = (*_level)[v]; 785 } 786 } 787 788 for (InEdgeIt e(_graph, n); e != INVALID; ++e) { 789 Value rem = (*_flow)[e]; 790 if (!_tolerance.positive(rem)) continue; 791 Node v = _graph.source(e); 792 if ((*_level)[v] < level) { 793 if (!_level->active(v) && v != _source) { 794 _level->activate(v); 795 } 796 if (!_tolerance.less(rem, excess)) { 797 _flow->set(e, (*_flow)[e] - excess); 798 _excess->set(v, (*_excess)[v] + excess); 799 excess = 0; 800 goto no_more_push; 801 } else { 802 excess -= rem; 803 _excess->set(v, (*_excess)[v] + rem); 804 _flow->set(e, 0); 805 } 806 } else if (new_level > (*_level)[v]) { 807 new_level = (*_level)[v]; 808 } 809 } 810 811 no_more_push: 812 813 _excess->set(n, excess); 814 815 if (excess != 0) { 816 if (new_level + 1 < _level->maxLevel()) { 817 _level->liftHighestActive(new_level + 1); 818 } else { 819 // Calculation error 820 _level->liftHighestActiveToTop(); 821 } 822 if (_level->emptyLevel(level)) { 823 // Calculation error 824 _level->liftToTop(level); 825 } 826 } else { 827 _level->deactivate(n); 828 } 829 830 } 831 } 832 833 /// \brief Runs the preflow algorithm. 834 /// 835 /// Runs the preflow algorithm. 836 /// \note pf.run() is just a shortcut of the following code. 837 /// \code 838 /// pf.init(); 839 /// pf.startFirstPhase(); 840 /// pf.startSecondPhase(); 841 /// \endcode 842 void run() { 843 init(); 844 startFirstPhase(); 845 startSecondPhase(); 846 } 847 848 /// \brief Runs the preflow algorithm to compute the minimum cut. 849 /// 850 /// Runs the preflow algorithm to compute the minimum cut. 851 /// \note pf.runMinCut() is just a shortcut of the following code. 852 /// \code 853 /// pf.init(); 854 /// pf.startFirstPhase(); 855 /// \endcode 856 void runMinCut() { 857 init(); 858 startFirstPhase(); 859 } 860 861 /// @} 862 863 /// \name Query Functions 864 /// The result of the %Dijkstra algorithm can be obtained using these 865 /// functions.\n 866 /// Before the use of these functions, 867 /// either run() or start() must be called. 110 868 111 ///Indicates the property of the starting flow map. 112 113 ///Indicates the property of the starting flow map. 114 /// 115 enum FlowEnum{ 116 ///indicates an unspecified edge map. \c flow will be 117 ///set to the constant zero flow in the beginning of 118 ///the algorithm in this case. 119 NO_FLOW, 120 ///constant zero flow 121 ZERO_FLOW, 122 ///any flow, i.e. the sum of the in-flows equals to 123 ///the sum of the out-flows in every node except the \c source and 124 ///the \c target. 125 GEN_FLOW, 126 ///any preflow, i.e. the sum of the in-flows is at 127 ///least the sum of the out-flows in every node except the \c source. 128 PRE_FLOW 129 }; 130 131 ///Indicates the state of the preflow algorithm. 132 133 ///Indicates the state of the preflow algorithm. 134 /// 135 enum StatusEnum { 136 ///before running the algorithm or 137 ///at an unspecified state. 138 AFTER_NOTHING, 139 ///right after running \ref phase1() 140 AFTER_PREFLOW_PHASE_1, 141 ///after running \ref phase2() 142 AFTER_PREFLOW_PHASE_2 143 }; 144 145 protected: 146 FlowEnum flow_prop; 147 StatusEnum status; // Do not needle this flag only if necessary. 148 149 public: 150 ///The constructor of the class. 151 152 ///The constructor of the class. 153 ///\param _gr The directed graph the algorithm runs on. 154 ///\param _s The source node. 155 ///\param _t The target node. 156 ///\param _cap The capacity of the edges. 157 ///\param _f The flow of the edges. 158 ///\param _sr Tol class. 159 ///Except the graph, all of these parameters can be reset by 160 ///calling \ref source, \ref target, \ref capacityMap and \ref 161 ///flowMap, resp. 162 Preflow(const Graph& _gr, Node _s, Node _t, 163 const CapacityMap& _cap, FlowMap& _f, 164 const Tol &_sr=Tol()) : 165 _g(&_gr), _source(_s), _target(_t), _capacity(&_cap), 166 _flow(&_f), _surely(_sr), 167 _node_num(countNodes(_gr)), level(_gr), excess(_gr,0), 168 flow_prop(NO_FLOW), status(AFTER_NOTHING) { 169 if ( _source==_target ) 170 throw InvalidArgument(); 171 } 172 173 ///Give a reference to the tolerance handler class 174 175 ///Give a reference to the tolerance handler class 176 ///\sa Tolerance 177 Tol &tolerance() { return _surely; } 178 179 ///Runs the preflow algorithm. 180 181 ///Runs the preflow algorithm. 182 /// 183 void run() { 184 phase1(flow_prop); 185 phase2(); 186 } 187 188 ///Runs the preflow algorithm. 189 190 ///Runs the preflow algorithm. 191 ///\pre The starting flow map must be 192 /// - a constant zero flow if \c fp is \c ZERO_FLOW, 193 /// - an arbitrary flow if \c fp is \c GEN_FLOW, 194 /// - an arbitrary preflow if \c fp is \c PRE_FLOW, 195 /// - any map if \c fp is NO_FLOW. 196 ///If the starting flow map is a flow or a preflow then 197 ///the algorithm terminates faster. 198 void run(FlowEnum fp) { 199 flow_prop=fp; 200 run(); 201 } 202 203 ///Runs the first phase of the preflow algorithm. 204 205 ///The preflow algorithm consists of two phases, this method runs 206 ///the first phase. After the first phase the maximum flow value 207 ///and a minimum value cut can already be computed, although a 208 ///maximum flow is not yet obtained. So after calling this method 209 ///\ref flowValue returns the value of a maximum flow and \ref 210 ///minCut returns a minimum cut. 211 ///\warning \ref minMinCut and \ref maxMinCut do not give minimum 212 ///value cuts unless calling \ref phase2. 213 ///\warning A real flow map (i.e. not \ref lemon::NullMap "NullMap") 214 ///is needed for this phase. 215 ///\pre The starting flow must be 216 ///- a constant zero flow if \c fp is \c ZERO_FLOW, 217 ///- an arbitary flow if \c fp is \c GEN_FLOW, 218 ///- an arbitary preflow if \c fp is \c PRE_FLOW, 219 ///- any map if \c fp is NO_FLOW. 220 void phase1(FlowEnum fp) 221 { 222 flow_prop=fp; 223 phase1(); 224 } 225 226 227 ///Runs the first phase of the preflow algorithm. 228 229 ///The preflow algorithm consists of two phases, this method runs 230 ///the first phase. After the first phase the maximum flow value 231 ///and a minimum value cut can already be computed, although a 232 ///maximum flow is not yet obtained. So after calling this method 233 ///\ref flowValue returns the value of a maximum flow and \ref 234 ///minCut returns a minimum cut. 235 ///\warning \ref minMinCut() and \ref maxMinCut() do not 236 ///give minimum value cuts unless calling \ref phase2(). 237 ///\warning A real flow map (i.e. not \ref lemon::NullMap "NullMap") 238 ///is needed for this phase. 239 void phase1() 240 { 241 int heur0=int(H0*_node_num); //time while running 'bound decrease' 242 int heur1=int(H1*_node_num); //time while running 'highest label' 243 int heur=heur1; //starting time interval (#of relabels) 244 int numrelabel=0; 245 246 bool what_heur=1; 247 //It is 0 in case 'bound decrease' and 1 in case 'highest label' 248 249 bool end=false; 250 //Needed for 'bound decrease', true means no active 251 //nodes are above bound b. 252 253 int k=_node_num-2; //bound on the highest level under n containing a node 254 int b=k; //bound on the highest level under n containing an active node 255 256 VecNode first(_node_num, INVALID); 257 NNMap next(*_g, INVALID); 258 259 NNMap left(*_g, INVALID); 260 NNMap right(*_g, INVALID); 261 VecNode level_list(_node_num,INVALID); 262 //List of the nodes in level i<n, set to n. 263 264 preflowPreproc(first, next, level_list, left, right); 265 266 //Push/relabel on the highest level active nodes. 267 while ( true ) { 268 if ( b == 0 ) { 269 if ( !what_heur && !end && k > 0 ) { 270 b=k; 271 end=true; 272 } else break; 273 } 274 275 if ( first[b]==INVALID ) --b; 276 else { 277 end=false; 278 Node w=first[b]; 279 first[b]=next[w]; 280 int newlevel=push(w, next, first); 281 if ( excess[w] != 0 ) { 282 relabel(w, newlevel, first, next, level_list, 283 left, right, b, k, what_heur); 284 } 285 286 ++numrelabel; 287 if ( numrelabel >= heur ) { 288 numrelabel=0; 289 if ( what_heur ) { 290 what_heur=0; 291 heur=heur0; 292 end=false; 293 } else { 294 what_heur=1; 295 heur=heur1; 296 b=k; 297 } 298 } 299 } 300 } 301 flow_prop=PRE_FLOW; 302 status=AFTER_PREFLOW_PHASE_1; 303 } 304 // Heuristics: 305 // 2 phase 306 // gap 307 // list 'level_list' on the nodes on level i implemented by hand 308 // stack 'active' on the active nodes on level i 309 // runs heuristic 'highest label' for H1*n relabels 310 // runs heuristic 'bound decrease' for H0*n relabels, 311 // starts with 'highest label' 312 // Parameters H0 and H1 are initialized to 20 and 1. 313 314 315 ///Runs the second phase of the preflow algorithm. 316 317 ///The preflow algorithm consists of two phases, this method runs 318 ///the second phase. After calling \ref phase1() and then 319 ///\ref phase2(), 320 /// \ref flowMap() return a maximum flow, \ref flowValue 321 ///returns the value of a maximum flow, \ref minCut returns a 322 ///minimum cut, while the methods \ref minMinCut and \ref 323 ///maxMinCut return the inclusionwise minimum and maximum cuts of 324 ///minimum value, resp. \pre \ref phase1 must be called before. 325 /// 326 /// \todo The inexact computation can cause positive excess on a set of 327 /// unpushable nodes. We may have to watch the empty level in this case 328 /// due to avoid the terrible long running time. 329 void phase2() 330 { 331 332 int k=_node_num-2; //bound on the highest level under n containing a node 333 int b=k; //bound on the highest level under n of an active node 334 335 336 VecNode first(_node_num, INVALID); 337 NNMap next(*_g, INVALID); 338 level.set(_source,0); 339 std::queue<Node> bfs_queue; 340 bfs_queue.push(_source); 341 342 while ( !bfs_queue.empty() ) { 343 344 Node v=bfs_queue.front(); 345 bfs_queue.pop(); 346 int l=level[v]+1; 347 348 for(InEdgeIt e(*_g,v); e!=INVALID; ++e) { 349 if ( !_surely.positive((*_capacity)[e] - (*_flow)[e])) continue; 350 Node u=_g->source(e); 351 if ( level[u] >= _node_num ) { 352 bfs_queue.push(u); 353 level.set(u, l); 354 if ( excess[u] != 0 ) { 355 next.set(u,first[l]); 356 first[l]=u; 357 } 358 } 359 } 360 361 for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) { 362 if ( !_surely.positive((*_flow)[e]) ) continue; 363 Node u=_g->target(e); 364 if ( level[u] >= _node_num ) { 365 bfs_queue.push(u); 366 level.set(u, l); 367 if ( excess[u] != 0 ) { 368 next.set(u,first[l]); 369 first[l]=u; 370 } 371 } 372 } 373 } 374 b=_node_num-2; 375 376 while ( true ) { 377 378 if ( b == 0 ) break; 379 if ( first[b]==INVALID ) --b; 380 else { 381 Node w=first[b]; 382 first[b]=next[w]; 383 int newlevel=push(w,next, first); 384 385 if ( newlevel == _node_num) { 386 excess.set(w, 0); 387 level.set(w,_node_num); 388 } 389 //relabel 390 if ( excess[w] != 0 ) { 391 level.set(w,++newlevel); 392 next.set(w,first[newlevel]); 393 first[newlevel]=w; 394 b=newlevel; 395 } 396 } 397 } // while(true) 398 flow_prop=GEN_FLOW; 399 status=AFTER_PREFLOW_PHASE_2; 400 } 401 402 /// Returns the value of the maximum flow. 403 869 ///@{ 870 871 /// \brief Returns the value of the maximum flow. 872 /// 404 873 /// Returns the value of the maximum flow by returning the excess 405 874 /// of the target node \c t. This value equals to the value of 406 /// the maximum flow already after running \ref phase1. 407 Num flowValue() const { 408 return excess[_target]; 409 } 410 411 412 ///Returns a minimum value cut. 413 414 ///Sets \c M to the characteristic vector of a minimum value 415 ///cut. This method can be called both after running \ref 416 ///phase1 and \ref phase2. It is much faster after 417 ///\ref phase1. \pre M should be a bool-valued node-map. \pre 418 ///If \ref minCut() is called after \ref phase2() then M should 419 ///be initialized to false. 420 template<typename _CutMap> 421 void minCut(_CutMap& M) const { 422 switch ( status ) { 423 case AFTER_PREFLOW_PHASE_1: 424 for(NodeIt v(*_g); v!=INVALID; ++v) { 425 if (level[v] < _node_num) { 426 M.set(v, false); 427 } else { 428 M.set(v, true); 429 } 430 } 431 break; 432 case AFTER_PREFLOW_PHASE_2: 433 minMinCut(M); 434 break; 435 case AFTER_NOTHING: 436 break; 437 } 438 } 439 440 ///Returns the inclusionwise minimum of the minimum value cuts. 441 442 ///Sets \c M to the characteristic vector of the minimum value cut 443 ///which is inclusionwise minimum. It is computed by processing a 444 ///bfs from the source node \c s in the residual graph. \pre M 445 ///should be a node map of bools initialized to false. \pre \ref 446 ///phase2 should already be run. 447 template<typename _CutMap> 448 void minMinCut(_CutMap& M) const { 449 450 std::queue<Node> queue; 451 M.set(_source,true); 452 queue.push(_source); 453 454 while (!queue.empty()) { 455 Node w=queue.front(); 456 queue.pop(); 457 458 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { 459 Node v=_g->target(e); 460 if (!M[v] && _surely.positive((*_capacity)[e] -(*_flow)[e]) ) { 461 queue.push(v); 462 M.set(v, true); 463 } 464 } 465 466 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { 467 Node v=_g->source(e); 468 if (!M[v] && _surely.positive((*_flow)[e]) ) { 469 queue.push(v); 470 M.set(v, true); 471 } 472 } 473 } 875 /// the maximum flow already after the first phase. 876 Value flowValue() const { 877 return (*_excess)[_target]; 878 } 879 880 /// \brief Returns true when the node is on the source side of minimum cut. 881 /// 882 /// Returns true when the node is on the source side of minimum 883 /// cut. This method can be called both after running \ref 884 /// startFirstPhase() and \ref startSecondPhase(). 885 bool minCut(const Node& node) const { 886 return ((*_level)[node] == _level->maxLevel()) == _phase; 887 } 888 889 /// \brief Returns a minimum value cut. 890 /// 891 /// Sets the \c cutMap to the characteristic vector of a minimum value 892 /// cut. This method can be called both after running \ref 893 /// startFirstPhase() and \ref startSecondPhase(). The result after second 894 /// phase could be changed slightly if inexact computation is used. 895 /// \pre The \c cutMap should be a bool-valued node-map. 896 template <typename CutMap> 897 void minCutMap(CutMap& cutMap) const { 898 for (NodeIt n(_graph); n != INVALID; ++n) { 899 cutMap.set(n, minCut(n)); 900 } 901 } 902 903 /// \brief Returns the flow on the edge. 904 /// 905 /// Sets the \c flowMap to the flow on the edges. This method can 906 /// be called after the second phase of algorithm. 907 Value flow(const Edge& edge) const { 908 return (*_flow)[edge]; 474 909 } 475 910 476 ///Returns the inclusionwise maximum of the minimum value cuts. 477 478 ///Sets \c M to the characteristic vector of the minimum value cut 479 ///which is inclusionwise maximum. It is computed by processing a 480 ///backward bfs from the target node \c t in the residual graph. 481 ///\pre \ref phase2() or run() should already be run. 482 template<typename _CutMap> 483 void maxMinCut(_CutMap& M) const { 484 485 for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true); 486 487 std::queue<Node> queue; 488 489 M.set(_target,false); 490 queue.push(_target); 491 492 while (!queue.empty()) { 493 Node w=queue.front(); 494 queue.pop(); 495 496 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { 497 Node v=_g->source(e); 498 if (M[v] && _surely.positive((*_capacity)[e] - (*_flow)[e]) ) { 499 queue.push(v); 500 M.set(v, false); 501 } 502 } 503 504 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { 505 Node v=_g->target(e); 506 if (M[v] && _surely.positive((*_flow)[e]) ) { 507 queue.push(v); 508 M.set(v, false); 509 } 510 } 511 } 512 } 513 514 ///Sets the source node to \c _s. 515 516 ///Sets the source node to \c _s. 517 /// 518 void source(Node _s) { 519 _source=_s; 520 if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW; 521 status=AFTER_NOTHING; 522 } 523 524 ///Returns the source node. 525 526 ///Returns the source node. 527 /// 528 Node source() const { 529 return _source; 530 } 531 532 ///Sets the target node to \c _t. 533 534 ///Sets the target node to \c _t. 535 /// 536 void target(Node _t) { 537 _target=_t; 538 if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW; 539 status=AFTER_NOTHING; 540 } 541 542 ///Returns the target node. 543 544 ///Returns the target node. 545 /// 546 Node target() const { 547 return _target; 548 } 549 550 /// Sets the edge map of the capacities to _cap. 551 552 /// Sets the edge map of the capacities to _cap. 553 /// 554 void capacityMap(const CapacityMap& _cap) { 555 _capacity=&_cap; 556 status=AFTER_NOTHING; 557 } 558 /// Returns a reference to capacity map. 559 560 /// Returns a reference to capacity map. 561 /// 562 const CapacityMap &capacityMap() const { 563 return *_capacity; 564 } 565 566 /// Sets the edge map of the flows to _flow. 567 568 /// Sets the edge map of the flows to _flow. 569 /// 570 void flowMap(FlowMap& _f) { 571 _flow=&_f; 572 flow_prop=NO_FLOW; 573 status=AFTER_NOTHING; 574 } 575 576 /// Returns a reference to flow map. 577 578 /// Returns a reference to flow map. 579 /// 580 const FlowMap &flowMap() const { 581 return *_flow; 582 } 583 584 private: 585 586 int push(Node w, NNMap& next, VecNode& first) { 587 588 int lev=level[w]; 589 Num exc=excess[w]; 590 int newlevel=_node_num; //bound on the next level of w 591 592 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { 593 if ( !_surely.positive((*_capacity)[e] - (*_flow)[e])) continue; 594 Node v=_g->target(e); 595 596 if( lev > level[v] ) { //Push is allowed now 597 598 if ( excess[v] == 0 && v!=_target && v!=_source ) { 599 next.set(v,first[level[v]]); 600 first[level[v]]=v; 601 } 602 603 Num cap=(*_capacity)[e]; 604 Num flo=(*_flow)[e]; 605 Num remcap=cap-flo; 606 607 if ( ! _surely.less(remcap, exc) ) { //A nonsaturating push. 608 609 _flow->set(e, flo+exc); 610 excess.set(v, excess[v]+exc); 611 exc=0; 612 break; 613 614 } else { //A saturating push. 615 _flow->set(e, cap); 616 excess.set(v, excess[v]+remcap); 617 exc-=remcap; 618 } 619 } else if ( newlevel > level[v] ) newlevel = level[v]; 620 } //for out edges wv 621 622 if ( exc != 0 ) { 623 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { 624 625 if ( !_surely.positive((*_flow)[e]) ) continue; 626 Node v=_g->source(e); 627 628 if( lev > level[v] ) { //Push is allowed now 629 630 if ( excess[v] == 0 && v!=_target && v!=_source ) { 631 next.set(v,first[level[v]]); 632 first[level[v]]=v; 633 } 634 635 Num flo=(*_flow)[e]; 636 637 if ( !_surely.less(flo, exc) ) { //A nonsaturating push. 638 639 _flow->set(e, flo-exc); 640 excess.set(v, excess[v]+exc); 641 exc=0; 642 break; 643 } else { //A saturating push. 644 645 excess.set(v, excess[v]+flo); 646 exc-=flo; 647 _flow->set(e,0); 648 } 649 } else if ( newlevel > level[v] ) newlevel = level[v]; 650 } //for in edges vw 651 652 } // if w still has excess after the out edge for cycle 653 654 excess.set(w, exc); 655 656 return newlevel; 657 } 658 659 660 661 void preflowPreproc(VecNode& first, NNMap& next, 662 VecNode& level_list, NNMap& left, NNMap& right) 663 { 664 for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num); 665 std::queue<Node> bfs_queue; 666 667 if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) { 668 //Reverse_bfs from t in the residual graph, 669 //to find the starting level. 670 level.set(_target,0); 671 bfs_queue.push(_target); 672 673 while ( !bfs_queue.empty() ) { 674 675 Node v=bfs_queue.front(); 676 bfs_queue.pop(); 677 int l=level[v]+1; 678 679 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) { 680 if ( !_surely.positive((*_capacity)[e] - (*_flow)[e] )) continue; 681 Node w=_g->source(e); 682 if ( level[w] == _node_num && w != _source ) { 683 bfs_queue.push(w); 684 Node z=level_list[l]; 685 if ( z!=INVALID ) left.set(z,w); 686 right.set(w,z); 687 level_list[l]=w; 688 level.set(w, l); 689 } 690 } 691 692 for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) { 693 if ( !_surely.positive((*_flow)[e]) ) continue; 694 Node w=_g->target(e); 695 if ( level[w] == _node_num && w != _source ) { 696 bfs_queue.push(w); 697 Node z=level_list[l]; 698 if ( z!=INVALID ) left.set(z,w); 699 right.set(w,z); 700 level_list[l]=w; 701 level.set(w, l); 702 } 703 } 704 } //while 705 } //if 706 707 708 switch (flow_prop) { 709 case NO_FLOW: 710 for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0); 711 case ZERO_FLOW: 712 for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0); 713 714 //Reverse_bfs from t, to find the starting level. 715 level.set(_target,0); 716 bfs_queue.push(_target); 717 718 while ( !bfs_queue.empty() ) { 719 720 Node v=bfs_queue.front(); 721 bfs_queue.pop(); 722 int l=level[v]+1; 723 724 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) { 725 Node w=_g->source(e); 726 if ( level[w] == _node_num && w != _source ) { 727 bfs_queue.push(w); 728 Node z=level_list[l]; 729 if ( z!=INVALID ) left.set(z,w); 730 right.set(w,z); 731 level_list[l]=w; 732 level.set(w, l); 733 } 734 } 735 } 736 737 //the starting flow 738 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { 739 Num c=(*_capacity)[e]; 740 if ( !_surely.positive(c) ) continue; 741 Node w=_g->target(e); 742 if ( level[w] < _node_num ) { 743 if ( excess[w] == 0 && w!=_target ) { //putting into the stack 744 next.set(w,first[level[w]]); 745 first[level[w]]=w; 746 } 747 _flow->set(e, c); 748 excess.set(w, excess[w]+c); 749 } 750 } 751 break; 752 753 case GEN_FLOW: 754 for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0); 755 { 756 Num exc=0; 757 for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e]; 758 for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e]; 759 if (!_surely.positive(exc)) { 760 exc = 0; 761 } 762 excess.set(_target,exc); 763 } 764 765 //the starting flow 766 for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e) { 767 Num rem=(*_capacity)[e]-(*_flow)[e]; 768 if ( !_surely.positive(rem) ) continue; 769 Node w=_g->target(e); 770 if ( level[w] < _node_num ) { 771 if ( excess[w] == 0 && w!=_target ) { //putting into the stack 772 next.set(w,first[level[w]]); 773 first[level[w]]=w; 774 } 775 _flow->set(e, (*_capacity)[e]); 776 excess.set(w, excess[w]+rem); 777 } 778 } 779 780 for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) { 781 if ( !_surely.positive((*_flow)[e]) ) continue; 782 Node w=_g->source(e); 783 if ( level[w] < _node_num ) { 784 if ( excess[w] == 0 && w!=_target ) { 785 next.set(w,first[level[w]]); 786 first[level[w]]=w; 787 } 788 excess.set(w, excess[w]+(*_flow)[e]); 789 _flow->set(e, 0); 790 } 791 } 792 break; 793 794 case PRE_FLOW: 795 //the starting flow 796 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { 797 Num rem=(*_capacity)[e]-(*_flow)[e]; 798 if ( !_surely.positive(rem) ) continue; 799 Node w=_g->target(e); 800 if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]); 801 } 802 803 for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { 804 if ( !_surely.positive((*_flow)[e]) ) continue; 805 Node w=_g->source(e); 806 if ( level[w] < _node_num ) _flow->set(e, 0); 807 } 808 809 //computing the excess 810 for(NodeIt w(*_g); w!=INVALID; ++w) { 811 Num exc=0; 812 for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e]; 813 for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e]; 814 if (!_surely.positive(exc)) { 815 exc = 0; 816 } 817 excess.set(w,exc); 818 819 //putting the active nodes into the stack 820 int lev=level[w]; 821 if ( exc != 0 && lev < _node_num && Node(w) != _target ) { 822 next.set(w,first[lev]); 823 first[lev]=w; 824 } 825 } 826 break; 827 } //switch 828 } //preflowPreproc 829 830 831 void relabel(Node w, int newlevel, VecNode& first, NNMap& next, 832 VecNode& level_list, NNMap& left, 833 NNMap& right, int& b, int& k, bool what_heur ) 834 { 835 836 int lev=level[w]; 837 838 Node right_n=right[w]; 839 Node left_n=left[w]; 840 841 //unlacing starts 842 if ( right_n!=INVALID ) { 843 if ( left_n!=INVALID ) { 844 right.set(left_n, right_n); 845 left.set(right_n, left_n); 846 } else { 847 level_list[lev]=right_n; 848 left.set(right_n, INVALID); 849 } 850 } else { 851 if ( left_n!=INVALID ) { 852 right.set(left_n, INVALID); 853 } else { 854 level_list[lev]=INVALID; 855 } 856 } 857 //unlacing ends 858 859 if ( level_list[lev]==INVALID ) { 860 861 //gapping starts 862 for (int i=lev; i!=k ; ) { 863 Node v=level_list[++i]; 864 while ( v!=INVALID ) { 865 level.set(v,_node_num); 866 v=right[v]; 867 } 868 level_list[i]=INVALID; 869 if ( !what_heur ) first[i]=INVALID; 870 } 871 872 level.set(w,_node_num); 873 b=lev-1; 874 k=b; 875 //gapping ends 876 877 } else { 878 879 if ( newlevel == _node_num ) level.set(w,_node_num); 880 else { 881 level.set(w,++newlevel); 882 next.set(w,first[newlevel]); 883 first[newlevel]=w; 884 if ( what_heur ) b=newlevel; 885 if ( k < newlevel ) ++k; //now k=newlevel 886 Node z=level_list[newlevel]; 887 if ( z!=INVALID ) left.set(z,w); 888 right.set(w,z); 889 left.set(w,INVALID); 890 level_list[newlevel]=w; 891 } 892 } 893 } //relabel 894 895 }; 896 897 ///\ingroup max_flow 898 ///\brief Function type interface for Preflow algorithm. 899 /// 900 ///Function type interface for Preflow algorithm. 901 ///\sa Preflow 902 template<class GR, class CM, class FM> 903 Preflow<GR,typename CM::Value,CM,FM> preflow(const GR &g, 904 typename GR::Node source, 905 typename GR::Node target, 906 const CM &cap, 907 FM &flow 908 ) 909 { 910 return Preflow<GR,typename CM::Value,CM,FM>(g,source,target,cap,flow); 911 } 912 913 } //namespace lemon 914 915 #endif //LEMON_PREFLOW_H 911 /// @} 912 }; 913 } 914 915 #endif
Note: See TracChangeset
for help on using the changeset viewer.