COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/preflow.h @ 471:a40985a922d0

Last change on this file since 471:a40985a922d0 was 471:a40985a922d0, checked in by marci, 17 years ago

misc

File size: 15.1 KB
Line 
1// -*- C++ -*-
2
3/*
4Heuristics:
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 
12Parameters H0 and H1 are initialized to 20 and 1.
13
14Constructors:
15
16Preflow(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
19Members:
20
21void run()
22
23Num flowValue() : returns the value of a maximum flow
24
25void 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
28void 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
31void 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_PREFLOW_H
37#define HUGO_PREFLOW_H
38
39#define H0 20
40#define H1 1
41
42#include <vector>
43#include <queue>
44#include <stack>
45
46#include <graph_wrapper.h>
47
48namespace hugo {
49
50  template <typename Graph, typename Num,
51            typename CapMap=typename Graph::template EdgeMap<Num>,
52            typename FlowMap=typename Graph::template EdgeMap<Num> >
53  class Preflow {
54   
55    typedef typename Graph::Node Node;
56    typedef typename Graph::NodeIt NodeIt;
57    typedef typename Graph::OutEdgeIt OutEdgeIt;
58    typedef typename Graph::InEdgeIt InEdgeIt;
59
60    typedef typename std::vector<std::stack<Node> > VecStack;
61    typedef typename Graph::template NodeMap<Node> NNMap;
62    typedef typename std::vector<Node> VecNode;
63
64    const Graph* g;
65    Node s;
66    Node t;
67    const CapMap* capacity; 
68    FlowMap* flow;
69    int n;      //the number of nodes of G
70    typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
71    typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
72    typedef typename ResGW::Edge ResGWEdge;
73    //typedef typename ResGW::template NodeMap<bool> ReachedMap;
74    typedef typename Graph::template NodeMap<int> ReachedMap;
75    ReachedMap level;
76    //level works as a bool map in augmenting path algorithms
77    //and is used by bfs for storing reached information.
78    //In preflow, it shows levels of nodes.
79    //typename Graph::template NodeMap<int> level;   
80    typename Graph::template NodeMap<Num> excess;
81
82  public:
83 
84    enum flowEnum{
85      ZERO_FLOW=0,
86      GEN_FLOW=1,
87      PREFLOW=2
88    };
89
90    Preflow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,
91            FlowMap& _flow) :
92      g(&_G), s(_s), t(_t), capacity(&_capacity),
93      flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {}
94
95    void run() {
96      preflow( ZERO_FLOW );
97    }
98   
99    void preflow( flowEnum fe ) {
100      preflowPhase0(fe);
101      preflowPhase1();
102    }
103
104    void preflowPhase0( flowEnum fe );
105
106    void preflowPhase1();
107
108    bool augmentOnShortestPath();
109
110    template<typename MutableGraph> bool augmentOnBlockingFlow();
111
112    bool augmentOnBlockingFlow2();
113
114    //Returns the maximum value of a flow.
115    Num flowValue() {
116      return excess[t];
117    }
118
119    //should be used only between preflowPhase0 and preflowPhase1
120    template<typename _CutMap>
121    void actMinCut(_CutMap& M) {
122      NodeIt v;
123      for(g->first(v); g->valid(v); g->next(v))
124      if ( level[v] < n ) {
125        M.set(v,false);
126      } else {
127        M.set(v,true);
128      }
129    }
130
131
132
133    /*
134      Returns the minimum min cut, by a bfs from s in the residual graph.
135    */
136    template<typename _CutMap>
137    void minMinCut(_CutMap& M) {
138   
139      std::queue<Node> queue;
140     
141      M.set(s,true);     
142      queue.push(s);
143
144      while (!queue.empty()) {
145        Node w=queue.front();
146        queue.pop();
147
148        OutEdgeIt e;
149        for(g->first(e,w) ; g->valid(e); g->next(e)) {
150          Node v=g->head(e);
151          if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
152            queue.push(v);
153            M.set(v, true);
154          }
155        }
156
157        InEdgeIt f;
158        for(g->first(f,w) ; g->valid(f); g->next(f)) {
159          Node v=g->tail(f);
160          if (!M[v] && (*flow)[f] > 0 ) {
161            queue.push(v);
162            M.set(v, true);
163          }
164        }
165      }
166    }
167
168
169 
170    /*
171      Returns the maximum min cut, by a reverse bfs
172      from t in the residual graph.
173    */
174   
175    template<typename _CutMap>
176    void maxMinCut(_CutMap& M) {
177
178      NodeIt v;
179      for(g->first(v) ; g->valid(v); g->next(v)) {
180        M.set(v, true);
181      }
182
183      std::queue<Node> queue;
184     
185      M.set(t,false);       
186      queue.push(t);
187
188      while (!queue.empty()) {
189        Node w=queue.front();
190        queue.pop();
191
192
193        InEdgeIt e;
194        for(g->first(e,w) ; g->valid(e); g->next(e)) {
195          Node v=g->tail(e);
196          if (M[v] && (*flow)[e] < (*capacity)[e] ) {
197            queue.push(v);
198            M.set(v, false);
199          }
200        }
201       
202        OutEdgeIt f;
203        for(g->first(f,w) ; g->valid(f); g->next(f)) {
204          Node v=g->head(f);
205          if (M[v] && (*flow)[f] > 0 ) {
206            queue.push(v);
207            M.set(v, false);
208          }
209        }
210      }
211    }
212
213
214    template<typename CutMap>
215    void minCut(CutMap& M) {
216      minMinCut(M);
217    }
218
219    void resetTarget(Node _t) {t=_t;}
220    void resetSource(Node _s) {s=_s;}
221   
222    void resetCap(const CapMap& _cap) {
223      capacity=&_cap;
224    }
225   
226    void resetFlow(FlowMap& _flow) {
227      flow=&_flow;
228    }
229
230
231  private:
232
233    int push(Node w, VecStack& active) {
234     
235      int lev=level[w];
236      Num exc=excess[w];
237      int newlevel=n;       //bound on the next level of w
238         
239      OutEdgeIt e;
240      for(g->first(e,w); g->valid(e); g->next(e)) {
241           
242        if ( (*flow)[e] >= (*capacity)[e] ) continue;
243        Node v=g->head(e);           
244           
245        if( lev > level[v] ) { //Push is allowed now
246         
247          if ( excess[v]<=0 && v!=t && v!=s ) {
248            int lev_v=level[v];
249            active[lev_v].push(v);
250          }
251         
252          Num cap=(*capacity)[e];
253          Num flo=(*flow)[e];
254          Num remcap=cap-flo;
255         
256          if ( remcap >= exc ) { //A nonsaturating push.
257           
258            flow->set(e, flo+exc);
259            excess.set(v, excess[v]+exc);
260            exc=0;
261            break;
262           
263          } else { //A saturating push.
264            flow->set(e, cap);
265            excess.set(v, excess[v]+remcap);
266            exc-=remcap;
267          }
268        } else if ( newlevel > level[v] ) newlevel = level[v];
269      } //for out edges wv
270     
271      if ( exc > 0 ) { 
272        InEdgeIt e;
273        for(g->first(e,w); g->valid(e); g->next(e)) {
274         
275          if( (*flow)[e] <= 0 ) continue;
276          Node v=g->tail(e);
277         
278          if( lev > level[v] ) { //Push is allowed now
279           
280            if ( excess[v]<=0 && v!=t && v!=s ) {
281              int lev_v=level[v];
282              active[lev_v].push(v);
283            }
284           
285            Num flo=(*flow)[e];
286           
287            if ( flo >= exc ) { //A nonsaturating push.
288             
289              flow->set(e, flo-exc);
290              excess.set(v, excess[v]+exc);
291              exc=0;
292              break;
293            } else {  //A saturating push.
294             
295              excess.set(v, excess[v]+flo);
296              exc-=flo;
297              flow->set(e,0);
298            } 
299          } else if ( newlevel > level[v] ) newlevel = level[v];
300        } //for in edges vw
301       
302      } // if w still has excess after the out edge for cycle
303     
304      excess.set(w, exc);
305     
306      return newlevel;
307     }
308
309
310    void preflowPreproc ( flowEnum fe, VecStack& active,
311                          VecNode& level_list, NNMap& left, NNMap& right ) {
312
313      std::queue<Node> bfs_queue;
314     
315      switch ( fe ) {
316      case ZERO_FLOW:
317        {
318          //Reverse_bfs from t, to find the starting level.
319          level.set(t,0);
320          bfs_queue.push(t);
321       
322          while (!bfs_queue.empty()) {
323           
324            Node v=bfs_queue.front();   
325            bfs_queue.pop();
326            int l=level[v]+1;
327           
328            InEdgeIt e;
329            for(g->first(e,v); g->valid(e); g->next(e)) {
330              Node w=g->tail(e);
331              if ( level[w] == n && w != s ) {
332                bfs_queue.push(w);
333                Node first=level_list[l];
334                if ( g->valid(first) ) left.set(first,w);
335                right.set(w,first);
336                level_list[l]=w;
337                level.set(w, l);
338              }
339            }
340          }
341         
342          //the starting flow
343          OutEdgeIt e;
344          for(g->first(e,s); g->valid(e); g->next(e))
345            {
346              Num c=(*capacity)[e];
347              if ( c <= 0 ) continue;
348              Node w=g->head(e);
349              if ( level[w] < n ) {       
350                if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
351                flow->set(e, c);
352                excess.set(w, excess[w]+c);
353              }
354            }
355          break;
356        }
357       
358      case GEN_FLOW:
359      case PREFLOW:
360        {
361          //Reverse_bfs from t in the residual graph,
362          //to find the starting level.
363          level.set(t,0);
364          bfs_queue.push(t);
365         
366          while (!bfs_queue.empty()) {
367           
368            Node v=bfs_queue.front();   
369            bfs_queue.pop();
370            int l=level[v]+1;
371           
372            InEdgeIt e;
373            for(g->first(e,v); g->valid(e); g->next(e)) {
374              if ( (*capacity)[e] <= (*flow)[e] ) continue;
375              Node w=g->tail(e);
376              if ( level[w] == n && w != s ) {
377                bfs_queue.push(w);
378                Node first=level_list[l];
379                if ( g->valid(first) ) left.set(first,w);
380                right.set(w,first);
381                level_list[l]=w;
382                level.set(w, l);
383              }
384            }
385           
386            OutEdgeIt f;
387            for(g->first(f,v); g->valid(f); g->next(f)) {
388              if ( 0 >= (*flow)[f] ) continue;
389              Node w=g->head(f);
390              if ( level[w] == n && w != s ) {
391                bfs_queue.push(w);
392                Node first=level_list[l];
393                if ( g->valid(first) ) left.set(first,w);
394                right.set(w,first);
395                level_list[l]=w;
396                level.set(w, l);
397              }
398            }
399          }
400         
401         
402          //the starting flow
403          OutEdgeIt e;
404          for(g->first(e,s); g->valid(e); g->next(e))
405            {
406              Num rem=(*capacity)[e]-(*flow)[e];
407              if ( rem <= 0 ) continue;
408              Node w=g->head(e);
409              if ( level[w] < n ) {       
410                if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
411                flow->set(e, (*capacity)[e]);
412                excess.set(w, excess[w]+rem);
413              }
414            }
415         
416          InEdgeIt f;
417          for(g->first(f,s); g->valid(f); g->next(f))
418            {
419              if ( (*flow)[f] <= 0 ) continue;
420              Node w=g->tail(f);
421              if ( level[w] < n ) {       
422                if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
423                excess.set(w, excess[w]+(*flow)[f]);
424                flow->set(f, 0);
425              }
426            } 
427          break;
428        } //case PREFLOW
429      }
430    } //preflowPreproc
431
432
433
434    void relabel(Node w, int newlevel, VecStack& active, 
435                 VecNode& level_list, NNMap& left,
436                 NNMap& right, int& b, int& k, bool what_heur ) {
437
438      Num lev=level[w];
439     
440      Node right_n=right[w];
441      Node left_n=left[w];
442     
443      //unlacing starts
444      if ( g->valid(right_n) ) {
445        if ( g->valid(left_n) ) {
446          right.set(left_n, right_n);
447          left.set(right_n, left_n);
448        } else {
449          level_list[lev]=right_n;   
450          left.set(right_n, INVALID);
451        }
452      } else {
453        if ( g->valid(left_n) ) {
454          right.set(left_n, INVALID);
455        } else {
456          level_list[lev]=INVALID;   
457        }
458      }
459      //unlacing ends
460               
461      if ( !g->valid(level_list[lev]) ) {
462             
463        //gapping starts
464        for (int i=lev; i!=k ; ) {
465          Node v=level_list[++i];
466          while ( g->valid(v) ) {
467            level.set(v,n);
468            v=right[v];
469          }
470          level_list[i]=INVALID;
471          if ( !what_heur ) {
472            while ( !active[i].empty() ) {
473              active[i].pop();    //FIXME: ezt szebben kene
474            }
475          }         
476        }
477       
478        level.set(w,n);
479        b=lev-1;
480        k=b;
481        //gapping ends
482       
483      } else {
484       
485        if ( newlevel == n ) level.set(w,n);
486        else {
487          level.set(w,++newlevel);
488          active[newlevel].push(w);
489          if ( what_heur ) b=newlevel;
490          if ( k < newlevel ) ++k;      //now k=newlevel
491          Node first=level_list[newlevel];
492          if ( g->valid(first) ) left.set(first,w);
493          right.set(w,first);
494          left.set(w,INVALID);
495          level_list[newlevel]=w;
496        }
497      }
498     
499    } //relabel
500   
501  };
502
503
504  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
505  void Preflow<Graph, Num, CapMap, FlowMap>::preflowPhase0( flowEnum fe )
506  {
507     
508      int heur0=(int)(H0*n);  //time while running 'bound decrease'
509      int heur1=(int)(H1*n);  //time while running 'highest label'
510      int heur=heur1;         //starting time interval (#of relabels)
511      int numrelabel=0;
512     
513      bool what_heur=1;       
514      //It is 0 in case 'bound decrease' and 1 in case 'highest label'
515
516      bool end=false;     
517      //Needed for 'bound decrease', true means no active nodes are above bound b.
518
519      int k=n-2;  //bound on the highest level under n containing a node
520      int b=k;    //bound on the highest level under n of an active node
521     
522      VecStack active(n);
523     
524      NNMap left(*g, INVALID);
525      NNMap right(*g, INVALID);
526      VecNode level_list(n,INVALID);
527      //List of the nodes in level i<n, set to n.
528
529      NodeIt v;
530      for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
531      //setting each node to level n
532     
533      switch ( fe ) {
534      case PREFLOW:
535        {
536          //counting the excess
537          NodeIt v;
538          for(g->first(v); g->valid(v); g->next(v)) {
539            Num exc=0;
540         
541            InEdgeIt e;
542            for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
543            OutEdgeIt f;
544            for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
545           
546            excess.set(v,exc);   
547           
548            //putting the active nodes into the stack
549            int lev=level[v];
550            if ( exc > 0 && lev < n && v != t ) active[lev].push(v);
551          }
552          break;
553        }
554      case GEN_FLOW:
555        {
556          //Counting the excess of t
557          Num exc=0;
558         
559          InEdgeIt e;
560          for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
561          OutEdgeIt f;
562          for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
563         
564          excess.set(t,exc);   
565         
566          break;
567        }
568      default:
569        break;
570      }
571     
572      preflowPreproc( fe, active, level_list, left, right );
573      //End of preprocessing
574     
575     
576      //Push/relabel on the highest level active nodes.
577      while ( true ) {
578        if ( b == 0 ) {
579          if ( !what_heur && !end && k > 0 ) {
580            b=k;
581            end=true;
582          } else break;
583        }
584       
585        if ( active[b].empty() ) --b;
586        else {
587          end=false; 
588          Node w=active[b].top();
589          active[b].pop();
590          int newlevel=push(w,active);
591          if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list,
592                                       left, right, b, k, what_heur);
593         
594          ++numrelabel;
595          if ( numrelabel >= heur ) {
596            numrelabel=0;
597            if ( what_heur ) {
598              what_heur=0;
599              heur=heur0;
600              end=false;
601            } else {
602              what_heur=1;
603              heur=heur1;
604              b=k;
605            }
606          }
607        }
608      }
609    }
610
611
612
613  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
614  void Preflow<Graph, Num, CapMap, FlowMap>::preflowPhase1()
615  {
616     
617      int k=n-2;  //bound on the highest level under n containing a node
618      int b=k;    //bound on the highest level under n of an active node
619     
620      VecStack active(n);
621      level.set(s,0);
622      std::queue<Node> bfs_queue;
623      bfs_queue.push(s);
624           
625      while (!bfs_queue.empty()) {
626       
627        Node v=bfs_queue.front();       
628        bfs_queue.pop();
629        int l=level[v]+1;
630             
631        InEdgeIt e;
632        for(g->first(e,v); g->valid(e); g->next(e)) {
633          if ( (*capacity)[e] <= (*flow)[e] ) continue;
634          Node u=g->tail(e);
635          if ( level[u] >= n ) {
636            bfs_queue.push(u);
637            level.set(u, l);
638            if ( excess[u] > 0 ) active[l].push(u);
639          }
640        }
641       
642        OutEdgeIt f;
643        for(g->first(f,v); g->valid(f); g->next(f)) {
644          if ( 0 >= (*flow)[f] ) continue;
645          Node u=g->head(f);
646          if ( level[u] >= n ) {
647            bfs_queue.push(u);
648            level.set(u, l);
649            if ( excess[u] > 0 ) active[l].push(u);
650          }
651        }
652      }
653      b=n-2;
654
655      while ( true ) {
656       
657        if ( b == 0 ) break;
658
659        if ( active[b].empty() ) --b;
660        else {
661          Node w=active[b].top();
662          active[b].pop();
663          int newlevel=push(w,active);   
664
665          //relabel
666          if ( excess[w] > 0 ) {
667            level.set(w,++newlevel);
668            active[newlevel].push(w);
669            b=newlevel;
670          }
671        }  // if stack[b] is nonempty
672      } // while(true)
673    }
674
675
676
677
678
679
680
681
682
683
684
685
686} //namespace hugo
687
688#endif //HUGO_PREFLOW_H
689
690
691
692
Note: See TracBrowser for help on using the repository browser.