Changeset 485:7f461ab4af1a in lemon0.x for src/work/jacint
 Timestamp:
 04/29/04 19:34:42 (17 years ago)
 Branch:
 default
 Phase:
 public
 Convert:
 svn:c9d7d8f590d60310b91f818b3a526b0e/lemon/trunk@642
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

src/work/jacint/max_flow.h
r480 r485 2 2 3 3 /* 4 Heuristics:5 2 phase6 gap7 list 'level_list' on the nodes on level i implemented by hand8 stack 'active' on the active nodes on level i9 runs heuristic 'highest label' for H1*n relabels10 runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'4 Heuristics: 5 2 phase 6 gap 7 list 'level_list' on the nodes on level i implemented by hand 8 stack 'active' on the active nodes on level i 9 runs heuristic 'highest label' for H1*n relabels 10 runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label' 11 11 12 Parameters H0 and H1 are initialized to 20 and 1.13 14 Constructors:15 16 Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if17 18 19 Members:20 21 void run()22 23 Num flowValue() : returns the value of a maximum flow24 25 void minMinCut(CutMap& M) : sets M to the characteristic vector of the26 27 28 void maxMinCut(CutMap& M) : sets M to the characteristic vector of the29 30 31 void minCut(CutMap& M) : sets M to the characteristic vector of32 12 Parameters H0 and H1 are initialized to 20 and 1. 13 14 Constructors: 15 16 Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if 17 FlowMap is not constant zero, and should be true if it is 18 19 Members: 20 21 void run() 22 23 Num flowValue() : returns the value of a maximum flow 24 25 void minMinCut(CutMap& M) : sets M to the characteristic vector of the 26 minimum min cut. M should be a map of bools initialized to false. ??Is it OK? 27 28 void maxMinCut(CutMap& M) : sets M to the characteristic vector of the 29 maximum min cut. M should be a map of bools initialized to false. 30 31 void minCut(CutMap& M) : sets M to the characteristic vector of 32 a min cut. M should be a map of bools initialized to false. 33 33 34 34 */ … … 53 53 namespace hugo { 54 54 55 ///\author Marton Makai, Jacint Szabo 55 56 template <typename Graph, typename Num, 56 57 typename CapMap=typename Graph::template EdgeMap<Num>, … … 98 99 flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {} 99 100 101 /// A max flow algorithm is run. 102 ///\pre the flow have to be 0 at the beginning. 100 103 void run() { 101 104 preflow( ZERO_FLOW ); 102 105 } 103 106 107 /// A preflow algorithm is run. The initial edgeset have to be a flow, 108 /// or from a preflow, according to \c fe. 104 109 void preflow( flowEnum fe ) { 105 110 preflowPhase0(fe); … … 107 112 } 108 113 114 /// Run the first phase of preflow, starting from a 0 flow, from a flow, 115 /// or from a preflow, according to \c fe. 109 116 void preflowPhase0( flowEnum fe ); 110 117 118 /// Second phase of preflow. 111 119 void preflowPhase1(); 112 120 121 /// Starting from a flow, this method searches for an augmenting path 122 /// according to the EdmondsKarp algorithm 123 /// and augments the flow on if any. 113 124 bool augmentOnShortestPath(); 114 125 126 /// Starting from a flow, this method searches for an augmenting blockin 127 /// flow according to Dinits' algorithm and augments the flow on if any. 128 /// The blocking flow is computed in a physically constructed 129 /// residual graph of type \c Mutablegraph. 115 130 template<typename MutableGraph> bool augmentOnBlockingFlow(); 116 131 132 /// The same as \c augmentOnBlockingFlow<MutableGraph> but the 133 /// residual graph is not constructed physically. 117 134 bool augmentOnBlockingFlow2(); 118 135 … … 127 144 } 128 145 129 //should be used only between preflowPhase0 and preflowPhase1 146 /// Should be used between preflowPhase0 and preflowPhase1. 147 ///\todo We have to make some status variable which shows the actual state 148 /// of the class. This enables us to determine which methods are valid 149 /// for MinCut computation 130 150 template<typename _CutMap> 131 151 void actMinCut(_CutMap& M) { 132 152 NodeIt v; 133 for(g>first(v); g>valid(v); g>next(v)) 134 135 M.set(v,false);136 137 M.set(v,true);138 139 }140 141 142 143 / *144 Returns the minimum min cut, bya bfs from s in the residual graph.145 */153 for(g>first(v); g>valid(v); g>next(v)) { 154 if ( level[v] < n ) { 155 M.set(v,false); 156 } else { 157 M.set(v,true); 158 } 159 } 160 } 161 162 163 /// The unique inclusionwise minimum cut is computed by 164 /// processing a bfs from s in the residual graph. 165 ///\pre flow have to be a max flow otherwise it will the whole nodeset. 146 166 template<typename _CutMap> 147 167 void minMinCut(_CutMap& M) { … … 177 197 178 198 179 180 /* 181 Returns the maximum min cut, by a reverse bfs 182 from t in the residual graph. 183 */ 184 199 /// The unique inclusionwise maximum cut is computed by 200 /// processing a reverse bfs from t in the residual graph. 201 ///\pre flow have to be a max flow otherwise it will be empty. 185 202 template<typename _CutMap> 186 203 void maxMinCut(_CutMap& M) { … … 222 239 223 240 241 /// A minimum cut is computed. 224 242 template<typename CutMap> 225 void minCut(CutMap& M) { 226 minMinCut(M); 227 } 228 243 void minCut(CutMap& M) { minMinCut(M); } 244 245 /// 246 void resetSource(Node _s) {s=_s;} 247 /// 229 248 void resetTarget(Node _t) {t=_t;} 230 void resetSource(Node _s) {s=_s;}231 249 232 void resetCap(const CapMap& _cap) { 233 capacity=&_cap; 234 } 250 /// capacitymap is changed. 251 void resetCap(const CapMap& _cap) { capacity=&_cap; } 235 252 236 void resetFlow(FlowMap& _flow) { 237 flow=&_flow; 238 } 253 /// flowmap is changed. 254 void resetFlow(FlowMap& _flow) { flow=&_flow; } 239 255 240 256 … … 315 331 316 332 return newlevel; 317 333 } 318 334 319 335 … … 321 337 VecNode& level_list, NNMap& left, NNMap& right ) { 322 338 323 324 325 326 327 {328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 bfs_queue.push(w);343 Node first=level_list[l];344 if ( g>valid(first) ) left.set(first,w);345 right.set(w,first);346 level_list[l]=w;347 level.set(w, l);348 349 350 351 352 353 354 355 356 357 358 359 360 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);361 flow>set(e, c);362 excess.set(w, excess[w]+c);363 364 365 366 }367 368 369 370 {371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 bfs_queue.push(w);388 Node first=level_list[l];389 if ( g>valid(first) ) left.set(first,w);390 right.set(w,first);391 level_list[l]=w;392 level.set(w, l);393 394 395 396 397 398 399 400 401 bfs_queue.push(w);402 Node first=level_list[l];403 if ( g>valid(first) ) left.set(first,w);404 right.set(w,first);405 level_list[l]=w;406 level.set(w, l);407 408 409 410 411 412 413 414 415 416 417 418 419 420 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);421 flow>set(e, (*capacity)[e]);422 excess.set(w, excess[w]+rem);423 424 425 426 427 428 429 430 431 432 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);433 excess.set(w, excess[w]+(*flow)[f]);434 flow>set(f, 0);435 436 437 438 } //case PREFLOW439 440 339 std::queue<Node> bfs_queue; 340 341 switch ( fe ) { 342 case ZERO_FLOW: 343 { 344 //Reverse_bfs from t, to find the starting level. 345 level.set(t,0); 346 bfs_queue.push(t); 347 348 while (!bfs_queue.empty()) { 349 350 Node v=bfs_queue.front(); 351 bfs_queue.pop(); 352 int l=level[v]+1; 353 354 InEdgeIt e; 355 for(g>first(e,v); g>valid(e); g>next(e)) { 356 Node w=g>tail(e); 357 if ( level[w] == n && w != s ) { 358 bfs_queue.push(w); 359 Node first=level_list[l]; 360 if ( g>valid(first) ) left.set(first,w); 361 right.set(w,first); 362 level_list[l]=w; 363 level.set(w, l); 364 } 365 } 366 } 367 368 //the starting flow 369 OutEdgeIt e; 370 for(g>first(e,s); g>valid(e); g>next(e)) 371 { 372 Num c=(*capacity)[e]; 373 if ( c <= 0 ) continue; 374 Node w=g>head(e); 375 if ( level[w] < n ) { 376 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w); 377 flow>set(e, c); 378 excess.set(w, excess[w]+c); 379 } 380 } 381 break; 382 } 383 384 case GEN_FLOW: 385 case PREFLOW: 386 { 387 //Reverse_bfs from t in the residual graph, 388 //to find the starting level. 389 level.set(t,0); 390 bfs_queue.push(t); 391 392 while (!bfs_queue.empty()) { 393 394 Node v=bfs_queue.front(); 395 bfs_queue.pop(); 396 int l=level[v]+1; 397 398 InEdgeIt e; 399 for(g>first(e,v); g>valid(e); g>next(e)) { 400 if ( (*capacity)[e] <= (*flow)[e] ) continue; 401 Node w=g>tail(e); 402 if ( level[w] == n && w != s ) { 403 bfs_queue.push(w); 404 Node first=level_list[l]; 405 if ( g>valid(first) ) left.set(first,w); 406 right.set(w,first); 407 level_list[l]=w; 408 level.set(w, l); 409 } 410 } 411 412 OutEdgeIt f; 413 for(g>first(f,v); g>valid(f); g>next(f)) { 414 if ( 0 >= (*flow)[f] ) continue; 415 Node w=g>head(f); 416 if ( level[w] == n && w != s ) { 417 bfs_queue.push(w); 418 Node first=level_list[l]; 419 if ( g>valid(first) ) left.set(first,w); 420 right.set(w,first); 421 level_list[l]=w; 422 level.set(w, l); 423 } 424 } 425 } 426 427 428 //the starting flow 429 OutEdgeIt e; 430 for(g>first(e,s); g>valid(e); g>next(e)) 431 { 432 Num rem=(*capacity)[e](*flow)[e]; 433 if ( rem <= 0 ) continue; 434 Node w=g>head(e); 435 if ( level[w] < n ) { 436 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w); 437 flow>set(e, (*capacity)[e]); 438 excess.set(w, excess[w]+rem); 439 } 440 } 441 442 InEdgeIt f; 443 for(g>first(f,s); g>valid(f); g>next(f)) 444 { 445 if ( (*flow)[f] <= 0 ) continue; 446 Node w=g>tail(f); 447 if ( level[w] < n ) { 448 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w); 449 excess.set(w, excess[w]+(*flow)[f]); 450 flow>set(f, 0); 451 } 452 } 453 break; 454 } //case PREFLOW 455 } 456 } //preflowPreproc 441 457 442 458 … … 522 538 } 523 539 int operator[](const typename MapGraphWrapper::Node& n) 524 525 // int get(const typename MapGraphWrapper::Node& n) const {526 // return dist[n]; }527 // bool get(const typename MapGraphWrapper::Edge& e) const {528 // return (dist.get(g>tail(e))<dist.get(g>head(e))); }540 { return dist[n]; } 541 // int get(const typename MapGraphWrapper::Node& n) const { 542 // return dist[n]; } 543 // bool get(const typename MapGraphWrapper::Edge& e) const { 544 // return (dist.get(g>tail(e))<dist.get(g>head(e))); } 529 545 bool operator[](const typename MapGraphWrapper::Edge& e) const { 530 546 return (dist[g>tail(e)]<dist[g>head(e)]); … … 539 555 { 540 556 541 542 543 544 557 int heur0=(int)(H0*n); //time while running 'bound decrease' 558 int heur1=(int)(H1*n); //time while running 'highest label' 559 int heur=heur1; //starting time interval (#of relabels) 560 int numrelabel=0; 545 561 546 bool what_heur=1; 547 //It is 0 in case 'bound decrease' and 1 in case 'highest label' 548 549 bool end=false; 550 //Needed for 'bound decrease', true means no active nodes are above bound b. 551 552 int k=n2; //bound on the highest level under n containing a node 553 int b=k; //bound on the highest level under n of an active node 554 555 VecStack active(n); 556 557 NNMap left(*g, INVALID); 558 NNMap right(*g, INVALID); 559 VecNode level_list(n,INVALID); 560 //List of the nodes in level i<n, set to n. 561 562 NodeIt v; 563 for(g>first(v); g>valid(v); g>next(v)) level.set(v,n); 564 //setting each node to level n 565 566 switch ( fe ) { 567 case PREFLOW: 568 { 569 //counting the excess 570 NodeIt v; 571 for(g>first(v); g>valid(v); g>next(v)) { 572 Num exc=0; 573 574 InEdgeIt e; 575 for(g>first(e,v); g>valid(e); g>next(e)) exc+=(*flow)[e]; 576 OutEdgeIt f; 577 for(g>first(f,v); g>valid(f); g>next(f)) exc=(*flow)[f]; 578 579 excess.set(v,exc); 580 581 //putting the active nodes into the stack 582 int lev=level[v]; 583 if ( exc > 0 && lev < n && v != t ) active[lev].push(v); 562 bool what_heur=1; 563 //It is 0 in case 'bound decrease' and 1 in case 'highest label' 564 565 bool end=false; 566 //Needed for 'bound decrease', true means no active nodes are above bound b. 567 568 int k=n2; //bound on the highest level under n containing a node 569 int b=k; //bound on the highest level under n of an active node 570 571 VecStack active(n); 572 573 NNMap left(*g, INVALID); 574 NNMap right(*g, INVALID); 575 VecNode level_list(n,INVALID); 576 //List of the nodes in level i<n, set to n. 577 578 NodeIt v; 579 for(g>first(v); g>valid(v); g>next(v)) level.set(v,n); 580 //setting each node to level n 581 582 switch ( fe ) { 583 case PREFLOW: 584 { 585 //counting the excess 586 NodeIt v; 587 for(g>first(v); g>valid(v); g>next(v)) { 588 Num exc=0; 589 590 InEdgeIt e; 591 for(g>first(e,v); g>valid(e); g>next(e)) exc+=(*flow)[e]; 592 OutEdgeIt f; 593 for(g>first(f,v); g>valid(f); g>next(f)) exc=(*flow)[f]; 594 595 excess.set(v,exc); 596 597 //putting the active nodes into the stack 598 int lev=level[v]; 599 if ( exc > 0 && lev < n && v != t ) active[lev].push(v); 600 } 601 break; 602 } 603 case GEN_FLOW: 604 { 605 //Counting the excess of t 606 Num exc=0; 607 608 InEdgeIt e; 609 for(g>first(e,t); g>valid(e); g>next(e)) exc+=(*flow)[e]; 610 OutEdgeIt f; 611 for(g>first(f,t); g>valid(f); g>next(f)) exc=(*flow)[f]; 612 613 excess.set(t,exc); 614 615 break; 616 } 617 default: 618 break; 619 } 620 621 preflowPreproc( fe, active, level_list, left, right ); 622 //End of preprocessing 623 624 625 //Push/relabel on the highest level active nodes. 626 while ( true ) { 627 if ( b == 0 ) { 628 if ( !what_heur && !end && k > 0 ) { 629 b=k; 630 end=true; 631 } else break; 632 } 633 634 if ( active[b].empty() ) b; 635 else { 636 end=false; 637 Node w=active[b].top(); 638 active[b].pop(); 639 int newlevel=push(w,active); 640 if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list, 641 left, right, b, k, what_heur); 642 643 ++numrelabel; 644 if ( numrelabel >= heur ) { 645 numrelabel=0; 646 if ( what_heur ) { 647 what_heur=0; 648 heur=heur0; 649 end=false; 650 } else { 651 what_heur=1; 652 heur=heur1; 653 b=k; 584 654 } 585 break; 586 } 587 case GEN_FLOW: 588 { 589 //Counting the excess of t 590 Num exc=0; 591 592 InEdgeIt e; 593 for(g>first(e,t); g>valid(e); g>next(e)) exc+=(*flow)[e]; 594 OutEdgeIt f; 595 for(g>first(f,t); g>valid(f); g>next(f)) exc=(*flow)[f]; 596 597 excess.set(t,exc); 598 599 break; 600 } 601 default: 602 break; 603 } 604 605 preflowPreproc( fe, active, level_list, left, right ); 606 //End of preprocessing 607 608 609 //Push/relabel on the highest level active nodes. 610 while ( true ) { 611 if ( b == 0 ) { 612 if ( !what_heur && !end && k > 0 ) { 613 b=k; 614 end=true; 615 } else break; 616 } 617 618 if ( active[b].empty() ) b; 619 else { 620 end=false; 621 Node w=active[b].top(); 622 active[b].pop(); 623 int newlevel=push(w,active); 624 if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list, 625 left, right, b, k, what_heur); 626 627 ++numrelabel; 628 if ( numrelabel >= heur ) { 629 numrelabel=0; 630 if ( what_heur ) { 631 what_heur=0; 632 heur=heur0; 633 end=false; 634 } else { 635 what_heur=1; 636 heur=heur1; 637 b=k; 638 } 639 } 640 } 655 } 641 656 } 642 } 657 } 658 } 643 659 644 660 … … 648 664 { 649 665 650 651 652 653 654 655 656 657 658 659 660 661 662 666 int k=n2; //bound on the highest level under n containing a node 667 int b=k; //bound on the highest level under n of an active node 668 669 VecStack active(n); 670 level.set(s,0); 671 std::queue<Node> bfs_queue; 672 bfs_queue.push(s); 673 674 while (!bfs_queue.empty()) { 675 676 Node v=bfs_queue.front(); 677 bfs_queue.pop(); 678 int l=level[v]+1; 663 679 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 680 InEdgeIt e; 681 for(g>first(e,v); g>valid(e); g>next(e)) { 682 if ( (*capacity)[e] <= (*flow)[e] ) continue; 683 Node u=g>tail(e); 684 if ( level[u] >= n ) { 685 bfs_queue.push(u); 686 level.set(u, l); 687 if ( excess[u] > 0 ) active[l].push(u); 688 } 689 } 690 691 OutEdgeIt f; 692 for(g>first(f,v); g>valid(f); g>next(f)) { 693 if ( 0 >= (*flow)[f] ) continue; 694 Node u=g>head(f); 695 if ( level[u] >= n ) { 696 bfs_queue.push(u); 697 level.set(u, l); 698 if ( excess[u] > 0 ) active[l].push(u); 699 } 700 } 701 } 702 b=n2; 703 704 while ( true ) { 705 706 if ( b == 0 ) break; 707 708 if ( active[b].empty() ) b; 709 else { 710 Node w=active[b].top(); 711 active[b].pop(); 712 int newlevel=push(w,active); 713 714 //relabel 715 if ( excess[w] > 0 ) { 716 level.set(w,++newlevel); 717 active[newlevel].push(w); 718 b=newlevel; 719 } 720 } // if stack[b] is nonempty 721 } // while(true) 722 } 707 723 708 724 … … 711 727 bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath() 712 728 { 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 729 ResGW res_graph(*g, *capacity, *flow); 730 bool _augment=false; 731 732 //ReachedMap level(res_graph); 733 FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); 734 BfsIterator<ResGW, ReachedMap> bfs(res_graph, level); 735 bfs.pushAndSetReached(s); 736 737 typename ResGW::template NodeMap<ResGWEdge> pred(res_graph); 738 pred.set(s, INVALID); 739 740 typename ResGW::template NodeMap<Num> free(res_graph); 741 742 //searching for augmenting path 743 while ( !bfs.finished() ) { 744 ResGWOutEdgeIt e=bfs; 745 if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) { 746 Node v=res_graph.tail(e); 747 Node w=res_graph.head(e); 748 pred.set(w, e); 749 if (res_graph.valid(pred[v])) { 750 free.set(w, std::min(free[v], res_graph.resCap(e))); 751 } else { 752 free.set(w, res_graph.resCap(e)); 753 } 754 if (res_graph.head(e)==t) { _augment=true; break; } 755 } 756 757 ++bfs; 758 } //end of searching augmenting path 759 760 if (_augment) { 761 Node n=t; 762 Num augment_value=free[t]; 763 while (res_graph.valid(pred[n])) { 764 ResGWEdge e=pred[n]; 765 res_graph.augment(e, augment_value); 766 n=res_graph.tail(e); 767 } 768 } 769 770 return _augment; 771 } 756 772 757 773 … … 767 783 bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow() 768 784 { 769 typedef MutableGraph MG; 770 bool _augment=false; 771 772 ResGW res_graph(*g, *capacity, *flow); 773 774 //bfs for distances on the residual graph 775 //ReachedMap level(res_graph); 776 FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); 777 BfsIterator<ResGW, ReachedMap> bfs(res_graph, level); 778 bfs.pushAndSetReached(s); 779 typename ResGW::template NodeMap<int> 780 dist(res_graph); //filled up with 0's 781 782 //F will contain the physical copy of the residual graph 783 //with the set of edges which are on shortest paths 784 MG F; 785 typename ResGW::template NodeMap<typename MG::Node> 786 res_graph_to_F(res_graph); 787 { 788 typename ResGW::NodeIt n; 789 for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) { 790 res_graph_to_F.set(n, F.addNode()); 791 } 792 } 793 794 typename MG::Node sF=res_graph_to_F[s]; 795 typename MG::Node tF=res_graph_to_F[t]; 796 typename MG::template EdgeMap<ResGWEdge> original_edge(F); 797 typename MG::template EdgeMap<Num> residual_capacity(F); 798 799 while ( !bfs.finished() ) { 800 ResGWOutEdgeIt e=bfs; 801 if (res_graph.valid(e)) { 802 if (bfs.isBNodeNewlyReached()) { 803 dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1); 785 typedef MutableGraph MG; 786 bool _augment=false; 787 788 ResGW res_graph(*g, *capacity, *flow); 789 790 //bfs for distances on the residual graph 791 //ReachedMap level(res_graph); 792 FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); 793 BfsIterator<ResGW, ReachedMap> bfs(res_graph, level); 794 bfs.pushAndSetReached(s); 795 typename ResGW::template NodeMap<int> 796 dist(res_graph); //filled up with 0's 797 798 //F will contain the physical copy of the residual graph 799 //with the set of edges which are on shortest paths 800 MG F; 801 typename ResGW::template NodeMap<typename MG::Node> 802 res_graph_to_F(res_graph); 803 { 804 typename ResGW::NodeIt n; 805 for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) { 806 res_graph_to_F.set(n, F.addNode()); 807 } 808 } 809 810 typename MG::Node sF=res_graph_to_F[s]; 811 typename MG::Node tF=res_graph_to_F[t]; 812 typename MG::template EdgeMap<ResGWEdge> original_edge(F); 813 typename MG::template EdgeMap<Num> residual_capacity(F); 814 815 while ( !bfs.finished() ) { 816 ResGWOutEdgeIt e=bfs; 817 if (res_graph.valid(e)) { 818 if (bfs.isBNodeNewlyReached()) { 819 dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1); 820 typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]); 821 original_edge.update(); 822 original_edge.set(f, e); 823 residual_capacity.update(); 824 residual_capacity.set(f, res_graph.resCap(e)); 825 } else { 826 if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) { 804 827 typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]); 805 828 original_edge.update(); … … 807 830 residual_capacity.update(); 808 831 residual_capacity.set(f, res_graph.resCap(e)); 832 } 833 } 834 } 835 ++bfs; 836 } //computing distances from s in the residual graph 837 838 bool __augment=true; 839 840 while (__augment) { 841 __augment=false; 842 //computing blocking flow with dfs 843 DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F); 844 typename MG::template NodeMap<typename MG::Edge> pred(F); 845 pred.set(sF, INVALID); 846 //invalid iterators for sources 847 848 typename MG::template NodeMap<Num> free(F); 849 850 dfs.pushAndSetReached(sF); 851 while (!dfs.finished()) { 852 ++dfs; 853 if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) { 854 if (dfs.isBNodeNewlyReached()) { 855 typename MG::Node v=F.aNode(dfs); 856 typename MG::Node w=F.bNode(dfs); 857 pred.set(w, dfs); 858 if (F.valid(pred[v])) { 859 free.set(w, std::min(free[v], residual_capacity[dfs])); 860 } else { 861 free.set(w, residual_capacity[dfs]); 862 } 863 if (w==tF) { 864 __augment=true; 865 _augment=true; 866 break; 867 } 868 809 869 } else { 810 if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) { 811 typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]); 812 original_edge.update(); 813 original_edge.set(f, e); 814 residual_capacity.update(); 815 residual_capacity.set(f, res_graph.resCap(e)); 816 } 870 F.erase(/*typename MG::OutEdgeIt*/(dfs)); 817 871 } 818 } 819 ++bfs; 820 } //computing distances from s in the residual graph 821 822 bool __augment=true; 823 824 while (__augment) { 825 __augment=false; 826 //computing blocking flow with dfs 827 DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F); 828 typename MG::template NodeMap<typename MG::Edge> pred(F); 829 pred.set(sF, INVALID); 830 //invalid iterators for sources 831 832 typename MG::template NodeMap<Num> free(F); 833 834 dfs.pushAndSetReached(sF); 835 while (!dfs.finished()) { 836 ++dfs; 837 if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) { 838 if (dfs.isBNodeNewlyReached()) { 839 typename MG::Node v=F.aNode(dfs); 840 typename MG::Node w=F.bNode(dfs); 841 pred.set(w, dfs); 842 if (F.valid(pred[v])) { 843 free.set(w, std::min(free[v], residual_capacity[dfs])); 844 } else { 845 free.set(w, residual_capacity[dfs]); 846 } 847 if (w==tF) { 848 __augment=true; 849 _augment=true; 850 break; 851 } 852 853 } else { 854 F.erase(/*typename MG::OutEdgeIt*/(dfs)); 855 } 856 } 857 } 858 859 if (__augment) { 860 typename MG::Node n=tF; 861 Num augment_value=free[tF]; 862 while (F.valid(pred[n])) { 863 typename MG::Edge e=pred[n]; 864 res_graph.augment(original_edge[e], augment_value); 865 n=F.tail(e); 866 if (residual_capacity[e]==augment_value) 867 F.erase(e); 868 else 869 residual_capacity.set(e, residual_capacity[e]augment_value); 870 } 871 } 872 873 } 872 } 873 } 874 875 if (__augment) { 876 typename MG::Node n=tF; 877 Num augment_value=free[tF]; 878 while (F.valid(pred[n])) { 879 typename MG::Edge e=pred[n]; 880 res_graph.augment(original_edge[e], augment_value); 881 n=F.tail(e); 882 if (residual_capacity[e]==augment_value) 883 F.erase(e); 884 else 885 residual_capacity.set(e, residual_capacity[e]augment_value); 886 } 887 } 888 889 } 874 890 875 876 891 return _augment; 892 } 877 893 878 894 … … 884 900 bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2() 885 901 { 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 902 bool _augment=false; 903 904 ResGW res_graph(*g, *capacity, *flow); 905 906 //ReachedMap level(res_graph); 907 FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); 908 BfsIterator<ResGW, ReachedMap> bfs(res_graph, level); 909 910 bfs.pushAndSetReached(s); 911 DistanceMap<ResGW> dist(res_graph); 912 while ( !bfs.finished() ) { 913 ResGWOutEdgeIt e=bfs; 914 if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) { 915 dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1); 916 } 917 ++bfs; 918 } //computing distances from s in the residual graph 903 919 904 920 //Subgraph containing the edges on some shortest paths 905 906 907 908 909 910 911 912 913 914 915 916 921 ConstMap<typename ResGW::Node, bool> true_map(true); 922 typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>, 923 DistanceMap<ResGW> > FilterResGW; 924 FilterResGW filter_res_graph(res_graph, true_map, dist); 925 926 //Subgraph, which is able to delete edges which are already 927 //met by the dfs 928 typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt> 929 first_out_edges(filter_res_graph); 930 typename FilterResGW::NodeIt v; 931 for(filter_res_graph.first(v); filter_res_graph.valid(v); 932 filter_res_graph.next(v)) 917 933 { 918 934 typename FilterResGW::OutEdgeIt e; … … 920 936 first_out_edges.set(v, e); 921 937 } 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 typename ErasingResGW::Node(946 947 948 949 950 );951 952 953 954 938 typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW:: 939 template NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW; 940 ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges); 941 942 bool __augment=true; 943 944 while (__augment) { 945 946 __augment=false; 947 //computing blocking flow with dfs 948 DfsIterator< ErasingResGW, 949 typename ErasingResGW::template NodeMap<bool> > 950 dfs(erasing_res_graph); 951 typename ErasingResGW:: 952 template NodeMap<typename ErasingResGW::OutEdgeIt> 953 pred(erasing_res_graph); 954 pred.set(s, INVALID); 955 //invalid iterators for sources 956 957 typename ErasingResGW::template NodeMap<Num> 958 free1(erasing_res_graph); 959 960 dfs.pushAndSetReached( 961 typename ErasingResGW::Node( 962 typename FilterResGW::Node( 963 typename ResGW::Node(s) 964 ) 965 ) 966 ); 967 while (!dfs.finished()) { 968 ++dfs; 969 if (erasing_res_graph.valid( 970 typename ErasingResGW::OutEdgeIt(dfs))) 955 971 { 956 972 if (dfs.isBNodeNewlyReached()) { … … 962 978 if (erasing_res_graph.valid(pred[v])) { 963 979 free1.set(w, std::min(free1[v], res_graph.resCap( 964 980 typename ErasingResGW::OutEdgeIt(dfs)))); 965 981 } else { 966 982 free1.set(w, res_graph.resCap( 967 typename ErasingResGW::OutEdgeIt(dfs)));983 typename ErasingResGW::OutEdgeIt(dfs))); 968 984 } 969 985 … … 977 993 } 978 994 } 979 980 981 982 983 // typename ResGW::NodeMap<Num> a(res_graph);984 // typename ResGW::Node b;985 // Num j=a[b];986 // typename FilterResGW::NodeMap<Num> a1(filter_res_graph);987 // typename FilterResGW::Node b1;988 // Num j1=a1[b1];989 // typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);990 // typename ErasingResGW::Node b2;991 // Num j2=a2[b2];992 993 994 995 996 997 998 999 } 1000 } 1001 1002 995 } 996 997 if (__augment) { 998 typename ErasingResGW::Node n=typename FilterResGW::Node(typename ResGW::Node(t)); 999 // typename ResGW::NodeMap<Num> a(res_graph); 1000 // typename ResGW::Node b; 1001 // Num j=a[b]; 1002 // typename FilterResGW::NodeMap<Num> a1(filter_res_graph); 1003 // typename FilterResGW::Node b1; 1004 // Num j1=a1[b1]; 1005 // typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph); 1006 // typename ErasingResGW::Node b2; 1007 // Num j2=a2[b2]; 1008 Num augment_value=free1[n]; 1009 while (erasing_res_graph.valid(pred[n])) { 1010 typename ErasingResGW::OutEdgeIt e=pred[n]; 1011 res_graph.augment(e, augment_value); 1012 n=erasing_res_graph.tail(e); 1013 if (res_graph.resCap(e)==0) 1014 erasing_res_graph.erase(e); 1015 } 1016 } 1017 1018 } //while (__augment) 1003 1019 1004 1005 1020 return _augment; 1021 } 1006 1022 1007 1023
Note: See TracChangeset
for help on using the changeset viewer.