37 ///\ingroup flowalgs |
37 ///\ingroup flowalgs |
38 ///\brief %Preflow algorithms class. |
38 ///\brief %Preflow algorithms class. |
39 /// |
39 /// |
40 ///This class provides an implementation of the \e preflow \e |
40 ///This class provides an implementation of the \e preflow \e |
41 ///algorithm producing a flow of maximum value in a directed |
41 ///algorithm producing a flow of maximum value in a directed |
42 ///graph. The preflow algorithms are the fastest known max flow algorithms |
42 ///graph. The preflow algorithms are the fastest known max flow algorithms. |
43 ///up to now. The \e source node, the \e target node, the \e |
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 |
44 ///capacity of the edges and the \e starting \e flow value of the |
45 ///edges should be passed to the algorithm through the |
45 ///edges should be passed to the algorithm through the |
46 ///constructor. It is possible to change these quantities using the |
46 ///constructor. It is possible to change these quantities using the |
47 ///functions \ref source, \ref target, \ref capacityMap and \ref |
47 ///functions \ref source, \ref target, \ref capacityMap and \ref |
48 ///flowMap. |
48 ///flowMap. |
57 /// |
57 /// |
58 ///\param Graph The directed graph type the algorithm runs on. |
58 ///\param Graph The directed graph type the algorithm runs on. |
59 ///\param Num The number type of the capacities and the flow values. |
59 ///\param Num The number type of the capacities and the flow values. |
60 ///\param CapacityMap The capacity map type. |
60 ///\param CapacityMap The capacity map type. |
61 ///\param FlowMap The flow map type. |
61 ///\param FlowMap The flow map type. |
|
62 ///\param Tolerance The tolerance type. |
62 /// |
63 /// |
63 ///\author Jacint Szabo |
64 ///\author Jacint Szabo |
64 ///\todo Second template parameter is superfluous |
65 ///\todo Second template parameter is superfluous |
65 ///\todo Using tolerance |
|
66 template <typename Graph, typename Num, |
66 template <typename Graph, typename Num, |
67 typename CapacityMap=typename Graph::template EdgeMap<Num>, |
67 typename CapacityMap=typename Graph::template EdgeMap<Num>, |
68 typename FlowMap=typename Graph::template EdgeMap<Num>, |
68 typename FlowMap=typename Graph::template EdgeMap<Num>, |
69 typename Tolerance=Tolerance<Num> > |
69 typename Tolerance=Tolerance<Num> > |
70 class Preflow { |
70 class Preflow { |
245 bool end=false; |
245 bool end=false; |
246 //Needed for 'bound decrease', true means no active |
246 //Needed for 'bound decrease', true means no active |
247 //nodes are above bound b. |
247 //nodes are above bound b. |
248 |
248 |
249 int k=_node_num-2; //bound on the highest level under n containing a node |
249 int k=_node_num-2; //bound on the highest level under n containing a node |
250 int b=k; //bound on the highest level under n of an active node |
250 int b=k; //bound on the highest level under n containing an active node |
251 |
251 |
252 VecNode first(_node_num, INVALID); |
252 VecNode first(_node_num, INVALID); |
253 NNMap next(*_g, INVALID); |
253 NNMap next(*_g, INVALID); |
254 |
254 |
255 NNMap left(*_g, INVALID); |
255 NNMap left(*_g, INVALID); |
272 else { |
272 else { |
273 end=false; |
273 end=false; |
274 Node w=first[b]; |
274 Node w=first[b]; |
275 first[b]=next[w]; |
275 first[b]=next[w]; |
276 int newlevel=push(w, next, first); |
276 int newlevel=push(w, next, first); |
277 if ( excess[w] > 0 ) relabel(w, newlevel, first, next, level_list, |
277 if ( _surely.positive(excess[w]) ) relabel(w, newlevel, first, next, level_list, |
278 left, right, b, k, what_heur); |
278 left, right, b, k, what_heur); |
279 |
279 |
280 ++numrelabel; |
280 ++numrelabel; |
281 if ( numrelabel >= heur ) { |
281 if ( numrelabel >= heur ) { |
282 numrelabel=0; |
282 numrelabel=0; |
334 Node v=bfs_queue.front(); |
334 Node v=bfs_queue.front(); |
335 bfs_queue.pop(); |
335 bfs_queue.pop(); |
336 int l=level[v]+1; |
336 int l=level[v]+1; |
337 |
337 |
338 for(InEdgeIt e(*_g,v); e!=INVALID; ++e) { |
338 for(InEdgeIt e(*_g,v); e!=INVALID; ++e) { |
339 if ( (*_capacity)[e] <= (*_flow)[e] ) continue; |
339 if ( !_surely.less((*_flow)[e], (*_capacity)[e]) ) continue; |
340 Node u=_g->source(e); |
340 Node u=_g->source(e); |
341 if ( level[u] >= _node_num ) { |
341 if ( level[u] >= _node_num ) { |
342 bfs_queue.push(u); |
342 bfs_queue.push(u); |
343 level.set(u, l); |
343 level.set(u, l); |
344 if ( excess[u] > 0 ) { |
344 if ( _surely.positive(excess[u]) ) { |
345 next.set(u,first[l]); |
345 next.set(u,first[l]); |
346 first[l]=u; |
346 first[l]=u; |
347 } |
347 } |
348 } |
348 } |
349 } |
349 } |
350 |
350 |
351 for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) { |
351 for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) { |
352 if ( 0 >= (*_flow)[e] ) continue; |
352 if ( !_surely.positive((*_flow)[e]) ) continue; |
353 Node u=_g->target(e); |
353 Node u=_g->target(e); |
354 if ( level[u] >= _node_num ) { |
354 if ( level[u] >= _node_num ) { |
355 bfs_queue.push(u); |
355 bfs_queue.push(u); |
356 level.set(u, l); |
356 level.set(u, l); |
357 if ( excess[u] > 0 ) { |
357 if ( _surely.positive(excess[u]) ) { |
358 next.set(u,first[l]); |
358 next.set(u,first[l]); |
359 first[l]=u; |
359 first[l]=u; |
360 } |
360 } |
361 } |
361 } |
362 } |
362 } |
371 Node w=first[b]; |
371 Node w=first[b]; |
372 first[b]=next[w]; |
372 first[b]=next[w]; |
373 int newlevel=push(w,next, first); |
373 int newlevel=push(w,next, first); |
374 |
374 |
375 //relabel |
375 //relabel |
376 if ( excess[w] > 0 ) { |
376 if ( _surely.positive(excess[w]) ) { |
377 level.set(w,++newlevel); |
377 level.set(w,++newlevel); |
378 next.set(w,first[newlevel]); |
378 next.set(w,first[newlevel]); |
379 first[newlevel]=w; |
379 first[newlevel]=w; |
380 b=newlevel; |
380 b=newlevel; |
381 } |
381 } |
441 Node w=queue.front(); |
441 Node w=queue.front(); |
442 queue.pop(); |
442 queue.pop(); |
443 |
443 |
444 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
444 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
445 Node v=_g->target(e); |
445 Node v=_g->target(e); |
446 if (!M[v] && (*_flow)[e] < (*_capacity)[e] ) { |
446 if (!M[v] && _surely.less((*_flow)[e] , (*_capacity)[e]) ) { |
447 queue.push(v); |
447 queue.push(v); |
448 M.set(v, true); |
448 M.set(v, true); |
449 } |
449 } |
450 } |
450 } |
451 |
451 |
452 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
452 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
453 Node v=_g->source(e); |
453 Node v=_g->source(e); |
454 if (!M[v] && (*_flow)[e] > 0 ) { |
454 if (!M[v] && _surely.positive((*_flow)[e]) ) { |
455 queue.push(v); |
455 queue.push(v); |
456 M.set(v, true); |
456 M.set(v, true); |
457 } |
457 } |
458 } |
458 } |
459 } |
459 } |
479 Node w=queue.front(); |
479 Node w=queue.front(); |
480 queue.pop(); |
480 queue.pop(); |
481 |
481 |
482 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
482 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
483 Node v=_g->source(e); |
483 Node v=_g->source(e); |
484 if (M[v] && (*_flow)[e] < (*_capacity)[e] ) { |
484 if (M[v] && _surely.less((*_flow)[e], (*_capacity)[e]) ) { |
485 queue.push(v); |
485 queue.push(v); |
486 M.set(v, false); |
486 M.set(v, false); |
487 } |
487 } |
488 } |
488 } |
489 |
489 |
490 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
490 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
491 Node v=_g->target(e); |
491 Node v=_g->target(e); |
492 if (M[v] && (*_flow)[e] > 0 ) { |
492 if (M[v] && _surely.positive((*_flow)[e]) ) { |
493 queue.push(v); |
493 queue.push(v); |
494 M.set(v, false); |
494 M.set(v, false); |
495 } |
495 } |
496 } |
496 } |
497 } |
497 } |
574 int lev=level[w]; |
574 int lev=level[w]; |
575 Num exc=excess[w]; |
575 Num exc=excess[w]; |
576 int newlevel=_node_num; //bound on the next level of w |
576 int newlevel=_node_num; //bound on the next level of w |
577 |
577 |
578 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
578 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
579 if ( (*_flow)[e] >= (*_capacity)[e] ) continue; |
579 if ( !_surely.less((*_flow)[e], (*_capacity)[e]) ) continue; |
580 Node v=_g->target(e); |
580 Node v=_g->target(e); |
581 |
581 |
582 if( lev > level[v] ) { //Push is allowed now |
582 if( lev > level[v] ) { //Push is allowed now |
583 |
583 |
584 if ( excess[v]<=0 && v!=_target && v!=_source ) { |
584 if ( !_surely.positive(excess[v]) && v!=_target && v!=_source ) { |
585 next.set(v,first[level[v]]); |
585 next.set(v,first[level[v]]); |
586 first[level[v]]=v; |
586 first[level[v]]=v; |
587 } |
587 } |
588 |
588 |
589 Num cap=(*_capacity)[e]; |
589 Num cap=(*_capacity)[e]; |
590 Num flo=(*_flow)[e]; |
590 Num flo=(*_flow)[e]; |
591 Num remcap=cap-flo; |
591 Num remcap=cap-flo; |
592 |
592 |
593 if ( remcap >= exc ) { //A nonsaturating push. |
593 if ( ! _surely.less(remcap, exc) ) { //A nonsaturating push. |
594 |
594 |
595 _flow->set(e, flo+exc); |
595 _flow->set(e, flo+exc); |
596 excess.set(v, excess[v]+exc); |
596 excess.set(v, excess[v]+exc); |
597 exc=0; |
597 exc=0; |
598 break; |
598 break; |
603 exc-=remcap; |
603 exc-=remcap; |
604 } |
604 } |
605 } else if ( newlevel > level[v] ) newlevel = level[v]; |
605 } else if ( newlevel > level[v] ) newlevel = level[v]; |
606 } //for out edges wv |
606 } //for out edges wv |
607 |
607 |
608 if ( exc > 0 ) { |
608 if ( _surely.positive(exc) ) { |
609 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
609 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
610 |
610 |
611 if( (*_flow)[e] <= 0 ) continue; |
611 if ( !_surely.positive((*_flow)[e]) ) continue; |
612 Node v=_g->source(e); |
612 Node v=_g->source(e); |
613 |
613 |
614 if( lev > level[v] ) { //Push is allowed now |
614 if( lev > level[v] ) { //Push is allowed now |
615 |
615 |
616 if ( excess[v]<=0 && v!=_target && v!=_source ) { |
616 if ( !_surely.positive(excess[v]) && v!=_target && v!=_source ) { |
617 next.set(v,first[level[v]]); |
617 next.set(v,first[level[v]]); |
618 first[level[v]]=v; |
618 first[level[v]]=v; |
619 } |
619 } |
620 |
620 |
621 Num flo=(*_flow)[e]; |
621 Num flo=(*_flow)[e]; |
622 |
622 |
623 if ( flo >= exc ) { //A nonsaturating push. |
623 if ( !_surely.less(flo, exc) ) { //A nonsaturating push. |
624 |
624 |
625 _flow->set(e, flo-exc); |
625 _flow->set(e, flo-exc); |
626 excess.set(v, excess[v]+exc); |
626 excess.set(v, excess[v]+exc); |
627 exc=0; |
627 exc=0; |
628 break; |
628 break; |
661 Node v=bfs_queue.front(); |
661 Node v=bfs_queue.front(); |
662 bfs_queue.pop(); |
662 bfs_queue.pop(); |
663 int l=level[v]+1; |
663 int l=level[v]+1; |
664 |
664 |
665 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) { |
665 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) { |
666 if ( (*_capacity)[e] <= (*_flow)[e] ) continue; |
666 if ( !_surely.less((*_flow)[e],(*_capacity)[e]) ) continue; |
667 Node w=_g->source(e); |
667 Node w=_g->source(e); |
668 if ( level[w] == _node_num && w != _source ) { |
668 if ( level[w] == _node_num && w != _source ) { |
669 bfs_queue.push(w); |
669 bfs_queue.push(w); |
670 Node z=level_list[l]; |
670 Node z=level_list[l]; |
671 if ( z!=INVALID ) left.set(z,w); |
671 if ( z!=INVALID ) left.set(z,w); |
674 level.set(w, l); |
674 level.set(w, l); |
675 } |
675 } |
676 } |
676 } |
677 |
677 |
678 for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) { |
678 for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) { |
679 if ( 0 >= (*_flow)[e] ) continue; |
679 if ( !_surely.positive((*_flow)[e]) ) continue; |
680 Node w=_g->target(e); |
680 Node w=_g->target(e); |
681 if ( level[w] == _node_num && w != _source ) { |
681 if ( level[w] == _node_num && w != _source ) { |
682 bfs_queue.push(w); |
682 bfs_queue.push(w); |
683 Node z=level_list[l]; |
683 Node z=level_list[l]; |
684 if ( z!=INVALID ) left.set(z,w); |
684 if ( z!=INVALID ) left.set(z,w); |
721 } |
721 } |
722 |
722 |
723 //the starting flow |
723 //the starting flow |
724 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { |
724 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { |
725 Num c=(*_capacity)[e]; |
725 Num c=(*_capacity)[e]; |
726 if ( c <= 0 ) continue; |
726 if ( !_surely.positive(c) ) continue; |
727 Node w=_g->target(e); |
727 Node w=_g->target(e); |
728 if ( level[w] < _node_num ) { |
728 if ( level[w] < _node_num ) { |
729 if ( excess[w] <= 0 && w!=_target ) { //putting into the stack |
729 if ( !_surely.positive(excess[w]) && w!=_target ) { //putting into the stack |
730 next.set(w,first[level[w]]); |
730 next.set(w,first[level[w]]); |
731 first[level[w]]=w; |
731 first[level[w]]=w; |
732 } |
732 } |
733 _flow->set(e, c); |
733 _flow->set(e, c); |
734 excess.set(w, excess[w]+c); |
734 excess.set(w, excess[w]+c); |
746 } |
746 } |
747 |
747 |
748 //the starting flow |
748 //the starting flow |
749 for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e) { |
749 for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e) { |
750 Num rem=(*_capacity)[e]-(*_flow)[e]; |
750 Num rem=(*_capacity)[e]-(*_flow)[e]; |
751 if ( rem <= 0 ) continue; |
751 if ( !_surely.positive(rem) ) continue; |
752 Node w=_g->target(e); |
752 Node w=_g->target(e); |
753 if ( level[w] < _node_num ) { |
753 if ( level[w] < _node_num ) { |
754 if ( excess[w] <= 0 && w!=_target ) { //putting into the stack |
754 if ( !_surely.positive(excess[w]) && w!=_target ) { //putting into the stack |
755 next.set(w,first[level[w]]); |
755 next.set(w,first[level[w]]); |
756 first[level[w]]=w; |
756 first[level[w]]=w; |
757 } |
757 } |
758 _flow->set(e, (*_capacity)[e]); |
758 _flow->set(e, (*_capacity)[e]); |
759 excess.set(w, excess[w]+rem); |
759 excess.set(w, excess[w]+rem); |
760 } |
760 } |
761 } |
761 } |
762 |
762 |
763 for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) { |
763 for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) { |
764 if ( (*_flow)[e] <= 0 ) continue; |
764 if ( !_surely.positive((*_flow)[e]) ) continue; |
765 Node w=_g->source(e); |
765 Node w=_g->source(e); |
766 if ( level[w] < _node_num ) { |
766 if ( level[w] < _node_num ) { |
767 if ( excess[w] <= 0 && w!=_target ) { |
767 if ( !_surely.positive(excess[w]) && w!=_target ) { |
768 next.set(w,first[level[w]]); |
768 next.set(w,first[level[w]]); |
769 first[level[w]]=w; |
769 first[level[w]]=w; |
770 } |
770 } |
771 excess.set(w, excess[w]+(*_flow)[e]); |
771 excess.set(w, excess[w]+(*_flow)[e]); |
772 _flow->set(e, 0); |
772 _flow->set(e, 0); |
776 |
776 |
777 case PRE_FLOW: |
777 case PRE_FLOW: |
778 //the starting flow |
778 //the starting flow |
779 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { |
779 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { |
780 Num rem=(*_capacity)[e]-(*_flow)[e]; |
780 Num rem=(*_capacity)[e]-(*_flow)[e]; |
781 if ( rem <= 0 ) continue; |
781 if ( !_surely.positive(rem) ) continue; |
782 Node w=_g->target(e); |
782 Node w=_g->target(e); |
783 if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]); |
783 if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]); |
784 } |
784 } |
785 |
785 |
786 for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { |
786 for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { |
787 if ( (*_flow)[e] <= 0 ) continue; |
787 if ( !_surely.positive((*_flow)[e]) ) continue; |
788 Node w=_g->source(e); |
788 Node w=_g->source(e); |
789 if ( level[w] < _node_num ) _flow->set(e, 0); |
789 if ( level[w] < _node_num ) _flow->set(e, 0); |
790 } |
790 } |
791 |
791 |
792 //computing the excess |
792 //computing the excess |
796 for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e]; |
796 for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e]; |
797 excess.set(w,exc); |
797 excess.set(w,exc); |
798 |
798 |
799 //putting the active nodes into the stack |
799 //putting the active nodes into the stack |
800 int lev=level[w]; |
800 int lev=level[w]; |
801 if ( exc > 0 && lev < _node_num && Node(w) != _target ) { |
801 if ( _surely.positive(exc) && lev < _node_num && Node(w) != _target ) { |
802 next.set(w,first[lev]); |
802 next.set(w,first[lev]); |
803 first[lev]=w; |
803 first[lev]=w; |
804 } |
804 } |
805 } |
805 } |
806 break; |
806 break; |