95 MaxFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity, |
96 MaxFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity, |
96 FlowMap& _flow) : |
97 FlowMap& _flow) : |
97 g(&_G), s(_s), t(_t), capacity(&_capacity), |
98 g(&_G), s(_s), t(_t), capacity(&_capacity), |
98 flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {} |
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 void run() { |
103 void run() { |
101 preflow( ZERO_FLOW ); |
104 preflow( ZERO_FLOW ); |
102 } |
105 } |
103 |
106 |
|
107 /// A preflow algorithm is run. The initial edge-set have to be a flow, |
|
108 /// or from a preflow, according to \c fe. |
104 void preflow( flowEnum fe ) { |
109 void preflow( flowEnum fe ) { |
105 preflowPhase0(fe); |
110 preflowPhase0(fe); |
106 preflowPhase1(); |
111 preflowPhase1(); |
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 void preflowPhase0( flowEnum fe ); |
116 void preflowPhase0( flowEnum fe ); |
110 |
117 |
|
118 /// Second phase of preflow. |
111 void preflowPhase1(); |
119 void preflowPhase1(); |
112 |
120 |
|
121 /// Starting from a flow, this method searches for an augmenting path |
|
122 /// according to the Edmonds-Karp algorithm |
|
123 /// and augments the flow on if any. |
113 bool augmentOnShortestPath(); |
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 template<typename MutableGraph> bool augmentOnBlockingFlow(); |
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 bool augmentOnBlockingFlow2(); |
134 bool augmentOnBlockingFlow2(); |
118 |
135 |
119 /// Returns the actual flow value. |
136 /// Returns the actual flow value. |
120 /// More precisely, it returns the negative excess of s, thus |
137 /// More precisely, it returns the negative excess of s, thus |
121 /// this works also for preflows. |
138 /// this works also for preflows. |
312 } // if w still has excess after the out edge for cycle |
328 } // if w still has excess after the out edge for cycle |
313 |
329 |
314 excess.set(w, exc); |
330 excess.set(w, exc); |
315 |
331 |
316 return newlevel; |
332 return newlevel; |
317 } |
333 } |
318 |
334 |
319 |
335 |
320 void preflowPreproc ( flowEnum fe, VecStack& active, |
336 void preflowPreproc ( flowEnum fe, VecStack& active, |
321 VecNode& level_list, NNMap& left, NNMap& right ) { |
337 VecNode& level_list, NNMap& left, NNMap& right ) { |
322 |
338 |
323 std::queue<Node> bfs_queue; |
339 std::queue<Node> bfs_queue; |
324 |
340 |
325 switch ( fe ) { |
341 switch ( fe ) { |
326 case ZERO_FLOW: |
342 case ZERO_FLOW: |
327 { |
343 { |
328 //Reverse_bfs from t, to find the starting level. |
344 //Reverse_bfs from t, to find the starting level. |
329 level.set(t,0); |
345 level.set(t,0); |
330 bfs_queue.push(t); |
346 bfs_queue.push(t); |
331 |
347 |
332 while (!bfs_queue.empty()) { |
348 while (!bfs_queue.empty()) { |
333 |
349 |
334 Node v=bfs_queue.front(); |
350 Node v=bfs_queue.front(); |
335 bfs_queue.pop(); |
351 bfs_queue.pop(); |
336 int l=level[v]+1; |
352 int l=level[v]+1; |
337 |
353 |
338 InEdgeIt e; |
354 InEdgeIt e; |
339 for(g->first(e,v); g->valid(e); g->next(e)) { |
355 for(g->first(e,v); g->valid(e); g->next(e)) { |
340 Node w=g->tail(e); |
356 Node w=g->tail(e); |
341 if ( level[w] == n && w != s ) { |
357 if ( level[w] == n && w != s ) { |
342 bfs_queue.push(w); |
358 bfs_queue.push(w); |
343 Node first=level_list[l]; |
359 Node first=level_list[l]; |
344 if ( g->valid(first) ) left.set(first,w); |
360 if ( g->valid(first) ) left.set(first,w); |
345 right.set(w,first); |
361 right.set(w,first); |
346 level_list[l]=w; |
362 level_list[l]=w; |
347 level.set(w, l); |
363 level.set(w, l); |
348 } |
364 } |
349 } |
365 } |
350 } |
366 } |
351 |
367 |
352 //the starting flow |
368 //the starting flow |
353 OutEdgeIt e; |
369 OutEdgeIt e; |
354 for(g->first(e,s); g->valid(e); g->next(e)) |
370 for(g->first(e,s); g->valid(e); g->next(e)) |
355 { |
371 { |
356 Num c=(*capacity)[e]; |
372 Num c=(*capacity)[e]; |
357 if ( c <= 0 ) continue; |
373 if ( c <= 0 ) continue; |
358 Node w=g->head(e); |
374 Node w=g->head(e); |
359 if ( level[w] < n ) { |
375 if ( level[w] < n ) { |
360 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w); |
376 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w); |
361 flow->set(e, c); |
377 flow->set(e, c); |
362 excess.set(w, excess[w]+c); |
378 excess.set(w, excess[w]+c); |
363 } |
379 } |
364 } |
380 } |
365 break; |
381 break; |
366 } |
382 } |
367 |
383 |
368 case GEN_FLOW: |
384 case GEN_FLOW: |
369 case PREFLOW: |
385 case PREFLOW: |
370 { |
386 { |
371 //Reverse_bfs from t in the residual graph, |
387 //Reverse_bfs from t in the residual graph, |
372 //to find the starting level. |
388 //to find the starting level. |
373 level.set(t,0); |
389 level.set(t,0); |
374 bfs_queue.push(t); |
390 bfs_queue.push(t); |
375 |
391 |
376 while (!bfs_queue.empty()) { |
392 while (!bfs_queue.empty()) { |
377 |
393 |
378 Node v=bfs_queue.front(); |
394 Node v=bfs_queue.front(); |
379 bfs_queue.pop(); |
395 bfs_queue.pop(); |
380 int l=level[v]+1; |
396 int l=level[v]+1; |
381 |
397 |
382 InEdgeIt e; |
398 InEdgeIt e; |
383 for(g->first(e,v); g->valid(e); g->next(e)) { |
399 for(g->first(e,v); g->valid(e); g->next(e)) { |
384 if ( (*capacity)[e] <= (*flow)[e] ) continue; |
400 if ( (*capacity)[e] <= (*flow)[e] ) continue; |
385 Node w=g->tail(e); |
401 Node w=g->tail(e); |
386 if ( level[w] == n && w != s ) { |
402 if ( level[w] == n && w != s ) { |
387 bfs_queue.push(w); |
403 bfs_queue.push(w); |
388 Node first=level_list[l]; |
404 Node first=level_list[l]; |
389 if ( g->valid(first) ) left.set(first,w); |
405 if ( g->valid(first) ) left.set(first,w); |
390 right.set(w,first); |
406 right.set(w,first); |
391 level_list[l]=w; |
407 level_list[l]=w; |
392 level.set(w, l); |
408 level.set(w, l); |
393 } |
409 } |
394 } |
410 } |
395 |
411 |
396 OutEdgeIt f; |
412 OutEdgeIt f; |
397 for(g->first(f,v); g->valid(f); g->next(f)) { |
413 for(g->first(f,v); g->valid(f); g->next(f)) { |
398 if ( 0 >= (*flow)[f] ) continue; |
414 if ( 0 >= (*flow)[f] ) continue; |
399 Node w=g->head(f); |
415 Node w=g->head(f); |
400 if ( level[w] == n && w != s ) { |
416 if ( level[w] == n && w != s ) { |
401 bfs_queue.push(w); |
417 bfs_queue.push(w); |
402 Node first=level_list[l]; |
418 Node first=level_list[l]; |
403 if ( g->valid(first) ) left.set(first,w); |
419 if ( g->valid(first) ) left.set(first,w); |
404 right.set(w,first); |
420 right.set(w,first); |
405 level_list[l]=w; |
421 level_list[l]=w; |
406 level.set(w, l); |
422 level.set(w, l); |
407 } |
423 } |
408 } |
424 } |
409 } |
425 } |
410 |
426 |
411 |
427 |
412 //the starting flow |
428 //the starting flow |
413 OutEdgeIt e; |
429 OutEdgeIt e; |
414 for(g->first(e,s); g->valid(e); g->next(e)) |
430 for(g->first(e,s); g->valid(e); g->next(e)) |
415 { |
431 { |
416 Num rem=(*capacity)[e]-(*flow)[e]; |
432 Num rem=(*capacity)[e]-(*flow)[e]; |
417 if ( rem <= 0 ) continue; |
433 if ( rem <= 0 ) continue; |
418 Node w=g->head(e); |
434 Node w=g->head(e); |
419 if ( level[w] < n ) { |
435 if ( level[w] < n ) { |
420 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w); |
436 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w); |
421 flow->set(e, (*capacity)[e]); |
437 flow->set(e, (*capacity)[e]); |
422 excess.set(w, excess[w]+rem); |
438 excess.set(w, excess[w]+rem); |
423 } |
439 } |
424 } |
440 } |
425 |
441 |
426 InEdgeIt f; |
442 InEdgeIt f; |
427 for(g->first(f,s); g->valid(f); g->next(f)) |
443 for(g->first(f,s); g->valid(f); g->next(f)) |
428 { |
444 { |
429 if ( (*flow)[f] <= 0 ) continue; |
445 if ( (*flow)[f] <= 0 ) continue; |
430 Node w=g->tail(f); |
446 Node w=g->tail(f); |
431 if ( level[w] < n ) { |
447 if ( level[w] < n ) { |
432 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w); |
448 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w); |
433 excess.set(w, excess[w]+(*flow)[f]); |
449 excess.set(w, excess[w]+(*flow)[f]); |
434 flow->set(f, 0); |
450 flow->set(f, 0); |
435 } |
451 } |
436 } |
452 } |
437 break; |
453 break; |
438 } //case PREFLOW |
454 } //case PREFLOW |
439 } |
455 } |
440 } //preflowPreproc |
456 } //preflowPreproc |
441 |
457 |
442 |
458 |
443 |
459 |
444 void relabel(Node w, int newlevel, VecStack& active, |
460 void relabel(Node w, int newlevel, VecStack& active, |
445 VecNode& level_list, NNMap& left, |
461 VecNode& level_list, NNMap& left, |
536 |
552 |
537 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
553 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
538 void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase0( flowEnum fe ) |
554 void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase0( flowEnum fe ) |
539 { |
555 { |
540 |
556 |
541 int heur0=(int)(H0*n); //time while running 'bound decrease' |
557 int heur0=(int)(H0*n); //time while running 'bound decrease' |
542 int heur1=(int)(H1*n); //time while running 'highest label' |
558 int heur1=(int)(H1*n); //time while running 'highest label' |
543 int heur=heur1; //starting time interval (#of relabels) |
559 int heur=heur1; //starting time interval (#of relabels) |
544 int numrelabel=0; |
560 int numrelabel=0; |
545 |
561 |
546 bool what_heur=1; |
562 bool what_heur=1; |
547 //It is 0 in case 'bound decrease' and 1 in case 'highest label' |
563 //It is 0 in case 'bound decrease' and 1 in case 'highest label' |
548 |
564 |
549 bool end=false; |
565 bool end=false; |
550 //Needed for 'bound decrease', true means no active nodes are above bound b. |
566 //Needed for 'bound decrease', true means no active nodes are above bound b. |
551 |
567 |
552 int k=n-2; //bound on the highest level under n containing a node |
568 int k=n-2; //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 |
569 int b=k; //bound on the highest level under n of an active node |
554 |
570 |
555 VecStack active(n); |
571 VecStack active(n); |
556 |
572 |
557 NNMap left(*g, INVALID); |
573 NNMap left(*g, INVALID); |
558 NNMap right(*g, INVALID); |
574 NNMap right(*g, INVALID); |
559 VecNode level_list(n,INVALID); |
575 VecNode level_list(n,INVALID); |
560 //List of the nodes in level i<n, set to n. |
576 //List of the nodes in level i<n, set to n. |
561 |
577 |
562 NodeIt v; |
578 NodeIt v; |
563 for(g->first(v); g->valid(v); g->next(v)) level.set(v,n); |
579 for(g->first(v); g->valid(v); g->next(v)) level.set(v,n); |
564 //setting each node to level n |
580 //setting each node to level n |
565 |
581 |
566 switch ( fe ) { |
582 switch ( fe ) { |
567 case PREFLOW: |
583 case PREFLOW: |
568 { |
584 { |
569 //counting the excess |
585 //counting the excess |
570 NodeIt v; |
586 NodeIt v; |
571 for(g->first(v); g->valid(v); g->next(v)) { |
587 for(g->first(v); g->valid(v); g->next(v)) { |
572 Num exc=0; |
588 Num exc=0; |
573 |
589 |
574 InEdgeIt e; |
590 InEdgeIt e; |
575 for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e]; |
591 for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e]; |
576 OutEdgeIt f; |
592 OutEdgeIt f; |
577 for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f]; |
593 for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f]; |
578 |
594 |
579 excess.set(v,exc); |
595 excess.set(v,exc); |
580 |
596 |
581 //putting the active nodes into the stack |
597 //putting the active nodes into the stack |
582 int lev=level[v]; |
598 int lev=level[v]; |
583 if ( exc > 0 && lev < n && v != t ) active[lev].push(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; |
655 } |
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 } |
|
641 } |
656 } |
642 } |
657 } |
|
658 } |
643 |
659 |
644 |
660 |
645 |
661 |
646 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
662 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
647 void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase1() |
663 void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase1() |
648 { |
664 { |
649 |
665 |
650 int k=n-2; //bound on the highest level under n containing a node |
666 int k=n-2; //bound on the highest level under n containing a node |
651 int b=k; //bound on the highest level under n of an active node |
667 int b=k; //bound on the highest level under n of an active node |
652 |
668 |
653 VecStack active(n); |
669 VecStack active(n); |
654 level.set(s,0); |
670 level.set(s,0); |
655 std::queue<Node> bfs_queue; |
671 std::queue<Node> bfs_queue; |
656 bfs_queue.push(s); |
672 bfs_queue.push(s); |
657 |
673 |
658 while (!bfs_queue.empty()) { |
674 while (!bfs_queue.empty()) { |
659 |
675 |
660 Node v=bfs_queue.front(); |
676 Node v=bfs_queue.front(); |
661 bfs_queue.pop(); |
677 bfs_queue.pop(); |
662 int l=level[v]+1; |
678 int l=level[v]+1; |
663 |
679 |
664 InEdgeIt e; |
680 InEdgeIt e; |
665 for(g->first(e,v); g->valid(e); g->next(e)) { |
681 for(g->first(e,v); g->valid(e); g->next(e)) { |
666 if ( (*capacity)[e] <= (*flow)[e] ) continue; |
682 if ( (*capacity)[e] <= (*flow)[e] ) continue; |
667 Node u=g->tail(e); |
683 Node u=g->tail(e); |
668 if ( level[u] >= n ) { |
684 if ( level[u] >= n ) { |
669 bfs_queue.push(u); |
685 bfs_queue.push(u); |
670 level.set(u, l); |
686 level.set(u, l); |
671 if ( excess[u] > 0 ) active[l].push(u); |
687 if ( excess[u] > 0 ) active[l].push(u); |
672 } |
688 } |
673 } |
689 } |
674 |
690 |
675 OutEdgeIt f; |
691 OutEdgeIt f; |
676 for(g->first(f,v); g->valid(f); g->next(f)) { |
692 for(g->first(f,v); g->valid(f); g->next(f)) { |
677 if ( 0 >= (*flow)[f] ) continue; |
693 if ( 0 >= (*flow)[f] ) continue; |
678 Node u=g->head(f); |
694 Node u=g->head(f); |
679 if ( level[u] >= n ) { |
695 if ( level[u] >= n ) { |
680 bfs_queue.push(u); |
696 bfs_queue.push(u); |
681 level.set(u, l); |
697 level.set(u, l); |
682 if ( excess[u] > 0 ) active[l].push(u); |
698 if ( excess[u] > 0 ) active[l].push(u); |
683 } |
699 } |
684 } |
700 } |
685 } |
701 } |
686 b=n-2; |
702 b=n-2; |
687 |
703 |
688 while ( true ) { |
704 while ( true ) { |
689 |
705 |
690 if ( b == 0 ) break; |
706 if ( b == 0 ) break; |
691 |
707 |
692 if ( active[b].empty() ) --b; |
708 if ( active[b].empty() ) --b; |
693 else { |
709 else { |
694 Node w=active[b].top(); |
710 Node w=active[b].top(); |
695 active[b].pop(); |
711 active[b].pop(); |
696 int newlevel=push(w,active); |
712 int newlevel=push(w,active); |
697 |
713 |
698 //relabel |
714 //relabel |
699 if ( excess[w] > 0 ) { |
715 if ( excess[w] > 0 ) { |
700 level.set(w,++newlevel); |
716 level.set(w,++newlevel); |
701 active[newlevel].push(w); |
717 active[newlevel].push(w); |
702 b=newlevel; |
718 b=newlevel; |
703 } |
719 } |
704 } // if stack[b] is nonempty |
720 } // if stack[b] is nonempty |
705 } // while(true) |
721 } // while(true) |
706 } |
722 } |
707 |
723 |
708 |
724 |
709 |
725 |
710 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
726 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
711 bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath() |
727 bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath() |
712 { |
728 { |
713 ResGW res_graph(*g, *capacity, *flow); |
729 ResGW res_graph(*g, *capacity, *flow); |
714 bool _augment=false; |
730 bool _augment=false; |
715 |
731 |
716 //ReachedMap level(res_graph); |
732 //ReachedMap level(res_graph); |
717 FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); |
733 FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); |
718 BfsIterator<ResGW, ReachedMap> bfs(res_graph, level); |
734 BfsIterator<ResGW, ReachedMap> bfs(res_graph, level); |
719 bfs.pushAndSetReached(s); |
735 bfs.pushAndSetReached(s); |
720 |
736 |
721 typename ResGW::template NodeMap<ResGWEdge> pred(res_graph); |
737 typename ResGW::template NodeMap<ResGWEdge> pred(res_graph); |
722 pred.set(s, INVALID); |
738 pred.set(s, INVALID); |
723 |
739 |
724 typename ResGW::template NodeMap<Num> free(res_graph); |
740 typename ResGW::template NodeMap<Num> free(res_graph); |
725 |
741 |
726 //searching for augmenting path |
742 //searching for augmenting path |
727 while ( !bfs.finished() ) { |
743 while ( !bfs.finished() ) { |
728 ResGWOutEdgeIt e=bfs; |
744 ResGWOutEdgeIt e=bfs; |
729 if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) { |
745 if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) { |
730 Node v=res_graph.tail(e); |
746 Node v=res_graph.tail(e); |
731 Node w=res_graph.head(e); |
747 Node w=res_graph.head(e); |
732 pred.set(w, e); |
748 pred.set(w, e); |
733 if (res_graph.valid(pred[v])) { |
749 if (res_graph.valid(pred[v])) { |
734 free.set(w, std::min(free[v], res_graph.resCap(e))); |
750 free.set(w, std::min(free[v], res_graph.resCap(e))); |
735 } else { |
751 } else { |
736 free.set(w, res_graph.resCap(e)); |
752 free.set(w, res_graph.resCap(e)); |
737 } |
753 } |
738 if (res_graph.head(e)==t) { _augment=true; break; } |
754 if (res_graph.head(e)==t) { _augment=true; break; } |
739 } |
755 } |
740 |
756 |
741 ++bfs; |
757 ++bfs; |
742 } //end of searching augmenting path |
758 } //end of searching augmenting path |
743 |
759 |
744 if (_augment) { |
760 if (_augment) { |
745 Node n=t; |
761 Node n=t; |
746 Num augment_value=free[t]; |
762 Num augment_value=free[t]; |
747 while (res_graph.valid(pred[n])) { |
763 while (res_graph.valid(pred[n])) { |
748 ResGWEdge e=pred[n]; |
764 ResGWEdge e=pred[n]; |
749 res_graph.augment(e, augment_value); |
765 res_graph.augment(e, augment_value); |
750 n=res_graph.tail(e); |
766 n=res_graph.tail(e); |
751 } |
767 } |
752 } |
768 } |
753 |
769 |
754 return _augment; |
770 return _augment; |
755 } |
771 } |
756 |
772 |
757 |
773 |
758 |
774 |
759 |
775 |
760 |
776 |
764 |
780 |
765 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
781 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
766 template<typename MutableGraph> |
782 template<typename MutableGraph> |
767 bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow() |
783 bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow() |
768 { |
784 { |
769 typedef MutableGraph MG; |
785 typedef MutableGraph MG; |
770 bool _augment=false; |
786 bool _augment=false; |
771 |
787 |
772 ResGW res_graph(*g, *capacity, *flow); |
788 ResGW res_graph(*g, *capacity, *flow); |
773 |
789 |
774 //bfs for distances on the residual graph |
790 //bfs for distances on the residual graph |
775 //ReachedMap level(res_graph); |
791 //ReachedMap level(res_graph); |
776 FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); |
792 FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); |
777 BfsIterator<ResGW, ReachedMap> bfs(res_graph, level); |
793 BfsIterator<ResGW, ReachedMap> bfs(res_graph, level); |
778 bfs.pushAndSetReached(s); |
794 bfs.pushAndSetReached(s); |
779 typename ResGW::template NodeMap<int> |
795 typename ResGW::template NodeMap<int> |
780 dist(res_graph); //filled up with 0's |
796 dist(res_graph); //filled up with 0's |
781 |
797 |
782 //F will contain the physical copy of the residual graph |
798 //F will contain the physical copy of the residual graph |
783 //with the set of edges which are on shortest paths |
799 //with the set of edges which are on shortest paths |
784 MG F; |
800 MG F; |
785 typename ResGW::template NodeMap<typename MG::Node> |
801 typename ResGW::template NodeMap<typename MG::Node> |
786 res_graph_to_F(res_graph); |
802 res_graph_to_F(res_graph); |
787 { |
803 { |
788 typename ResGW::NodeIt n; |
804 typename ResGW::NodeIt n; |
789 for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) { |
805 for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) { |
790 res_graph_to_F.set(n, F.addNode()); |
806 res_graph_to_F.set(n, F.addNode()); |
791 } |
807 } |
792 } |
808 } |
793 |
809 |
794 typename MG::Node sF=res_graph_to_F[s]; |
810 typename MG::Node sF=res_graph_to_F[s]; |
795 typename MG::Node tF=res_graph_to_F[t]; |
811 typename MG::Node tF=res_graph_to_F[t]; |
796 typename MG::template EdgeMap<ResGWEdge> original_edge(F); |
812 typename MG::template EdgeMap<ResGWEdge> original_edge(F); |
797 typename MG::template EdgeMap<Num> residual_capacity(F); |
813 typename MG::template EdgeMap<Num> residual_capacity(F); |
798 |
814 |
799 while ( !bfs.finished() ) { |
815 while ( !bfs.finished() ) { |
800 ResGWOutEdgeIt e=bfs; |
816 ResGWOutEdgeIt e=bfs; |
801 if (res_graph.valid(e)) { |
817 if (res_graph.valid(e)) { |
802 if (bfs.isBNodeNewlyReached()) { |
818 if (bfs.isBNodeNewlyReached()) { |
803 dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1); |
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 typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]); |
827 typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]); |
805 original_edge.update(); |
828 original_edge.update(); |
806 original_edge.set(f, e); |
829 original_edge.set(f, e); |
807 residual_capacity.update(); |
830 residual_capacity.update(); |
808 residual_capacity.set(f, res_graph.resCap(e)); |
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 } else { |
869 } else { |
810 if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) { |
870 F.erase(/*typename MG::OutEdgeIt*/(dfs)); |
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 } |
|
817 } |
871 } |
818 } |
872 } |
819 ++bfs; |
873 } |
820 } //computing distances from s in the residual graph |
874 |
821 |
875 if (__augment) { |
822 bool __augment=true; |
876 typename MG::Node n=tF; |
823 |
877 Num augment_value=free[tF]; |
824 while (__augment) { |
878 while (F.valid(pred[n])) { |
825 __augment=false; |
879 typename MG::Edge e=pred[n]; |
826 //computing blocking flow with dfs |
880 res_graph.augment(original_edge[e], augment_value); |
827 DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F); |
881 n=F.tail(e); |
828 typename MG::template NodeMap<typename MG::Edge> pred(F); |
882 if (residual_capacity[e]==augment_value) |
829 pred.set(sF, INVALID); |
883 F.erase(e); |
830 //invalid iterators for sources |
884 else |
831 |
885 residual_capacity.set(e, residual_capacity[e]-augment_value); |
832 typename MG::template NodeMap<Num> free(F); |
886 } |
833 |
887 } |
834 dfs.pushAndSetReached(sF); |
888 |
835 while (!dfs.finished()) { |
889 } |
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 } |
|
874 |
890 |
875 return _augment; |
891 return _augment; |
876 } |
892 } |
877 |
893 |
878 |
894 |
879 |
895 |
880 |
896 |
881 |
897 |
882 |
898 |
883 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
899 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
884 bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2() |
900 bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2() |
885 { |
901 { |
886 bool _augment=false; |
902 bool _augment=false; |
887 |
903 |
888 ResGW res_graph(*g, *capacity, *flow); |
904 ResGW res_graph(*g, *capacity, *flow); |
889 |
905 |
890 //ReachedMap level(res_graph); |
906 //ReachedMap level(res_graph); |
891 FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); |
907 FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); |
892 BfsIterator<ResGW, ReachedMap> bfs(res_graph, level); |
908 BfsIterator<ResGW, ReachedMap> bfs(res_graph, level); |
893 |
909 |
894 bfs.pushAndSetReached(s); |
910 bfs.pushAndSetReached(s); |
895 DistanceMap<ResGW> dist(res_graph); |
911 DistanceMap<ResGW> dist(res_graph); |
896 while ( !bfs.finished() ) { |
912 while ( !bfs.finished() ) { |
897 ResGWOutEdgeIt e=bfs; |
913 ResGWOutEdgeIt e=bfs; |
898 if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) { |
914 if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) { |
899 dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1); |
915 dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1); |
900 } |
916 } |
901 ++bfs; |
917 ++bfs; |
902 } //computing distances from s in the residual graph |
918 } //computing distances from s in the residual graph |
903 |
919 |
904 //Subgraph containing the edges on some shortest paths |
920 //Subgraph containing the edges on some shortest paths |
905 ConstMap<typename ResGW::Node, bool> true_map(true); |
921 ConstMap<typename ResGW::Node, bool> true_map(true); |
906 typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>, |
922 typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>, |
907 DistanceMap<ResGW> > FilterResGW; |
923 DistanceMap<ResGW> > FilterResGW; |
908 FilterResGW filter_res_graph(res_graph, true_map, dist); |
924 FilterResGW filter_res_graph(res_graph, true_map, dist); |
909 |
925 |
910 //Subgraph, which is able to delete edges which are already |
926 //Subgraph, which is able to delete edges which are already |
911 //met by the dfs |
927 //met by the dfs |
912 typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt> |
928 typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt> |
913 first_out_edges(filter_res_graph); |
929 first_out_edges(filter_res_graph); |
914 typename FilterResGW::NodeIt v; |
930 typename FilterResGW::NodeIt v; |
915 for(filter_res_graph.first(v); filter_res_graph.valid(v); |
931 for(filter_res_graph.first(v); filter_res_graph.valid(v); |
916 filter_res_graph.next(v)) |
932 filter_res_graph.next(v)) |
917 { |
933 { |
918 typename FilterResGW::OutEdgeIt e; |
934 typename FilterResGW::OutEdgeIt e; |
919 filter_res_graph.first(e, v); |
935 filter_res_graph.first(e, v); |
920 first_out_edges.set(v, e); |
936 first_out_edges.set(v, e); |
921 } |
937 } |
922 typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW:: |
938 typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW:: |
923 template NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW; |
939 template NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW; |
924 ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges); |
940 ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges); |
925 |
941 |
926 bool __augment=true; |
942 bool __augment=true; |
927 |
943 |
928 while (__augment) { |
944 while (__augment) { |
929 |
945 |
930 __augment=false; |
946 __augment=false; |
931 //computing blocking flow with dfs |
947 //computing blocking flow with dfs |
932 DfsIterator< ErasingResGW, |
948 DfsIterator< ErasingResGW, |
933 typename ErasingResGW::template NodeMap<bool> > |
949 typename ErasingResGW::template NodeMap<bool> > |
934 dfs(erasing_res_graph); |
950 dfs(erasing_res_graph); |
935 typename ErasingResGW:: |
951 typename ErasingResGW:: |
936 template NodeMap<typename ErasingResGW::OutEdgeIt> |
952 template NodeMap<typename ErasingResGW::OutEdgeIt> |
937 pred(erasing_res_graph); |
953 pred(erasing_res_graph); |
938 pred.set(s, INVALID); |
954 pred.set(s, INVALID); |
939 //invalid iterators for sources |
955 //invalid iterators for sources |
940 |
956 |
941 typename ErasingResGW::template NodeMap<Num> |
957 typename ErasingResGW::template NodeMap<Num> |
942 free1(erasing_res_graph); |
958 free1(erasing_res_graph); |
943 |
959 |
944 dfs.pushAndSetReached( |
960 dfs.pushAndSetReached( |
945 typename ErasingResGW::Node( |
961 typename ErasingResGW::Node( |
946 typename FilterResGW::Node( |
962 typename FilterResGW::Node( |
947 typename ResGW::Node(s) |
963 typename ResGW::Node(s) |
948 ) |
964 ) |
949 ) |
965 ) |
950 ); |
966 ); |
951 while (!dfs.finished()) { |
967 while (!dfs.finished()) { |
952 ++dfs; |
968 ++dfs; |
953 if (erasing_res_graph.valid( |
969 if (erasing_res_graph.valid( |
954 typename ErasingResGW::OutEdgeIt(dfs))) |
970 typename ErasingResGW::OutEdgeIt(dfs))) |
955 { |
971 { |
956 if (dfs.isBNodeNewlyReached()) { |
972 if (dfs.isBNodeNewlyReached()) { |
957 |
973 |
958 typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs); |
974 typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs); |
959 typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs); |
975 typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs); |
960 |
976 |
961 pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs)); |
977 pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs)); |
962 if (erasing_res_graph.valid(pred[v])) { |
978 if (erasing_res_graph.valid(pred[v])) { |
963 free1.set(w, std::min(free1[v], res_graph.resCap( |
979 free1.set(w, std::min(free1[v], res_graph.resCap( |
964 typename ErasingResGW::OutEdgeIt(dfs)))); |
980 typename ErasingResGW::OutEdgeIt(dfs)))); |
965 } else { |
981 } else { |
966 free1.set(w, res_graph.resCap( |
982 free1.set(w, res_graph.resCap( |
967 typename ErasingResGW::OutEdgeIt(dfs))); |
983 typename ErasingResGW::OutEdgeIt(dfs))); |
968 } |
984 } |
969 |
985 |
970 if (w==t) { |
986 if (w==t) { |
971 __augment=true; |
987 __augment=true; |
972 _augment=true; |
988 _augment=true; |