COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/preflow.h @ 372:e6a156fc186d

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