COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/max_flow.h @ 586:04fdffd38e89

Last change on this file since 586:04fdffd38e89 was 586:04fdffd38e89, checked in by Alpar Juttner, 16 years ago

doc

File size: 27.2 KB
Line 
1// -*- C++ -*-
2
3/*
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 
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
34*/
35
36#ifndef HUGO_MAX_FLOW_H
37#define HUGO_MAX_FLOW_H
38
39#define H0 20
40#define H1 1
41
42#include <vector>
43#include <queue>
44#include <stack>
45
46#include <hugo/graph_wrapper.h>
47#include <bfs_iterator.h>
48#include <hugo/invalid.h>
49#include <hugo/maps.h>
50#include <for_each_macros.h>
51
52/// \file
53/// \brief Dimacs file format reader.
54
55namespace hugo {
56
57
58//  ///\author Marton Makai, Jacint Szabo
59  /// A class for computing max flows and related quantities.
60  template <typename Graph, typename Num,
61            typename CapMap=typename Graph::template EdgeMap<Num>,
62            typename FlowMap=typename Graph::template EdgeMap<Num> >
63  class MaxFlow {
64   
65    typedef typename Graph::Node Node;
66    typedef typename Graph::NodeIt NodeIt;
67    typedef typename Graph::OutEdgeIt OutEdgeIt;
68    typedef typename Graph::InEdgeIt InEdgeIt;
69
70    typedef typename std::vector<std::stack<Node> > VecStack;
71    typedef typename Graph::template NodeMap<Node> NNMap;
72    typedef typename std::vector<Node> VecNode;
73
74    const Graph* g;
75    Node s;
76    Node t;
77    const CapMap* capacity; 
78    FlowMap* flow;
79    int n;      //the number of nodes of G
80    typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
81    typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
82    typedef typename ResGW::Edge ResGWEdge;
83    //typedef typename ResGW::template NodeMap<bool> ReachedMap;
84    typedef typename Graph::template NodeMap<int> ReachedMap;
85    ReachedMap level;
86    //level works as a bool map in augmenting path algorithms
87    //and is used by bfs for storing reached information.
88    //In preflow, it shows levels of nodes.
89    //typename Graph::template NodeMap<int> level;   
90    typename Graph::template NodeMap<Num> excess;
91//   protected:
92//     MaxFlow() { }
93//     void set(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,
94//           FlowMap& _flow)
95//       {
96//      g=&_G;
97//      s=_s;
98//      t=_t;
99//      capacity=&_capacity;
100//      flow=&_flow;
101//      n=_G.nodeNum;
102//      level.set (_G); //kellene vmi ilyesmi fv
103//      excess(_G,0); //itt is
104//       }
105
106  public:
107 
108    ///\todo Document this.
109    ///\todo Maybe, it should be PRE_FLOW instead.
110    ///- \c ZERO_FLOW means something,
111    ///- \c GEN_FLOW means something else,
112    ///- \c PREFLOW is something different.
113    enum flowEnum{
114      ZERO_FLOW=0,
115      GEN_FLOW=1,
116      PREFLOW=2
117    };
118
119    MaxFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,
120            FlowMap& _flow) :
121      g(&_G), s(_s), t(_t), capacity(&_capacity),
122      flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {}
123
124    /// A max flow algorithm is run.
125    ///\pre the flow have to be 0 at the beginning.
126    void run() {
127      preflow(ZERO_FLOW);
128    }
129   
130    /// A preflow algorithm is run.
131    ///\pre The initial edge-map have to be a
132    /// zero flow if \c fe is \c ZERO_FLOW,
133    /// a flow if \c fe is \c GEN_FLOW,
134    /// and a pre-flow it is \c PREFLOW.
135    void preflow(flowEnum fe) {
136      preflowPhase0(fe);
137      preflowPhase1();
138    }
139
140    /// Run the first phase of preflow, starting from a 0 flow, from a flow,
141    /// or from a preflow, according to \c fe.
142    void preflowPhase0( flowEnum fe );
143
144    /// Second phase of preflow.
145    void preflowPhase1();
146
147    /// Starting from a flow, this method searches for an augmenting path
148    /// according to the Edmonds-Karp algorithm
149    /// and augments the flow on if any.
150    /// The return value shows if the augmentation was succesful.
151    bool augmentOnShortestPath();
152
153    /// Starting from a flow, this method searches for an augmenting blockin
154    /// flow according to Dinits' algorithm and augments the flow on if any.
155    /// The blocking flow is computed in a physically constructed
156    /// residual graph of type \c Mutablegraph.
157    /// The return value show sif the augmentation was succesful.
158    template<typename MutableGraph> bool augmentOnBlockingFlow();
159
160    /// The same as \c augmentOnBlockingFlow<MutableGraph> but the
161    /// residual graph is not constructed physically.
162    /// The return value shows if the augmentation was succesful.
163    bool augmentOnBlockingFlow2();
164
165    /// Returns the actual flow value.
166    /// More precisely, it returns the negative excess of s, thus
167    /// this works also for preflows.
168    Num flowValue() {
169      Num a=0;
170      FOR_EACH_INC_LOC(OutEdgeIt, e, *g, s) a+=(*flow)[e];
171      FOR_EACH_INC_LOC(InEdgeIt, e, *g, s) a-=(*flow)[e];
172      return a;
173    }
174
175    /// Should be used between preflowPhase0 and preflowPhase1.
176    ///\todo We have to make some status variable which shows the actual state
177    /// of the class. This enables us to determine which methods are valid
178    /// for MinCut computation
179    template<typename _CutMap>
180    void actMinCut(_CutMap& M) {
181      NodeIt v;
182      for(g->first(v); g->valid(v); g->next(v)) {
183        if ( level[v] < n ) {
184          M.set(v,false);
185        } else {
186          M.set(v,true);
187        }
188      }
189    }
190
191    /// The unique inclusionwise minimum cut is computed by
192    /// processing a bfs from s in the residual graph.
193    ///\pre flow have to be a max flow otherwise it will the whole node-set.
194    template<typename _CutMap>
195    void minMinCut(_CutMap& M) {
196   
197      std::queue<Node> queue;
198     
199      M.set(s,true);     
200      queue.push(s);
201
202      while (!queue.empty()) {
203        Node w=queue.front();
204        queue.pop();
205
206        OutEdgeIt e;
207        for(g->first(e,w) ; g->valid(e); g->next(e)) {
208          Node v=g->head(e);
209          if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
210            queue.push(v);
211            M.set(v, true);
212          }
213        }
214
215        InEdgeIt f;
216        for(g->first(f,w) ; g->valid(f); g->next(f)) {
217          Node v=g->tail(f);
218          if (!M[v] && (*flow)[f] > 0 ) {
219            queue.push(v);
220            M.set(v, true);
221          }
222        }
223      }
224    }
225
226
227    /// The unique inclusionwise maximum cut is computed by
228    /// processing a reverse bfs from t in the residual graph.
229    ///\pre flow have to be a max flow otherwise it will be empty.
230    template<typename _CutMap>
231    void maxMinCut(_CutMap& M) {
232
233      NodeIt v;
234      for(g->first(v) ; g->valid(v); g->next(v)) {
235        M.set(v, true);
236      }
237
238      std::queue<Node> queue;
239     
240      M.set(t,false);       
241      queue.push(t);
242
243      while (!queue.empty()) {
244        Node w=queue.front();
245        queue.pop();
246
247
248        InEdgeIt e;
249        for(g->first(e,w) ; g->valid(e); g->next(e)) {
250          Node v=g->tail(e);
251          if (M[v] && (*flow)[e] < (*capacity)[e] ) {
252            queue.push(v);
253            M.set(v, false);
254          }
255        }
256       
257        OutEdgeIt f;
258        for(g->first(f,w) ; g->valid(f); g->next(f)) {
259          Node v=g->head(f);
260          if (M[v] && (*flow)[f] > 0 ) {
261            queue.push(v);
262            M.set(v, false);
263          }
264        }
265      }
266    }
267
268
269    /// A minimum cut is computed.
270    template<typename CutMap>
271    void minCut(CutMap& M) { minMinCut(M); }
272
273    ///
274    void resetSource(Node _s) { s=_s; }
275    ///
276    void resetTarget(Node _t) { t=_t; }
277   
278    /// capacity-map is changed.
279    void resetCap(const CapMap& _cap) { capacity=&_cap; }
280   
281    /// flow-map is changed.
282    void resetFlow(FlowMap& _flow) { flow=&_flow; }
283
284
285  private:
286
287    int push(Node w, VecStack& active) {
288     
289      int lev=level[w];
290      Num exc=excess[w];
291      int newlevel=n;       //bound on the next level of w
292         
293      OutEdgeIt e;
294      for(g->first(e,w); g->valid(e); g->next(e)) {
295           
296        if ( (*flow)[e] >= (*capacity)[e] ) continue;
297        Node v=g->head(e);           
298           
299        if( lev > level[v] ) { //Push is allowed now
300         
301          if ( excess[v]<=0 && v!=t && v!=s ) {
302            int lev_v=level[v];
303            active[lev_v].push(v);
304          }
305         
306          Num cap=(*capacity)[e];
307          Num flo=(*flow)[e];
308          Num remcap=cap-flo;
309         
310          if ( remcap >= exc ) { //A nonsaturating push.
311           
312            flow->set(e, flo+exc);
313            excess.set(v, excess[v]+exc);
314            exc=0;
315            break;
316           
317          } else { //A saturating push.
318            flow->set(e, cap);
319            excess.set(v, excess[v]+remcap);
320            exc-=remcap;
321          }
322        } else if ( newlevel > level[v] ) newlevel = level[v];
323      } //for out edges wv
324     
325      if ( exc > 0 ) { 
326        InEdgeIt e;
327        for(g->first(e,w); g->valid(e); g->next(e)) {
328         
329          if( (*flow)[e] <= 0 ) continue;
330          Node v=g->tail(e);
331         
332          if( lev > level[v] ) { //Push is allowed now
333           
334            if ( excess[v]<=0 && v!=t && v!=s ) {
335              int lev_v=level[v];
336              active[lev_v].push(v);
337            }
338           
339            Num flo=(*flow)[e];
340           
341            if ( flo >= exc ) { //A nonsaturating push.
342             
343              flow->set(e, flo-exc);
344              excess.set(v, excess[v]+exc);
345              exc=0;
346              break;
347            } else {  //A saturating push.
348             
349              excess.set(v, excess[v]+flo);
350              exc-=flo;
351              flow->set(e,0);
352            } 
353          } else if ( newlevel > level[v] ) newlevel = level[v];
354        } //for in edges vw
355       
356      } // if w still has excess after the out edge for cycle
357     
358      excess.set(w, exc);
359     
360      return newlevel;
361    }
362
363
364    void preflowPreproc ( flowEnum fe, VecStack& active,
365                          VecNode& level_list, NNMap& left, NNMap& right ) {
366
367                            std::queue<Node> bfs_queue;
368     
369                            switch ( fe ) {
370                            case ZERO_FLOW:
371                              {
372                                //Reverse_bfs from t, to find the starting level.
373                                level.set(t,0);
374                                bfs_queue.push(t);
375       
376                                while (!bfs_queue.empty()) {
377           
378                                  Node v=bfs_queue.front();     
379                                  bfs_queue.pop();
380                                  int l=level[v]+1;
381           
382                                  InEdgeIt e;
383                                  for(g->first(e,v); g->valid(e); g->next(e)) {
384                                    Node w=g->tail(e);
385                                    if ( level[w] == n && w != s ) {
386                                      bfs_queue.push(w);
387                                      Node first=level_list[l];
388                                      if ( g->valid(first) ) left.set(first,w);
389                                      right.set(w,first);
390                                      level_list[l]=w;
391                                      level.set(w, l);
392                                    }
393                                  }
394                                }
395         
396                                //the starting flow
397                                OutEdgeIt e;
398                                for(g->first(e,s); g->valid(e); g->next(e))
399                                  {
400                                    Num c=(*capacity)[e];
401                                    if ( c <= 0 ) continue;
402                                    Node w=g->head(e);
403                                    if ( level[w] < n ) {         
404                                      if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
405                                      flow->set(e, c);
406                                      excess.set(w, excess[w]+c);
407                                    }
408                                  }
409                                break;
410                              }
411       
412                            case GEN_FLOW:
413                            case PREFLOW:
414                              {
415                                //Reverse_bfs from t in the residual graph,
416                                //to find the starting level.
417                                level.set(t,0);
418                                bfs_queue.push(t);
419         
420                                while (!bfs_queue.empty()) {
421           
422                                  Node v=bfs_queue.front();     
423                                  bfs_queue.pop();
424                                  int l=level[v]+1;
425           
426                                  InEdgeIt e;
427                                  for(g->first(e,v); g->valid(e); g->next(e)) {
428                                    if ( (*capacity)[e] <= (*flow)[e] ) continue;
429                                    Node w=g->tail(e);
430                                    if ( level[w] == n && w != s ) {
431                                      bfs_queue.push(w);
432                                      Node first=level_list[l];
433                                      if ( g->valid(first) ) left.set(first,w);
434                                      right.set(w,first);
435                                      level_list[l]=w;
436                                      level.set(w, l);
437                                    }
438                                  }
439           
440                                  OutEdgeIt f;
441                                  for(g->first(f,v); g->valid(f); g->next(f)) {
442                                    if ( 0 >= (*flow)[f] ) continue;
443                                    Node w=g->head(f);
444                                    if ( level[w] == n && w != s ) {
445                                      bfs_queue.push(w);
446                                      Node first=level_list[l];
447                                      if ( g->valid(first) ) left.set(first,w);
448                                      right.set(w,first);
449                                      level_list[l]=w;
450                                      level.set(w, l);
451                                    }
452                                  }
453                                }
454         
455         
456                                //the starting flow
457                                OutEdgeIt e;
458                                for(g->first(e,s); g->valid(e); g->next(e))
459                                  {
460                                    Num rem=(*capacity)[e]-(*flow)[e];
461                                    if ( rem <= 0 ) continue;
462                                    Node w=g->head(e);
463                                    if ( level[w] < n ) {         
464                                      if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
465                                      flow->set(e, (*capacity)[e]);
466                                      excess.set(w, excess[w]+rem);
467                                    }
468                                  }
469         
470                                InEdgeIt f;
471                                for(g->first(f,s); g->valid(f); g->next(f))
472                                  {
473                                    if ( (*flow)[f] <= 0 ) continue;
474                                    Node w=g->tail(f);
475                                    if ( level[w] < n ) {         
476                                      if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
477                                      excess.set(w, excess[w]+(*flow)[f]);
478                                      flow->set(f, 0);
479                                    }
480                                  } 
481                                break;
482                              } //case PREFLOW
483                            }
484                          } //preflowPreproc
485
486
487
488    void relabel(Node w, int newlevel, VecStack& active, 
489                 VecNode& level_list, NNMap& left,
490                 NNMap& right, int& b, int& k, bool what_heur )
491    {
492
493      Num lev=level[w];
494     
495      Node right_n=right[w];
496      Node left_n=left[w];
497     
498      //unlacing starts
499      if ( g->valid(right_n) ) {
500        if ( g->valid(left_n) ) {
501          right.set(left_n, right_n);
502          left.set(right_n, left_n);
503        } else {
504          level_list[lev]=right_n;   
505          left.set(right_n, INVALID);
506        }
507      } else {
508        if ( g->valid(left_n) ) {
509          right.set(left_n, INVALID);
510        } else {
511          level_list[lev]=INVALID;   
512        }
513      }
514      //unlacing ends
515               
516      if ( !g->valid(level_list[lev]) ) {
517             
518        //gapping starts
519        for (int i=lev; i!=k ; ) {
520          Node v=level_list[++i];
521          while ( g->valid(v) ) {
522            level.set(v,n);
523            v=right[v];
524          }
525          level_list[i]=INVALID;
526          if ( !what_heur ) {
527            while ( !active[i].empty() ) {
528              active[i].pop();    //FIXME: ezt szebben kene
529            }
530          }         
531        }
532       
533        level.set(w,n);
534        b=lev-1;
535        k=b;
536        //gapping ends
537       
538      } else {
539       
540        if ( newlevel == n ) level.set(w,n);
541        else {
542          level.set(w,++newlevel);
543          active[newlevel].push(w);
544          if ( what_heur ) b=newlevel;
545          if ( k < newlevel ) ++k;      //now k=newlevel
546          Node first=level_list[newlevel];
547          if ( g->valid(first) ) left.set(first,w);
548          right.set(w,first);
549          left.set(w,INVALID);
550          level_list[newlevel]=w;
551        }
552      }
553     
554    } //relabel
555
556
557    template<typename MapGraphWrapper>
558    class DistanceMap {
559    protected:
560      const MapGraphWrapper* g;
561      typename MapGraphWrapper::template NodeMap<int> dist;
562    public:
563      DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { }
564      void set(const typename MapGraphWrapper::Node& n, int a) {
565        dist.set(n, a);
566      }
567      int operator[](const typename MapGraphWrapper::Node& n)
568      { return dist[n]; }
569      //       int get(const typename MapGraphWrapper::Node& n) const {
570      //        return dist[n]; }
571      //       bool get(const typename MapGraphWrapper::Edge& e) const {
572      //        return (dist.get(g->tail(e))<dist.get(g->head(e))); }
573      bool operator[](const typename MapGraphWrapper::Edge& e) const {
574        return (dist[g->tail(e)]<dist[g->head(e)]);
575      }
576    };
577   
578  };
579
580
581  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
582  void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase0( flowEnum fe )
583  {
584     
585    int heur0=(int)(H0*n);  //time while running 'bound decrease'
586    int heur1=(int)(H1*n);  //time while running 'highest label'
587    int heur=heur1;         //starting time interval (#of relabels)
588    int numrelabel=0;
589     
590    bool what_heur=1;       
591    //It is 0 in case 'bound decrease' and 1 in case 'highest label'
592
593    bool end=false;     
594    //Needed for 'bound decrease', true means no active nodes are above bound b.
595
596    int k=n-2;  //bound on the highest level under n containing a node
597    int b=k;    //bound on the highest level under n of an active node
598     
599    VecStack active(n);
600     
601    NNMap left(*g, INVALID);
602    NNMap right(*g, INVALID);
603    VecNode level_list(n,INVALID);
604    //List of the nodes in level i<n, set to n.
605
606    NodeIt v;
607    for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
608    //setting each node to level n
609     
610    switch ( fe ) {
611    case PREFLOW:
612      {
613        //counting the excess
614        NodeIt v;
615        for(g->first(v); g->valid(v); g->next(v)) {
616          Num exc=0;
617         
618          InEdgeIt e;
619          for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
620          OutEdgeIt f;
621          for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
622           
623          excess.set(v,exc);     
624           
625          //putting the active nodes into the stack
626          int lev=level[v];
627          if ( exc > 0 && lev < n && v != t ) active[lev].push(v);
628        }
629        break;
630      }
631    case GEN_FLOW:
632      {
633        //Counting the excess of t
634        Num exc=0;
635         
636        InEdgeIt e;
637        for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
638        OutEdgeIt f;
639        for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
640         
641        excess.set(t,exc);     
642         
643        break;
644      }
645    default:
646      break;
647    }
648     
649    preflowPreproc( fe, active, level_list, left, right );
650    //End of preprocessing
651     
652     
653    //Push/relabel on the highest level active nodes.
654    while ( true ) {
655      if ( b == 0 ) {
656        if ( !what_heur && !end && k > 0 ) {
657          b=k;
658          end=true;
659        } else break;
660      }
661       
662      if ( active[b].empty() ) --b;
663      else {
664        end=false; 
665        Node w=active[b].top();
666        active[b].pop();
667        int newlevel=push(w,active);
668        if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list,
669                                     left, right, b, k, what_heur);
670         
671        ++numrelabel;
672        if ( numrelabel >= heur ) {
673          numrelabel=0;
674          if ( what_heur ) {
675            what_heur=0;
676            heur=heur0;
677            end=false;
678          } else {
679            what_heur=1;
680            heur=heur1;
681            b=k;
682          }
683        }
684      }
685    }
686  }
687
688
689
690  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
691  void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase1()
692  {
693     
694    int k=n-2;  //bound on the highest level under n containing a node
695    int b=k;    //bound on the highest level under n of an active node
696     
697    VecStack active(n);
698    level.set(s,0);
699    std::queue<Node> bfs_queue;
700    bfs_queue.push(s);
701           
702    while (!bfs_queue.empty()) {
703       
704      Node v=bfs_queue.front();
705      bfs_queue.pop();
706      int l=level[v]+1;
707             
708      InEdgeIt e;
709      for(g->first(e,v); g->valid(e); g->next(e)) {
710        if ( (*capacity)[e] <= (*flow)[e] ) continue;
711        Node u=g->tail(e);
712        if ( level[u] >= n ) {
713          bfs_queue.push(u);
714          level.set(u, l);
715          if ( excess[u] > 0 ) active[l].push(u);
716        }
717      }
718       
719      OutEdgeIt f;
720      for(g->first(f,v); g->valid(f); g->next(f)) {
721        if ( 0 >= (*flow)[f] ) continue;
722        Node u=g->head(f);
723        if ( level[u] >= n ) {
724          bfs_queue.push(u);
725          level.set(u, l);
726          if ( excess[u] > 0 ) active[l].push(u);
727        }
728      }
729    }
730    b=n-2;
731
732    while ( true ) {
733       
734      if ( b == 0 ) break;
735
736      if ( active[b].empty() ) --b;
737      else {
738        Node w=active[b].top();
739        active[b].pop();
740        int newlevel=push(w,active);     
741
742        //relabel
743        if ( excess[w] > 0 ) {
744          level.set(w,++newlevel);
745          active[newlevel].push(w);
746          b=newlevel;
747        }
748      }  // if stack[b] is nonempty
749    } // while(true)
750  }
751
752
753
754  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
755  bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath()
756  {
757    ResGW res_graph(*g, *capacity, *flow);
758    bool _augment=false;
759     
760    //ReachedMap level(res_graph);
761    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
762    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
763    bfs.pushAndSetReached(s);
764       
765    typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
766    pred.set(s, INVALID);
767     
768    typename ResGW::template NodeMap<Num> free(res_graph);
769       
770    //searching for augmenting path
771    while ( !bfs.finished() ) {
772      ResGWOutEdgeIt e=bfs;
773      if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
774        Node v=res_graph.tail(e);
775        Node w=res_graph.head(e);
776        pred.set(w, e);
777        if (res_graph.valid(pred[v])) {
778          free.set(w, std::min(free[v], res_graph.resCap(e)));
779        } else {
780          free.set(w, res_graph.resCap(e));
781        }
782        if (res_graph.head(e)==t) { _augment=true; break; }
783      }
784       
785      ++bfs;
786    } //end of searching augmenting path
787
788    if (_augment) {
789      Node n=t;
790      Num augment_value=free[t];
791      while (res_graph.valid(pred[n])) {
792        ResGWEdge e=pred[n];
793        res_graph.augment(e, augment_value);
794        n=res_graph.tail(e);
795      }
796    }
797
798    return _augment;
799  }
800
801
802
803
804
805
806
807
808
809  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
810  template<typename MutableGraph>
811  bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow()
812  {     
813    typedef MutableGraph MG;
814    bool _augment=false;
815
816    ResGW res_graph(*g, *capacity, *flow);
817
818    //bfs for distances on the residual graph
819    //ReachedMap level(res_graph);
820    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
821    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
822    bfs.pushAndSetReached(s);
823    typename ResGW::template NodeMap<int>
824      dist(res_graph); //filled up with 0's
825
826    //F will contain the physical copy of the residual graph
827    //with the set of edges which are on shortest paths
828    MG F;
829    typename ResGW::template NodeMap<typename MG::Node>
830      res_graph_to_F(res_graph);
831    {
832      typename ResGW::NodeIt n;
833      for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
834        res_graph_to_F.set(n, F.addNode());
835      }
836    }
837
838    typename MG::Node sF=res_graph_to_F[s];
839    typename MG::Node tF=res_graph_to_F[t];
840    typename MG::template EdgeMap<ResGWEdge> original_edge(F);
841    typename MG::template EdgeMap<Num> residual_capacity(F);
842
843    while ( !bfs.finished() ) {
844      ResGWOutEdgeIt e=bfs;
845      if (res_graph.valid(e)) {
846        if (bfs.isBNodeNewlyReached()) {
847          dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
848          typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
849          original_edge.update();
850          original_edge.set(f, e);
851          residual_capacity.update();
852          residual_capacity.set(f, res_graph.resCap(e));
853        } else {
854          if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
855            typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
856            original_edge.update();
857            original_edge.set(f, e);
858            residual_capacity.update();
859            residual_capacity.set(f, res_graph.resCap(e));
860          }
861        }
862      }
863      ++bfs;
864    } //computing distances from s in the residual graph
865
866    bool __augment=true;
867
868    while (__augment) {
869      __augment=false;
870      //computing blocking flow with dfs
871      DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F);
872      typename MG::template NodeMap<typename MG::Edge> pred(F);
873      pred.set(sF, INVALID);
874      //invalid iterators for sources
875
876      typename MG::template NodeMap<Num> free(F);
877
878      dfs.pushAndSetReached(sF);     
879      while (!dfs.finished()) {
880        ++dfs;
881        if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
882          if (dfs.isBNodeNewlyReached()) {
883            typename MG::Node v=F.aNode(dfs);
884            typename MG::Node w=F.bNode(dfs);
885            pred.set(w, dfs);
886            if (F.valid(pred[v])) {
887              free.set(w, std::min(free[v], residual_capacity[dfs]));
888            } else {
889              free.set(w, residual_capacity[dfs]);
890            }
891            if (w==tF) {
892              __augment=true;
893              _augment=true;
894              break;
895            }
896             
897          } else {
898            F.erase(/*typename MG::OutEdgeIt*/(dfs));
899          }
900        }
901      }
902
903      if (__augment) {
904        typename MG::Node n=tF;
905        Num augment_value=free[tF];
906        while (F.valid(pred[n])) {
907          typename MG::Edge e=pred[n];
908          res_graph.augment(original_edge[e], augment_value);
909          n=F.tail(e);
910          if (residual_capacity[e]==augment_value)
911            F.erase(e);
912          else
913            residual_capacity.set(e, residual_capacity[e]-augment_value);
914        }
915      }
916       
917    }
918           
919    return _augment;
920  }
921
922
923
924
925
926
927  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
928  bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2()
929  {
930    bool _augment=false;
931
932    ResGW res_graph(*g, *capacity, *flow);
933     
934    //ReachedMap level(res_graph);
935    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
936    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
937
938    bfs.pushAndSetReached(s);
939    DistanceMap<ResGW> dist(res_graph);
940    while ( !bfs.finished() ) {
941      ResGWOutEdgeIt e=bfs;
942      if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
943        dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
944      }
945      ++bfs;
946    } //computing distances from s in the residual graph
947
948      //Subgraph containing the edges on some shortest paths
949    ConstMap<typename ResGW::Node, bool> true_map(true);
950    typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>,
951      DistanceMap<ResGW> > FilterResGW;
952    FilterResGW filter_res_graph(res_graph, true_map, dist);
953
954    //Subgraph, which is able to delete edges which are already
955    //met by the dfs
956    typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt>
957      first_out_edges(filter_res_graph);
958    typename FilterResGW::NodeIt v;
959    for(filter_res_graph.first(v); filter_res_graph.valid(v);
960        filter_res_graph.next(v))
961      {
962        typename FilterResGW::OutEdgeIt e;
963        filter_res_graph.first(e, v);
964        first_out_edges.set(v, e);
965      }
966    typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
967      template NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW;
968    ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);
969
970    bool __augment=true;
971
972    while (__augment) {
973
974      __augment=false;
975      //computing blocking flow with dfs
976      DfsIterator< ErasingResGW,
977        typename ErasingResGW::template NodeMap<bool> >
978        dfs(erasing_res_graph);
979      typename ErasingResGW::
980        template NodeMap<typename ErasingResGW::OutEdgeIt>
981        pred(erasing_res_graph);
982      pred.set(s, INVALID);
983      //invalid iterators for sources
984
985      typename ErasingResGW::template NodeMap<Num>
986        free1(erasing_res_graph);
987
988      dfs.pushAndSetReached(
989                            typename ErasingResGW::Node(
990                                                        typename FilterResGW::Node(
991                                                                                   typename ResGW::Node(s)
992                                                                                   )
993                                                        )
994                            );
995      while (!dfs.finished()) {
996        ++dfs;
997        if (erasing_res_graph.valid(
998                                    typename ErasingResGW::OutEdgeIt(dfs)))
999          {
1000            if (dfs.isBNodeNewlyReached()) {
1001         
1002              typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs);
1003              typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs);
1004
1005              pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs));
1006              if (erasing_res_graph.valid(pred[v])) {
1007                free1.set(w, std::min(free1[v], res_graph.resCap(
1008                                                                 typename ErasingResGW::OutEdgeIt(dfs))));
1009              } else {
1010                free1.set(w, res_graph.resCap(
1011                                              typename ErasingResGW::OutEdgeIt(dfs)));
1012              }
1013             
1014              if (w==t) {
1015                __augment=true;
1016                _augment=true;
1017                break;
1018              }
1019            } else {
1020              erasing_res_graph.erase(dfs);
1021            }
1022          }
1023      }
1024
1025      if (__augment) {
1026        typename ErasingResGW::Node n=typename FilterResGW::Node(typename ResGW::Node(t));
1027        //        typename ResGW::NodeMap<Num> a(res_graph);
1028        //        typename ResGW::Node b;
1029        //        Num j=a[b];
1030        //        typename FilterResGW::NodeMap<Num> a1(filter_res_graph);
1031        //        typename FilterResGW::Node b1;
1032        //        Num j1=a1[b1];
1033        //        typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);
1034        //        typename ErasingResGW::Node b2;
1035        //        Num j2=a2[b2];
1036        Num augment_value=free1[n];
1037        while (erasing_res_graph.valid(pred[n])) {
1038          typename ErasingResGW::OutEdgeIt e=pred[n];
1039          res_graph.augment(e, augment_value);
1040          n=erasing_res_graph.tail(e);
1041          if (res_graph.resCap(e)==0)
1042            erasing_res_graph.erase(e);
1043        }
1044      }
1045     
1046    } //while (__augment)
1047           
1048    return _augment;
1049  }
1050
1051
1052
1053
1054} //namespace hugo
1055
1056#endif //HUGO_MAX_FLOW_H
1057
1058
1059
1060
Note: See TracBrowser for help on using the repository browser.