00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #ifndef LEMON_PREFLOW_H
00020 #define LEMON_PREFLOW_H
00021
00022 #include <vector>
00023 #include <queue>
00024
00025 #include <lemon/error.h>
00026 #include <lemon/invalid.h>
00027 #include <lemon/tolerance.h>
00028 #include <lemon/maps.h>
00029 #include <lemon/graph_utils.h>
00030
00034
00035 namespace lemon {
00036
00065 template <typename Graph, typename Num,
00066 typename CapacityMap=typename Graph::template EdgeMap<Num>,
00067 typename FlowMap=typename Graph::template EdgeMap<Num>,
00068 typename TOL=Tolerance<Num> >
00069 class Preflow {
00070 protected:
00071 typedef typename Graph::Node Node;
00072 typedef typename Graph::NodeIt NodeIt;
00073 typedef typename Graph::EdgeIt EdgeIt;
00074 typedef typename Graph::OutEdgeIt OutEdgeIt;
00075 typedef typename Graph::InEdgeIt InEdgeIt;
00076
00077 typedef typename Graph::template NodeMap<Node> NNMap;
00078 typedef typename std::vector<Node> VecNode;
00079
00080 const Graph* _g;
00081 Node _source;
00082 Node _target;
00083 const CapacityMap* _capacity;
00084 FlowMap* _flow;
00085
00086 TOL surely;
00087
00088 int _node_num;
00089
00090 typename Graph::template NodeMap<int> level;
00091 typename Graph::template NodeMap<Num> excess;
00092
00093
00094 static const int H0=20;
00095 static const int H1=1;
00096
00097 public:
00098
00100
00102 class InvalidArgument : public lemon::LogicError {
00103 public:
00104 virtual const char* exceptionName() const {
00105 return "lemon::Preflow::InvalidArgument";
00106 }
00107 };
00108
00109
00111
00114 enum FlowEnum{
00118 NO_FLOW,
00120 ZERO_FLOW,
00124 GEN_FLOW,
00127 PRE_FLOW
00128 };
00129
00131
00134 enum StatusEnum {
00137 AFTER_NOTHING,
00139 AFTER_PREFLOW_PHASE_1,
00141 AFTER_PREFLOW_PHASE_2
00142 };
00143
00144 protected:
00145 FlowEnum flow_prop;
00146 StatusEnum status;
00147
00148 public:
00150
00161 Preflow(const Graph& _gr, Node _s, Node _t,
00162 const CapacityMap& _cap, FlowMap& _f,
00163 const TOL &tol=TOL()) :
00164 _g(&_gr), _source(_s), _target(_t), _capacity(&_cap),
00165 _flow(&_f), surely(tol),
00166 _node_num(countNodes(_gr)), level(_gr), excess(_gr,0),
00167 flow_prop(NO_FLOW), status(AFTER_NOTHING) {
00168 if ( _source==_target )
00169 throw InvalidArgument();
00170 }
00171
00173
00176 TOL &tolerance() { return surely; }
00177
00179
00182 void run() {
00183 phase1(flow_prop);
00184 phase2();
00185 }
00186
00188
00197 void run(FlowEnum fp) {
00198 flow_prop=fp;
00199 run();
00200 }
00201
00203
00217 void phase1(FlowEnum fp)
00218 {
00219 flow_prop=fp;
00220 phase1();
00221 }
00222
00223
00225
00234 void phase1()
00235 {
00236 int heur0=(int)(H0*_node_num);
00237 int heur1=(int)(H1*_node_num);
00238 int heur=heur1;
00239 int numrelabel=0;
00240
00241 bool what_heur=1;
00242
00243
00244 bool end=false;
00245
00246
00247
00248 int k=_node_num-2;
00249 int b=k;
00250
00251 VecNode first(_node_num, INVALID);
00252 NNMap next(*_g, INVALID);
00253
00254 NNMap left(*_g, INVALID);
00255 NNMap right(*_g, INVALID);
00256 VecNode level_list(_node_num,INVALID);
00257
00258
00259 preflowPreproc(first, next, level_list, left, right);
00260
00261
00262 while ( true ) {
00263 if ( b == 0 ) {
00264 if ( !what_heur && !end && k > 0 ) {
00265 b=k;
00266 end=true;
00267 } else break;
00268 }
00269
00270 if ( first[b]==INVALID ) --b;
00271 else {
00272 end=false;
00273 Node w=first[b];
00274 first[b]=next[w];
00275 int newlevel=push(w, next, first);
00276 if ( excess[w] > 0 ) relabel(w, newlevel, first, next, level_list,
00277 left, right, b, k, what_heur);
00278
00279 ++numrelabel;
00280 if ( numrelabel >= heur ) {
00281 numrelabel=0;
00282 if ( what_heur ) {
00283 what_heur=0;
00284 heur=heur0;
00285 end=false;
00286 } else {
00287 what_heur=1;
00288 heur=heur1;
00289 b=k;
00290 }
00291 }
00292 }
00293 }
00294 flow_prop=PRE_FLOW;
00295 status=AFTER_PREFLOW_PHASE_1;
00296 }
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00309
00318 void phase2()
00319 {
00320
00321 int k=_node_num-2;
00322 int b=k;
00323
00324
00325 VecNode first(_node_num, INVALID);
00326 NNMap next(*_g, INVALID);
00327 level.set(_source,0);
00328 std::queue<Node> bfs_queue;
00329 bfs_queue.push(_source);
00330
00331 while ( !bfs_queue.empty() ) {
00332
00333 Node v=bfs_queue.front();
00334 bfs_queue.pop();
00335 int l=level[v]+1;
00336
00337 for(InEdgeIt e(*_g,v); e!=INVALID; ++e) {
00338 if ( (*_capacity)[e] <= (*_flow)[e] ) continue;
00339 Node u=_g->source(e);
00340 if ( level[u] >= _node_num ) {
00341 bfs_queue.push(u);
00342 level.set(u, l);
00343 if ( excess[u] > 0 ) {
00344 next.set(u,first[l]);
00345 first[l]=u;
00346 }
00347 }
00348 }
00349
00350 for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) {
00351 if ( 0 >= (*_flow)[e] ) continue;
00352 Node u=_g->target(e);
00353 if ( level[u] >= _node_num ) {
00354 bfs_queue.push(u);
00355 level.set(u, l);
00356 if ( excess[u] > 0 ) {
00357 next.set(u,first[l]);
00358 first[l]=u;
00359 }
00360 }
00361 }
00362 }
00363 b=_node_num-2;
00364
00365 while ( true ) {
00366
00367 if ( b == 0 ) break;
00368 if ( first[b]==INVALID ) --b;
00369 else {
00370 Node w=first[b];
00371 first[b]=next[w];
00372 int newlevel=push(w,next, first);
00373
00374
00375 if ( excess[w] > 0 ) {
00376 level.set(w,++newlevel);
00377 next.set(w,first[newlevel]);
00378 first[newlevel]=w;
00379 b=newlevel;
00380 }
00381 }
00382 }
00383 flow_prop=GEN_FLOW;
00384 status=AFTER_PREFLOW_PHASE_2;
00385 }
00386
00388
00392 Num flowValue() const {
00393 return excess[_target];
00394 }
00395
00396
00398
00405 template<typename _CutMap>
00406 void minCut(_CutMap& M) const {
00407 switch ( status ) {
00408 case AFTER_PREFLOW_PHASE_1:
00409 for(NodeIt v(*_g); v!=INVALID; ++v) {
00410 if (level[v] < _node_num) {
00411 M.set(v, false);
00412 } else {
00413 M.set(v, true);
00414 }
00415 }
00416 break;
00417 case AFTER_PREFLOW_PHASE_2:
00418 minMinCut(M);
00419 break;
00420 case AFTER_NOTHING:
00421 break;
00422 }
00423 }
00424
00426
00432 template<typename _CutMap>
00433 void minMinCut(_CutMap& M) const {
00434
00435 std::queue<Node> queue;
00436 M.set(_source,true);
00437 queue.push(_source);
00438
00439 while (!queue.empty()) {
00440 Node w=queue.front();
00441 queue.pop();
00442
00443 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
00444 Node v=_g->target(e);
00445 if (!M[v] && (*_flow)[e] < (*_capacity)[e] ) {
00446 queue.push(v);
00447 M.set(v, true);
00448 }
00449 }
00450
00451 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
00452 Node v=_g->source(e);
00453 if (!M[v] && (*_flow)[e] > 0 ) {
00454 queue.push(v);
00455 M.set(v, true);
00456 }
00457 }
00458 }
00459 }
00460
00462
00467 template<typename _CutMap>
00468 void maxMinCut(_CutMap& M) const {
00469
00470 for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true);
00471
00472 std::queue<Node> queue;
00473
00474 M.set(_target,false);
00475 queue.push(_target);
00476
00477 while (!queue.empty()) {
00478 Node w=queue.front();
00479 queue.pop();
00480
00481 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
00482 Node v=_g->source(e);
00483 if (M[v] && (*_flow)[e] < (*_capacity)[e] ) {
00484 queue.push(v);
00485 M.set(v, false);
00486 }
00487 }
00488
00489 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
00490 Node v=_g->target(e);
00491 if (M[v] && (*_flow)[e] > 0 ) {
00492 queue.push(v);
00493 M.set(v, false);
00494 }
00495 }
00496 }
00497 }
00498
00500
00503 void source(Node _s) {
00504 _source=_s;
00505 if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
00506 status=AFTER_NOTHING;
00507 }
00508
00510
00513 Node source() const {
00514 return _source;
00515 }
00516
00518
00521 void target(Node _t) {
00522 _target=_t;
00523 if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
00524 status=AFTER_NOTHING;
00525 }
00526
00528
00531 Node target() const {
00532 return _target;
00533 }
00534
00536
00539 void capacityMap(const CapacityMap& _cap) {
00540 _capacity=&_cap;
00541 status=AFTER_NOTHING;
00542 }
00544
00547 const CapacityMap &capacityMap() const {
00548 return *_capacity;
00549 }
00550
00552
00555 void flowMap(FlowMap& _f) {
00556 _flow=&_f;
00557 flow_prop=NO_FLOW;
00558 status=AFTER_NOTHING;
00559 }
00560
00562
00565 const FlowMap &flowMap() const {
00566 return *_flow;
00567 }
00568
00569 private:
00570
00571 int push(Node w, NNMap& next, VecNode& first) {
00572
00573 int lev=level[w];
00574 Num exc=excess[w];
00575 int newlevel=_node_num;
00576
00577 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
00578 if ( (*_flow)[e] >= (*_capacity)[e] ) continue;
00579 Node v=_g->target(e);
00580
00581 if( lev > level[v] ) {
00582
00583 if ( excess[v]<=0 && v!=_target && v!=_source ) {
00584 next.set(v,first[level[v]]);
00585 first[level[v]]=v;
00586 }
00587
00588 Num cap=(*_capacity)[e];
00589 Num flo=(*_flow)[e];
00590 Num remcap=cap-flo;
00591
00592 if ( remcap >= exc ) {
00593
00594 _flow->set(e, flo+exc);
00595 excess.set(v, excess[v]+exc);
00596 exc=0;
00597 break;
00598
00599 } else {
00600 _flow->set(e, cap);
00601 excess.set(v, excess[v]+remcap);
00602 exc-=remcap;
00603 }
00604 } else if ( newlevel > level[v] ) newlevel = level[v];
00605 }
00606
00607 if ( exc > 0 ) {
00608 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
00609
00610 if( (*_flow)[e] <= 0 ) continue;
00611 Node v=_g->source(e);
00612
00613 if( lev > level[v] ) {
00614
00615 if ( excess[v]<=0 && v!=_target && v!=_source ) {
00616 next.set(v,first[level[v]]);
00617 first[level[v]]=v;
00618 }
00619
00620 Num flo=(*_flow)[e];
00621
00622 if ( flo >= exc ) {
00623
00624 _flow->set(e, flo-exc);
00625 excess.set(v, excess[v]+exc);
00626 exc=0;
00627 break;
00628 } else {
00629
00630 excess.set(v, excess[v]+flo);
00631 exc-=flo;
00632 _flow->set(e,0);
00633 }
00634 } else if ( newlevel > level[v] ) newlevel = level[v];
00635 }
00636
00637 }
00638
00639 excess.set(w, exc);
00640
00641 return newlevel;
00642 }
00643
00644
00645
00646 void preflowPreproc(VecNode& first, NNMap& next,
00647 VecNode& level_list, NNMap& left, NNMap& right)
00648 {
00649 for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num);
00650 std::queue<Node> bfs_queue;
00651
00652 if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
00653
00654
00655 level.set(_target,0);
00656 bfs_queue.push(_target);
00657
00658 while ( !bfs_queue.empty() ) {
00659
00660 Node v=bfs_queue.front();
00661 bfs_queue.pop();
00662 int l=level[v]+1;
00663
00664 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
00665 if ( (*_capacity)[e] <= (*_flow)[e] ) continue;
00666 Node w=_g->source(e);
00667 if ( level[w] == _node_num && w != _source ) {
00668 bfs_queue.push(w);
00669 Node z=level_list[l];
00670 if ( z!=INVALID ) left.set(z,w);
00671 right.set(w,z);
00672 level_list[l]=w;
00673 level.set(w, l);
00674 }
00675 }
00676
00677 for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
00678 if ( 0 >= (*_flow)[e] ) continue;
00679 Node w=_g->target(e);
00680 if ( level[w] == _node_num && w != _source ) {
00681 bfs_queue.push(w);
00682 Node z=level_list[l];
00683 if ( z!=INVALID ) left.set(z,w);
00684 right.set(w,z);
00685 level_list[l]=w;
00686 level.set(w, l);
00687 }
00688 }
00689 }
00690 }
00691
00692
00693 switch (flow_prop) {
00694 case NO_FLOW:
00695 for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0);
00696 case ZERO_FLOW:
00697 for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
00698
00699
00700 level.set(_target,0);
00701 bfs_queue.push(_target);
00702
00703 while ( !bfs_queue.empty() ) {
00704
00705 Node v=bfs_queue.front();
00706 bfs_queue.pop();
00707 int l=level[v]+1;
00708
00709 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
00710 Node w=_g->source(e);
00711 if ( level[w] == _node_num && w != _source ) {
00712 bfs_queue.push(w);
00713 Node z=level_list[l];
00714 if ( z!=INVALID ) left.set(z,w);
00715 right.set(w,z);
00716 level_list[l]=w;
00717 level.set(w, l);
00718 }
00719 }
00720 }
00721
00722
00723 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
00724 Num c=(*_capacity)[e];
00725 if ( c <= 0 ) continue;
00726 Node w=_g->target(e);
00727 if ( level[w] < _node_num ) {
00728 if ( excess[w] <= 0 && w!=_target ) {
00729 next.set(w,first[level[w]]);
00730 first[level[w]]=w;
00731 }
00732 _flow->set(e, c);
00733 excess.set(w, excess[w]+c);
00734 }
00735 }
00736 break;
00737
00738 case GEN_FLOW:
00739 for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
00740 {
00741 Num exc=0;
00742 for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e];
00743 for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e];
00744 excess.set(_target,exc);
00745 }
00746
00747
00748 for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e) {
00749 Num rem=(*_capacity)[e]-(*_flow)[e];
00750 if ( rem <= 0 ) continue;
00751 Node w=_g->target(e);
00752 if ( level[w] < _node_num ) {
00753 if ( excess[w] <= 0 && w!=_target ) {
00754 next.set(w,first[level[w]]);
00755 first[level[w]]=w;
00756 }
00757 _flow->set(e, (*_capacity)[e]);
00758 excess.set(w, excess[w]+rem);
00759 }
00760 }
00761
00762 for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) {
00763 if ( (*_flow)[e] <= 0 ) continue;
00764 Node w=_g->source(e);
00765 if ( level[w] < _node_num ) {
00766 if ( excess[w] <= 0 && w!=_target ) {
00767 next.set(w,first[level[w]]);
00768 first[level[w]]=w;
00769 }
00770 excess.set(w, excess[w]+(*_flow)[e]);
00771 _flow->set(e, 0);
00772 }
00773 }
00774 break;
00775
00776 case PRE_FLOW:
00777
00778 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
00779 Num rem=(*_capacity)[e]-(*_flow)[e];
00780 if ( rem <= 0 ) continue;
00781 Node w=_g->target(e);
00782 if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]);
00783 }
00784
00785 for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
00786 if ( (*_flow)[e] <= 0 ) continue;
00787 Node w=_g->source(e);
00788 if ( level[w] < _node_num ) _flow->set(e, 0);
00789 }
00790
00791
00792 for(NodeIt w(*_g); w!=INVALID; ++w) {
00793 Num exc=0;
00794 for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e];
00795 for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e];
00796 excess.set(w,exc);
00797
00798
00799 int lev=level[w];
00800 if ( exc > 0 && lev < _node_num && Node(w) != _target ) {
00801 next.set(w,first[lev]);
00802 first[lev]=w;
00803 }
00804 }
00805 break;
00806 }
00807 }
00808
00809
00810 void relabel(Node w, int newlevel, VecNode& first, NNMap& next,
00811 VecNode& level_list, NNMap& left,
00812 NNMap& right, int& b, int& k, bool what_heur )
00813 {
00814
00815 int lev=level[w];
00816
00817 Node right_n=right[w];
00818 Node left_n=left[w];
00819
00820
00821 if ( right_n!=INVALID ) {
00822 if ( left_n!=INVALID ) {
00823 right.set(left_n, right_n);
00824 left.set(right_n, left_n);
00825 } else {
00826 level_list[lev]=right_n;
00827 left.set(right_n, INVALID);
00828 }
00829 } else {
00830 if ( left_n!=INVALID ) {
00831 right.set(left_n, INVALID);
00832 } else {
00833 level_list[lev]=INVALID;
00834 }
00835 }
00836
00837
00838 if ( level_list[lev]==INVALID ) {
00839
00840
00841 for (int i=lev; i!=k ; ) {
00842 Node v=level_list[++i];
00843 while ( v!=INVALID ) {
00844 level.set(v,_node_num);
00845 v=right[v];
00846 }
00847 level_list[i]=INVALID;
00848 if ( !what_heur ) first[i]=INVALID;
00849 }
00850
00851 level.set(w,_node_num);
00852 b=lev-1;
00853 k=b;
00854
00855
00856 } else {
00857
00858 if ( newlevel == _node_num ) level.set(w,_node_num);
00859 else {
00860 level.set(w,++newlevel);
00861 next.set(w,first[newlevel]);
00862 first[newlevel]=w;
00863 if ( what_heur ) b=newlevel;
00864 if ( k < newlevel ) ++k;
00865 Node z=level_list[newlevel];
00866 if ( z!=INVALID ) left.set(z,w);
00867 right.set(w,z);
00868 left.set(w,INVALID);
00869 level_list[newlevel]=w;
00870 }
00871 }
00872 }
00873
00874 };
00875
00881 template<class GR, class CM, class FM>
00882 Preflow<GR,typename CM::Value,CM,FM> preflow(const GR &g,
00883 typename GR::Node source,
00884 typename GR::Node target,
00885 const CM &cap,
00886 FM &flow
00887 )
00888 {
00889 return Preflow<GR,typename CM::Value,CM,FM>(g,source,target,cap,flow);
00890 }
00891
00892 }
00893
00894 #endif //LEMON_PREFLOW_H