COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/max_flow.h @ 602:580b329c2a0c

Last change on this file since 602:580b329c2a0c was 602:580b329c2a0c, checked in by marci, 17 years ago

bfs_iterator -> bfs_dfs.h, some docs

File size: 26.9 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_dfs.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
368      std::queue<Node> bfs_queue;
369     
370      switch ( fe ) {
371      case ZERO_FLOW:
372        {
373          //Reverse_bfs from t, to find the starting level.
374          level.set(t,0);
375          bfs_queue.push(t);
376       
377          while (!bfs_queue.empty()) {
378           
379            Node v=bfs_queue.front();   
380            bfs_queue.pop();
381            int l=level[v]+1;
382           
383            InEdgeIt e;
384            for(g->first(e,v); g->valid(e); g->next(e)) {
385              Node w=g->tail(e);
386              if ( level[w] == n && w != s ) {
387                bfs_queue.push(w);
388                Node first=level_list[l];
389                if ( g->valid(first) ) left.set(first,w);
390                right.set(w,first);
391                level_list[l]=w;
392                level.set(w, l);
393              }
394            }
395          }
396         
397          //the starting flow
398          OutEdgeIt e;
399          for(g->first(e,s); g->valid(e); g->next(e))
400            {
401              Num c=(*capacity)[e];
402              if ( c <= 0 ) continue;
403              Node w=g->head(e);
404              if ( level[w] < n ) {       
405                if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
406                flow->set(e, c);
407                excess.set(w, excess[w]+c);
408              }
409            }
410          break;
411        }
412       
413      case GEN_FLOW:
414      case PREFLOW:
415        {
416          //Reverse_bfs from t in the residual graph,
417          //to find the starting level.
418          level.set(t,0);
419          bfs_queue.push(t);
420         
421          while (!bfs_queue.empty()) {
422           
423            Node v=bfs_queue.front();   
424            bfs_queue.pop();
425            int l=level[v]+1;
426           
427            InEdgeIt e;
428            for(g->first(e,v); g->valid(e); g->next(e)) {
429              if ( (*capacity)[e] <= (*flow)[e] ) continue;
430              Node w=g->tail(e);
431              if ( level[w] == n && w != s ) {
432                bfs_queue.push(w);
433                Node first=level_list[l];
434                if ( g->valid(first) ) left.set(first,w);
435                right.set(w,first);
436                level_list[l]=w;
437                level.set(w, l);
438              }
439            }
440           
441            OutEdgeIt f;
442            for(g->first(f,v); g->valid(f); g->next(f)) {
443              if ( 0 >= (*flow)[f] ) continue;
444              Node w=g->head(f);
445              if ( level[w] == n && w != s ) {
446                bfs_queue.push(w);
447                Node first=level_list[l];
448                if ( g->valid(first) ) left.set(first,w);
449                right.set(w,first);
450                level_list[l]=w;
451                level.set(w, l);
452              }
453            }
454          }
455         
456         
457          //the starting flow
458          OutEdgeIt e;
459          for(g->first(e,s); g->valid(e); g->next(e))
460            {
461              Num rem=(*capacity)[e]-(*flow)[e];
462              if ( rem <= 0 ) continue;
463              Node w=g->head(e);
464              if ( level[w] < n ) {       
465                if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
466                flow->set(e, (*capacity)[e]);
467                excess.set(w, excess[w]+rem);
468              }
469            }
470         
471          InEdgeIt f;
472          for(g->first(f,s); g->valid(f); g->next(f))
473            {
474              if ( (*flow)[f] <= 0 ) continue;
475              Node w=g->tail(f);
476              if ( level[w] < n ) {       
477                if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
478                excess.set(w, excess[w]+(*flow)[f]);
479                flow->set(f, 0);
480              }
481            } 
482          break;
483        } //case PREFLOW
484      }
485    } //preflowPreproc
486
487
488
489    void relabel(Node w, int newlevel, VecStack& active, 
490                 VecNode& level_list, NNMap& left,
491                 NNMap& right, int& b, int& k, bool what_heur )
492    {
493
494      Num lev=level[w];
495     
496      Node right_n=right[w];
497      Node left_n=left[w];
498     
499      //unlacing starts
500      if ( g->valid(right_n) ) {
501        if ( g->valid(left_n) ) {
502          right.set(left_n, right_n);
503          left.set(right_n, left_n);
504        } else {
505          level_list[lev]=right_n;   
506          left.set(right_n, INVALID);
507        }
508      } else {
509        if ( g->valid(left_n) ) {
510          right.set(left_n, INVALID);
511        } else {
512          level_list[lev]=INVALID;   
513        }
514      }
515      //unlacing ends
516               
517      if ( !g->valid(level_list[lev]) ) {
518             
519        //gapping starts
520        for (int i=lev; i!=k ; ) {
521          Node v=level_list[++i];
522          while ( g->valid(v) ) {
523            level.set(v,n);
524            v=right[v];
525          }
526          level_list[i]=INVALID;
527          if ( !what_heur ) {
528            while ( !active[i].empty() ) {
529              active[i].pop();    //FIXME: ezt szebben kene
530            }
531          }         
532        }
533       
534        level.set(w,n);
535        b=lev-1;
536        k=b;
537        //gapping ends
538       
539      } else {
540       
541        if ( newlevel == n ) level.set(w,n);
542        else {
543          level.set(w,++newlevel);
544          active[newlevel].push(w);
545          if ( what_heur ) b=newlevel;
546          if ( k < newlevel ) ++k;      //now k=newlevel
547          Node first=level_list[newlevel];
548          if ( g->valid(first) ) left.set(first,w);
549          right.set(w,first);
550          left.set(w,INVALID);
551          level_list[newlevel]=w;
552        }
553      }
554     
555    } //relabel
556
557
558    template<typename MapGraphWrapper>
559    class DistanceMap {
560    protected:
561      const MapGraphWrapper* g;
562      typename MapGraphWrapper::template NodeMap<int> dist;
563    public:
564      DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { }
565      void set(const typename MapGraphWrapper::Node& n, int a) {
566        dist.set(n, a);
567      }
568      int operator[](const typename MapGraphWrapper::Node& n)
569      { return dist[n]; }
570      //       int get(const typename MapGraphWrapper::Node& n) const {
571      //        return dist[n]; }
572      //       bool get(const typename MapGraphWrapper::Edge& e) const {
573      //        return (dist.get(g->tail(e))<dist.get(g->head(e))); }
574      bool operator[](const typename MapGraphWrapper::Edge& e) const {
575        return (dist[g->tail(e)]<dist[g->head(e)]);
576      }
577    };
578   
579  };
580
581
582  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
583  void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase0( flowEnum fe )
584  {
585     
586    int heur0=(int)(H0*n);  //time while running 'bound decrease'
587    int heur1=(int)(H1*n);  //time while running 'highest label'
588    int heur=heur1;         //starting time interval (#of relabels)
589    int numrelabel=0;
590     
591    bool what_heur=1;       
592    //It is 0 in case 'bound decrease' and 1 in case 'highest label'
593
594    bool end=false;     
595    //Needed for 'bound decrease', true means no active nodes are above bound b.
596
597    int k=n-2;  //bound on the highest level under n containing a node
598    int b=k;    //bound on the highest level under n of an active node
599     
600    VecStack active(n);
601     
602    NNMap left(*g, INVALID);
603    NNMap right(*g, INVALID);
604    VecNode level_list(n,INVALID);
605    //List of the nodes in level i<n, set to n.
606
607    NodeIt v;
608    for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
609    //setting each node to level n
610     
611    switch ( fe ) {
612    case PREFLOW:
613      {
614        //counting the excess
615        NodeIt v;
616        for(g->first(v); g->valid(v); g->next(v)) {
617          Num exc=0;
618         
619          InEdgeIt e;
620          for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
621          OutEdgeIt f;
622          for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
623           
624          excess.set(v,exc);     
625           
626          //putting the active nodes into the stack
627          int lev=level[v];
628          if ( exc > 0 && lev < n && v != t ) active[lev].push(v);
629        }
630        break;
631      }
632    case GEN_FLOW:
633      {
634        //Counting the excess of t
635        Num exc=0;
636         
637        InEdgeIt e;
638        for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
639        OutEdgeIt f;
640        for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
641         
642        excess.set(t,exc);     
643         
644        break;
645      }
646    default:
647      break;
648    }
649     
650    preflowPreproc( fe, active, level_list, left, right );
651    //End of preprocessing
652     
653     
654    //Push/relabel on the highest level active nodes.
655    while ( true ) {
656      if ( b == 0 ) {
657        if ( !what_heur && !end && k > 0 ) {
658          b=k;
659          end=true;
660        } else break;
661      }
662       
663      if ( active[b].empty() ) --b;
664      else {
665        end=false; 
666        Node w=active[b].top();
667        active[b].pop();
668        int newlevel=push(w,active);
669        if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list,
670                                     left, right, b, k, what_heur);
671         
672        ++numrelabel;
673        if ( numrelabel >= heur ) {
674          numrelabel=0;
675          if ( what_heur ) {
676            what_heur=0;
677            heur=heur0;
678            end=false;
679          } else {
680            what_heur=1;
681            heur=heur1;
682            b=k;
683          }
684        }
685      }
686    }
687  }
688
689
690
691  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
692  void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase1()
693  {
694     
695    int k=n-2;  //bound on the highest level under n containing a node
696    int b=k;    //bound on the highest level under n of an active node
697     
698    VecStack active(n);
699    level.set(s,0);
700    std::queue<Node> bfs_queue;
701    bfs_queue.push(s);
702           
703    while (!bfs_queue.empty()) {
704       
705      Node v=bfs_queue.front();
706      bfs_queue.pop();
707      int l=level[v]+1;
708             
709      InEdgeIt e;
710      for(g->first(e,v); g->valid(e); g->next(e)) {
711        if ( (*capacity)[e] <= (*flow)[e] ) continue;
712        Node u=g->tail(e);
713        if ( level[u] >= n ) {
714          bfs_queue.push(u);
715          level.set(u, l);
716          if ( excess[u] > 0 ) active[l].push(u);
717        }
718      }
719       
720      OutEdgeIt f;
721      for(g->first(f,v); g->valid(f); g->next(f)) {
722        if ( 0 >= (*flow)[f] ) continue;
723        Node u=g->head(f);
724        if ( level[u] >= n ) {
725          bfs_queue.push(u);
726          level.set(u, l);
727          if ( excess[u] > 0 ) active[l].push(u);
728        }
729      }
730    }
731    b=n-2;
732
733    while ( true ) {
734       
735      if ( b == 0 ) break;
736
737      if ( active[b].empty() ) --b;
738      else {
739        Node w=active[b].top();
740        active[b].pop();
741        int newlevel=push(w,active);     
742
743        //relabel
744        if ( excess[w] > 0 ) {
745          level.set(w,++newlevel);
746          active[newlevel].push(w);
747          b=newlevel;
748        }
749      }  // if stack[b] is nonempty
750    } // while(true)
751  }
752
753
754
755  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
756  bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath()
757  {
758    ResGW res_graph(*g, *capacity, *flow);
759    bool _augment=false;
760     
761    //ReachedMap level(res_graph);
762    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
763    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
764    bfs.pushAndSetReached(s);
765       
766    typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
767    pred.set(s, INVALID);
768     
769    typename ResGW::template NodeMap<Num> free(res_graph);
770       
771    //searching for augmenting path
772    while ( !bfs.finished() ) {
773      ResGWOutEdgeIt e=bfs;
774      if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
775        Node v=res_graph.tail(e);
776        Node w=res_graph.head(e);
777        pred.set(w, e);
778        if (res_graph.valid(pred[v])) {
779          free.set(w, std::min(free[v], res_graph.resCap(e)));
780        } else {
781          free.set(w, res_graph.resCap(e));
782        }
783        if (res_graph.head(e)==t) { _augment=true; break; }
784      }
785       
786      ++bfs;
787    } //end of searching augmenting path
788
789    if (_augment) {
790      Node n=t;
791      Num augment_value=free[t];
792      while (res_graph.valid(pred[n])) {
793        ResGWEdge e=pred[n];
794        res_graph.augment(e, augment_value);
795        n=res_graph.tail(e);
796      }
797    }
798
799    return _augment;
800  }
801
802
803
804
805
806
807
808
809
810  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
811  template<typename MutableGraph>
812  bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow()
813  {     
814    typedef MutableGraph MG;
815    bool _augment=false;
816
817    ResGW res_graph(*g, *capacity, *flow);
818
819    //bfs for distances on the residual graph
820    //ReachedMap level(res_graph);
821    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
822    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
823    bfs.pushAndSetReached(s);
824    typename ResGW::template NodeMap<int>
825      dist(res_graph); //filled up with 0's
826
827    //F will contain the physical copy of the residual graph
828    //with the set of edges which are on shortest paths
829    MG F;
830    typename ResGW::template NodeMap<typename MG::Node>
831      res_graph_to_F(res_graph);
832    {
833      typename ResGW::NodeIt n;
834      for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
835        res_graph_to_F.set(n, F.addNode());
836      }
837    }
838
839    typename MG::Node sF=res_graph_to_F[s];
840    typename MG::Node tF=res_graph_to_F[t];
841    typename MG::template EdgeMap<ResGWEdge> original_edge(F);
842    typename MG::template EdgeMap<Num> residual_capacity(F);
843
844    while ( !bfs.finished() ) {
845      ResGWOutEdgeIt e=bfs;
846      if (res_graph.valid(e)) {
847        if (bfs.isBNodeNewlyReached()) {
848          dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
849          typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
850          original_edge.update();
851          original_edge.set(f, e);
852          residual_capacity.update();
853          residual_capacity.set(f, res_graph.resCap(e));
854        } else {
855          if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
856            typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
857            original_edge.update();
858            original_edge.set(f, e);
859            residual_capacity.update();
860            residual_capacity.set(f, res_graph.resCap(e));
861          }
862        }
863      }
864      ++bfs;
865    } //computing distances from s in the residual graph
866
867    bool __augment=true;
868
869    while (__augment) {
870      __augment=false;
871      //computing blocking flow with dfs
872      DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F);
873      typename MG::template NodeMap<typename MG::Edge> pred(F);
874      pred.set(sF, INVALID);
875      //invalid iterators for sources
876
877      typename MG::template NodeMap<Num> free(F);
878
879      dfs.pushAndSetReached(sF);     
880      while (!dfs.finished()) {
881        ++dfs;
882        if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
883          if (dfs.isBNodeNewlyReached()) {
884            typename MG::Node v=F.aNode(dfs);
885            typename MG::Node w=F.bNode(dfs);
886            pred.set(w, dfs);
887            if (F.valid(pred[v])) {
888              free.set(w, std::min(free[v], residual_capacity[dfs]));
889            } else {
890              free.set(w, residual_capacity[dfs]);
891            }
892            if (w==tF) {
893              __augment=true;
894              _augment=true;
895              break;
896            }
897             
898          } else {
899            F.erase(/*typename MG::OutEdgeIt*/(dfs));
900          }
901        }
902      }
903
904      if (__augment) {
905        typename MG::Node n=tF;
906        Num augment_value=free[tF];
907        while (F.valid(pred[n])) {
908          typename MG::Edge e=pred[n];
909          res_graph.augment(original_edge[e], augment_value);
910          n=F.tail(e);
911          if (residual_capacity[e]==augment_value)
912            F.erase(e);
913          else
914            residual_capacity.set(e, residual_capacity[e]-augment_value);
915        }
916      }
917       
918    }
919           
920    return _augment;
921  }
922
923
924
925
926
927
928  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
929  bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2()
930  {
931    bool _augment=false;
932
933    ResGW res_graph(*g, *capacity, *flow);
934     
935    //ReachedMap level(res_graph);
936    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
937    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
938
939    bfs.pushAndSetReached(s);
940    DistanceMap<ResGW> dist(res_graph);
941    while ( !bfs.finished() ) {
942      ResGWOutEdgeIt e=bfs;
943      if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
944        dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
945      }
946      ++bfs;
947    } //computing distances from s in the residual graph
948
949      //Subgraph containing the edges on some shortest paths
950    ConstMap<typename ResGW::Node, bool> true_map(true);
951    typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>,
952      DistanceMap<ResGW> > FilterResGW;
953    FilterResGW filter_res_graph(res_graph, true_map, dist);
954
955    //Subgraph, which is able to delete edges which are already
956    //met by the dfs
957    typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt>
958      first_out_edges(filter_res_graph);
959    typename FilterResGW::NodeIt v;
960    for(filter_res_graph.first(v); filter_res_graph.valid(v);
961        filter_res_graph.next(v))
962      {
963        typename FilterResGW::OutEdgeIt e;
964        filter_res_graph.first(e, v);
965        first_out_edges.set(v, e);
966      }
967    typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
968      template NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW;
969    ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);
970
971    bool __augment=true;
972
973    while (__augment) {
974
975      __augment=false;
976      //computing blocking flow with dfs
977      DfsIterator< ErasingResGW,
978        typename ErasingResGW::template NodeMap<bool> >
979        dfs(erasing_res_graph);
980      typename ErasingResGW::
981        template NodeMap<typename ErasingResGW::OutEdgeIt>
982        pred(erasing_res_graph);
983      pred.set(s, INVALID);
984      //invalid iterators for sources
985
986      typename ErasingResGW::template NodeMap<Num>
987        free1(erasing_res_graph);
988
989      dfs.pushAndSetReached(
990                            typename ErasingResGW::Node(
991                                                        typename FilterResGW::Node(
992                                                                                   typename ResGW::Node(s)
993                                                                                   )
994                                                        )
995                            );
996      while (!dfs.finished()) {
997        ++dfs;
998        if (erasing_res_graph.valid(
999                                    typename ErasingResGW::OutEdgeIt(dfs)))
1000          {
1001            if (dfs.isBNodeNewlyReached()) {
1002         
1003              typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs);
1004              typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs);
1005
1006              pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs));
1007              if (erasing_res_graph.valid(pred[v])) {
1008                free1.set(w, std::min(free1[v], res_graph.resCap(
1009                                                                 typename ErasingResGW::OutEdgeIt(dfs))));
1010              } else {
1011                free1.set(w, res_graph.resCap(
1012                                              typename ErasingResGW::OutEdgeIt(dfs)));
1013              }
1014             
1015              if (w==t) {
1016                __augment=true;
1017                _augment=true;
1018                break;
1019              }
1020            } else {
1021              erasing_res_graph.erase(dfs);
1022            }
1023          }
1024      }
1025
1026      if (__augment) {
1027        typename ErasingResGW::Node n=typename FilterResGW::Node(typename ResGW::Node(t));
1028        //        typename ResGW::NodeMap<Num> a(res_graph);
1029        //        typename ResGW::Node b;
1030        //        Num j=a[b];
1031        //        typename FilterResGW::NodeMap<Num> a1(filter_res_graph);
1032        //        typename FilterResGW::Node b1;
1033        //        Num j1=a1[b1];
1034        //        typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);
1035        //        typename ErasingResGW::Node b2;
1036        //        Num j2=a2[b2];
1037        Num augment_value=free1[n];
1038        while (erasing_res_graph.valid(pred[n])) {
1039          typename ErasingResGW::OutEdgeIt e=pred[n];
1040          res_graph.augment(e, augment_value);
1041          n=erasing_res_graph.tail(e);
1042          if (res_graph.resCap(e)==0)
1043            erasing_res_graph.erase(e);
1044        }
1045      }
1046     
1047    } //while (__augment)
1048           
1049    return _augment;
1050  }
1051
1052
1053
1054
1055} //namespace hugo
1056
1057#endif //HUGO_MAX_FLOW_H
1058
1059
1060
1061
Note: See TracBrowser for help on using the repository browser.