00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #ifndef LEMON_PREFLOW_H
00018 #define LEMON_PREFLOW_H
00019
00020 #include <vector>
00021 #include <queue>
00022
00023 #include <lemon/invalid.h>
00024 #include <lemon/maps.h>
00025 #include <lemon/graph_utils.h>
00026
00030
00031 namespace lemon {
00032
00035
00037
00062 template <typename Graph, typename Num,
00063 typename CapMap=typename Graph::template EdgeMap<Num>,
00064 typename FlowMap=typename Graph::template EdgeMap<Num> >
00065 class Preflow {
00066 protected:
00067 typedef typename Graph::Node Node;
00068 typedef typename Graph::NodeIt NodeIt;
00069 typedef typename Graph::EdgeIt EdgeIt;
00070 typedef typename Graph::OutEdgeIt OutEdgeIt;
00071 typedef typename Graph::InEdgeIt InEdgeIt;
00072
00073 typedef typename Graph::template NodeMap<Node> NNMap;
00074 typedef typename std::vector<Node> VecNode;
00075
00076 const Graph* g;
00077 Node s;
00078 Node t;
00079 const CapMap* capacity;
00080 FlowMap* flow;
00081 int n;
00082
00083 typename Graph::template NodeMap<int> level;
00084 typename Graph::template NodeMap<Num> excess;
00085
00086
00087 static const int H0=20;
00088 static const int H1=1;
00089
00090 public:
00091
00093
00105 enum FlowEnum{
00106 NO_FLOW,
00107 ZERO_FLOW,
00108 GEN_FLOW,
00109 PRE_FLOW
00110 };
00111
00113
00119 enum StatusEnum {
00120 AFTER_NOTHING,
00121 AFTER_PREFLOW_PHASE_1,
00122 AFTER_PREFLOW_PHASE_2
00123 };
00124
00125 protected:
00126 FlowEnum flow_prop;
00127 StatusEnum status;
00128
00129 public:
00131
00141 Preflow(const Graph& _G, Node _s, Node _t,
00142 const CapMap& _capacity, FlowMap& _flow) :
00143 g(&_G), s(_s), t(_t), capacity(&_capacity),
00144 flow(&_flow), n(countNodes(_G)), level(_G), excess(_G,0),
00145 flow_prop(NO_FLOW), status(AFTER_NOTHING) { }
00146
00147
00148
00150
00153 void run() {
00154 phase1(flow_prop);
00155 phase2();
00156 }
00157
00159
00168 void run(FlowEnum fp) {
00169 flow_prop=fp;
00170 run();
00171 }
00172
00174
00188 void phase1(FlowEnum fp)
00189 {
00190 flow_prop=fp;
00191 phase1();
00192 }
00193
00194
00196
00205 void phase1()
00206 {
00207 int heur0=(int)(H0*n);
00208 int heur1=(int)(H1*n);
00209 int heur=heur1;
00210 int numrelabel=0;
00211
00212 bool what_heur=1;
00213
00214
00215 bool end=false;
00216
00217
00218
00219 int k=n-2;
00220 int b=k;
00221
00222 VecNode first(n, INVALID);
00223 NNMap next(*g, INVALID);
00224
00225 NNMap left(*g, INVALID);
00226 NNMap right(*g, INVALID);
00227 VecNode level_list(n,INVALID);
00228
00229
00230 preflowPreproc(first, next, level_list, left, right);
00231
00232
00233 while ( true ) {
00234 if ( b == 0 ) {
00235 if ( !what_heur && !end && k > 0 ) {
00236 b=k;
00237 end=true;
00238 } else break;
00239 }
00240
00241 if ( first[b]==INVALID ) --b;
00242 else {
00243 end=false;
00244 Node w=first[b];
00245 first[b]=next[w];
00246 int newlevel=push(w, next, first);
00247 if ( excess[w] > 0 ) relabel(w, newlevel, first, next, level_list,
00248 left, right, b, k, what_heur);
00249
00250 ++numrelabel;
00251 if ( numrelabel >= heur ) {
00252 numrelabel=0;
00253 if ( what_heur ) {
00254 what_heur=0;
00255 heur=heur0;
00256 end=false;
00257 } else {
00258 what_heur=1;
00259 heur=heur1;
00260 b=k;
00261 }
00262 }
00263 }
00264 }
00265 flow_prop=PRE_FLOW;
00266 status=AFTER_PREFLOW_PHASE_1;
00267 }
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00279
00287 void phase2()
00288 {
00289
00290 int k=n-2;
00291 int b=k;
00292
00293
00294 VecNode first(n, INVALID);
00295 NNMap next(*g, INVALID);
00296 level.set(s,0);
00297 std::queue<Node> bfs_queue;
00298 bfs_queue.push(s);
00299
00300 while ( !bfs_queue.empty() ) {
00301
00302 Node v=bfs_queue.front();
00303 bfs_queue.pop();
00304 int l=level[v]+1;
00305
00306 for(InEdgeIt e(*g,v); e!=INVALID; ++e) {
00307 if ( (*capacity)[e] <= (*flow)[e] ) continue;
00308 Node u=g->source(e);
00309 if ( level[u] >= n ) {
00310 bfs_queue.push(u);
00311 level.set(u, l);
00312 if ( excess[u] > 0 ) {
00313 next.set(u,first[l]);
00314 first[l]=u;
00315 }
00316 }
00317 }
00318
00319 for(OutEdgeIt e(*g,v); e!=INVALID; ++e) {
00320 if ( 0 >= (*flow)[e] ) continue;
00321 Node u=g->target(e);
00322 if ( level[u] >= n ) {
00323 bfs_queue.push(u);
00324 level.set(u, l);
00325 if ( excess[u] > 0 ) {
00326 next.set(u,first[l]);
00327 first[l]=u;
00328 }
00329 }
00330 }
00331 }
00332 b=n-2;
00333
00334 while ( true ) {
00335
00336 if ( b == 0 ) break;
00337 if ( first[b]==INVALID ) --b;
00338 else {
00339 Node w=first[b];
00340 first[b]=next[w];
00341 int newlevel=push(w,next, first);
00342
00343
00344 if ( excess[w] > 0 ) {
00345 level.set(w,++newlevel);
00346 next.set(w,first[newlevel]);
00347 first[newlevel]=w;
00348 b=newlevel;
00349 }
00350 }
00351 }
00352 flow_prop=GEN_FLOW;
00353 status=AFTER_PREFLOW_PHASE_2;
00354 }
00355
00357
00361 Num flowValue() const {
00362 return excess[t];
00363 }
00364
00365
00367
00374 template<typename _CutMap>
00375 void minCut(_CutMap& M) const {
00376 switch ( status ) {
00377 case AFTER_PREFLOW_PHASE_1:
00378 for(NodeIt v(*g); v!=INVALID; ++v) {
00379 if (level[v] < n) {
00380 M.set(v, false);
00381 } else {
00382 M.set(v, true);
00383 }
00384 }
00385 break;
00386 case AFTER_PREFLOW_PHASE_2:
00387 minMinCut(M);
00388 break;
00389 case AFTER_NOTHING:
00390 break;
00391 }
00392 }
00393
00395
00401 template<typename _CutMap>
00402 void minMinCut(_CutMap& M) const {
00403
00404 std::queue<Node> queue;
00405 M.set(s,true);
00406 queue.push(s);
00407
00408 while (!queue.empty()) {
00409 Node w=queue.front();
00410 queue.pop();
00411
00412 for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
00413 Node v=g->target(e);
00414 if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
00415 queue.push(v);
00416 M.set(v, true);
00417 }
00418 }
00419
00420 for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
00421 Node v=g->source(e);
00422 if (!M[v] && (*flow)[e] > 0 ) {
00423 queue.push(v);
00424 M.set(v, true);
00425 }
00426 }
00427 }
00428 }
00429
00431
00436 template<typename _CutMap>
00437 void maxMinCut(_CutMap& M) const {
00438
00439 for(NodeIt v(*g) ; v!=INVALID; ++v) M.set(v, true);
00440
00441 std::queue<Node> queue;
00442
00443 M.set(t,false);
00444 queue.push(t);
00445
00446 while (!queue.empty()) {
00447 Node w=queue.front();
00448 queue.pop();
00449
00450 for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
00451 Node v=g->source(e);
00452 if (M[v] && (*flow)[e] < (*capacity)[e] ) {
00453 queue.push(v);
00454 M.set(v, false);
00455 }
00456 }
00457
00458 for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
00459 Node v=g->target(e);
00460 if (M[v] && (*flow)[e] > 0 ) {
00461 queue.push(v);
00462 M.set(v, false);
00463 }
00464 }
00465 }
00466 }
00467
00469
00472 void setSource(Node _s) {
00473 s=_s;
00474 if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
00475 status=AFTER_NOTHING;
00476 }
00477
00479
00482 void setTarget(Node _t) {
00483 t=_t;
00484 if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
00485 status=AFTER_NOTHING;
00486 }
00487
00489
00492 void setCap(const CapMap& _cap) {
00493 capacity=&_cap;
00494 status=AFTER_NOTHING;
00495 }
00496
00498
00501 void setFlow(FlowMap& _flow) {
00502 flow=&_flow;
00503 flow_prop=NO_FLOW;
00504 status=AFTER_NOTHING;
00505 }
00506
00507
00508 private:
00509
00510 int push(Node w, NNMap& next, VecNode& first) {
00511
00512 int lev=level[w];
00513 Num exc=excess[w];
00514 int newlevel=n;
00515
00516 for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
00517 if ( (*flow)[e] >= (*capacity)[e] ) continue;
00518 Node v=g->target(e);
00519
00520 if( lev > level[v] ) {
00521
00522 if ( excess[v]<=0 && v!=t && v!=s ) {
00523 next.set(v,first[level[v]]);
00524 first[level[v]]=v;
00525 }
00526
00527 Num cap=(*capacity)[e];
00528 Num flo=(*flow)[e];
00529 Num remcap=cap-flo;
00530
00531 if ( remcap >= exc ) {
00532
00533 flow->set(e, flo+exc);
00534 excess.set(v, excess[v]+exc);
00535 exc=0;
00536 break;
00537
00538 } else {
00539 flow->set(e, cap);
00540 excess.set(v, excess[v]+remcap);
00541 exc-=remcap;
00542 }
00543 } else if ( newlevel > level[v] ) newlevel = level[v];
00544 }
00545
00546 if ( exc > 0 ) {
00547 for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
00548
00549 if( (*flow)[e] <= 0 ) continue;
00550 Node v=g->source(e);
00551
00552 if( lev > level[v] ) {
00553
00554 if ( excess[v]<=0 && v!=t && v!=s ) {
00555 next.set(v,first[level[v]]);
00556 first[level[v]]=v;
00557 }
00558
00559 Num flo=(*flow)[e];
00560
00561 if ( flo >= exc ) {
00562
00563 flow->set(e, flo-exc);
00564 excess.set(v, excess[v]+exc);
00565 exc=0;
00566 break;
00567 } else {
00568
00569 excess.set(v, excess[v]+flo);
00570 exc-=flo;
00571 flow->set(e,0);
00572 }
00573 } else if ( newlevel > level[v] ) newlevel = level[v];
00574 }
00575
00576 }
00577
00578 excess.set(w, exc);
00579
00580 return newlevel;
00581 }
00582
00583
00584
00585 void preflowPreproc(VecNode& first, NNMap& next,
00586 VecNode& level_list, NNMap& left, NNMap& right)
00587 {
00588 for(NodeIt v(*g); v!=INVALID; ++v) level.set(v,n);
00589 std::queue<Node> bfs_queue;
00590
00591 if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
00592
00593
00594 level.set(t,0);
00595 bfs_queue.push(t);
00596
00597 while ( !bfs_queue.empty() ) {
00598
00599 Node v=bfs_queue.front();
00600 bfs_queue.pop();
00601 int l=level[v]+1;
00602
00603 for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) {
00604 if ( (*capacity)[e] <= (*flow)[e] ) continue;
00605 Node w=g->source(e);
00606 if ( level[w] == n && w != s ) {
00607 bfs_queue.push(w);
00608 Node z=level_list[l];
00609 if ( z!=INVALID ) left.set(z,w);
00610 right.set(w,z);
00611 level_list[l]=w;
00612 level.set(w, l);
00613 }
00614 }
00615
00616 for(OutEdgeIt e(*g,v) ; e!=INVALID; ++e) {
00617 if ( 0 >= (*flow)[e] ) continue;
00618 Node w=g->target(e);
00619 if ( level[w] == n && w != s ) {
00620 bfs_queue.push(w);
00621 Node z=level_list[l];
00622 if ( z!=INVALID ) left.set(z,w);
00623 right.set(w,z);
00624 level_list[l]=w;
00625 level.set(w, l);
00626 }
00627 }
00628 }
00629 }
00630
00631
00632 switch (flow_prop) {
00633 case NO_FLOW:
00634 for(EdgeIt e(*g); e!=INVALID; ++e) flow->set(e,0);
00635 case ZERO_FLOW:
00636 for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0);
00637
00638
00639 level.set(t,0);
00640 bfs_queue.push(t);
00641
00642 while ( !bfs_queue.empty() ) {
00643
00644 Node v=bfs_queue.front();
00645 bfs_queue.pop();
00646 int l=level[v]+1;
00647
00648 for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) {
00649 Node w=g->source(e);
00650 if ( level[w] == n && w != s ) {
00651 bfs_queue.push(w);
00652 Node z=level_list[l];
00653 if ( z!=INVALID ) left.set(z,w);
00654 right.set(w,z);
00655 level_list[l]=w;
00656 level.set(w, l);
00657 }
00658 }
00659 }
00660
00661
00662 for(OutEdgeIt e(*g,s) ; e!=INVALID; ++e) {
00663 Num c=(*capacity)[e];
00664 if ( c <= 0 ) continue;
00665 Node w=g->target(e);
00666 if ( level[w] < n ) {
00667 if ( excess[w] <= 0 && w!=t ) {
00668 next.set(w,first[level[w]]);
00669 first[level[w]]=w;
00670 }
00671 flow->set(e, c);
00672 excess.set(w, excess[w]+c);
00673 }
00674 }
00675 break;
00676
00677 case GEN_FLOW:
00678 for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0);
00679 {
00680 Num exc=0;
00681 for(InEdgeIt e(*g,t) ; e!=INVALID; ++e) exc+=(*flow)[e];
00682 for(OutEdgeIt e(*g,t) ; e!=INVALID; ++e) exc-=(*flow)[e];
00683 excess.set(t,exc);
00684 }
00685
00686
00687 for(OutEdgeIt e(*g,s); e!=INVALID; ++e) {
00688 Num rem=(*capacity)[e]-(*flow)[e];
00689 if ( rem <= 0 ) continue;
00690 Node w=g->target(e);
00691 if ( level[w] < n ) {
00692 if ( excess[w] <= 0 && w!=t ) {
00693 next.set(w,first[level[w]]);
00694 first[level[w]]=w;
00695 }
00696 flow->set(e, (*capacity)[e]);
00697 excess.set(w, excess[w]+rem);
00698 }
00699 }
00700
00701 for(InEdgeIt e(*g,s); e!=INVALID; ++e) {
00702 if ( (*flow)[e] <= 0 ) continue;
00703 Node w=g->source(e);
00704 if ( level[w] < n ) {
00705 if ( excess[w] <= 0 && w!=t ) {
00706 next.set(w,first[level[w]]);
00707 first[level[w]]=w;
00708 }
00709 excess.set(w, excess[w]+(*flow)[e]);
00710 flow->set(e, 0);
00711 }
00712 }
00713 break;
00714
00715 case PRE_FLOW:
00716
00717 for(OutEdgeIt e(*g,s) ; e!=INVALID; ++e) {
00718 Num rem=(*capacity)[e]-(*flow)[e];
00719 if ( rem <= 0 ) continue;
00720 Node w=g->target(e);
00721 if ( level[w] < n ) flow->set(e, (*capacity)[e]);
00722 }
00723
00724 for(InEdgeIt e(*g,s) ; e!=INVALID; ++e) {
00725 if ( (*flow)[e] <= 0 ) continue;
00726 Node w=g->source(e);
00727 if ( level[w] < n ) flow->set(e, 0);
00728 }
00729
00730
00731 for(NodeIt w(*g); w!=INVALID; ++w) {
00732 Num exc=0;
00733 for(InEdgeIt e(*g,w); e!=INVALID; ++e) exc+=(*flow)[e];
00734 for(OutEdgeIt e(*g,w); e!=INVALID; ++e) exc-=(*flow)[e];
00735 excess.set(w,exc);
00736
00737
00738 int lev=level[w];
00739 if ( exc > 0 && lev < n && Node(w) != t ) {
00740 next.set(w,first[lev]);
00741 first[lev]=w;
00742 }
00743 }
00744 break;
00745 }
00746 }
00747
00748
00749 void relabel(Node w, int newlevel, VecNode& first, NNMap& next,
00750 VecNode& level_list, NNMap& left,
00751 NNMap& right, int& b, int& k, bool what_heur )
00752 {
00753
00754 int lev=level[w];
00755
00756 Node right_n=right[w];
00757 Node left_n=left[w];
00758
00759
00760 if ( right_n!=INVALID ) {
00761 if ( left_n!=INVALID ) {
00762 right.set(left_n, right_n);
00763 left.set(right_n, left_n);
00764 } else {
00765 level_list[lev]=right_n;
00766 left.set(right_n, INVALID);
00767 }
00768 } else {
00769 if ( left_n!=INVALID ) {
00770 right.set(left_n, INVALID);
00771 } else {
00772 level_list[lev]=INVALID;
00773 }
00774 }
00775
00776
00777 if ( level_list[lev]==INVALID ) {
00778
00779
00780 for (int i=lev; i!=k ; ) {
00781 Node v=level_list[++i];
00782 while ( v!=INVALID ) {
00783 level.set(v,n);
00784 v=right[v];
00785 }
00786 level_list[i]=INVALID;
00787 if ( !what_heur ) first[i]=INVALID;
00788 }
00789
00790 level.set(w,n);
00791 b=lev-1;
00792 k=b;
00793
00794
00795 } else {
00796
00797 if ( newlevel == n ) level.set(w,n);
00798 else {
00799 level.set(w,++newlevel);
00800 next.set(w,first[newlevel]);
00801 first[newlevel]=w;
00802 if ( what_heur ) b=newlevel;
00803 if ( k < newlevel ) ++k;
00804 Node z=level_list[newlevel];
00805 if ( z!=INVALID ) left.set(z,w);
00806 right.set(w,z);
00807 left.set(w,INVALID);
00808 level_list[newlevel]=w;
00809 }
00810 }
00811 }
00812
00813 };
00814 }
00815
00816 #endif
00817
00818
00819
00820