COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/preflow.h @ 389:770cc1f4861f

Last change on this file since 389:770cc1f4861f was 389:770cc1f4861f, checked in by marci, 16 years ago

modifications for better compatibility with gcc 3.4.0

File size: 13.3 KB
Line 
1// -*- C++ -*-
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
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
15 runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
16 
17Parameters H0 and H1 are initialized to 20 and 10.
18
19Constructors:
20
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
23
24Members:
25
26void run()
27
28T flowValue() : returns the value of a maximum flow
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
39FIXME reset
40
41*/
42
43#ifndef HUGO_PREFLOW_H
44#define HUGO_PREFLOW_H
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,
55            typename CapMap=typename Graph::template EdgeMap<T>,
56            typename FlowMap=typename Graph::template EdgeMap<T> >
57  class Preflow {
58   
59    typedef typename Graph::Node Node;
60    typedef typename Graph::Edge Edge;
61    typedef typename Graph::NodeIt NodeIt;
62    typedef typename Graph::OutEdgeIt OutEdgeIt;
63    typedef typename Graph::InEdgeIt InEdgeIt;
64   
65    const Graph& G;
66    Node s;
67    Node t;
68    const CapMap& capacity; 
69    FlowMap& flow;
70    T value;
71    bool constzero;
72
73  public:
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   
79    void run() {
80     
81      value=0;                //for the subsequent runs
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::template NodeMap<int> level(G,n);     
103      typename Graph::template NodeMap<T> excess(G);
104
105      std::vector<Node> active(n-1,INVALID);
106      typename Graph::template NodeMap<Node> next(G,INVALID);
107      //Stack of the active nodes in level i < n.
108      //We use it in both phases.
109
110      typename Graph::template NodeMap<Node> left(G,INVALID);
111      typename Graph::template NodeMap<Node> right(G,INVALID);
112      std::vector<Node> level_list(n,INVALID);
113      /*
114        List of the nodes in level i<n.
115      */
116
117
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            }
142          }
143        }
144
145        //the starting flow
146        OutEdgeIt e;
147        for(G.first(e,s); G.valid(e); G.next(e))
148        {
149          T c=capacity[e];
150          if ( c == 0 ) continue;
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;
156            }
157            flow.set(e, c);
158            excess.set(w, excess[w]+c);
159          }
160        }
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
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);
288            std::queue<Node> bfs_queue;
289            bfs_queue.push(s);
290           
291            while (!bfs_queue.empty()) {
292             
293              Node v=bfs_queue.front();
294              bfs_queue.pop();
295              int l=level[v]+1;
296             
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 ) {
302                  bfs_queue.push(u);
303                  level.set(u, l);
304                  if ( excess[u] > 0 ) {
305                    next.set(u,active[l]);
306                    active[l]=u;
307                  }
308                }
309              }
310           
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 ) {
316                  bfs_queue.push(u);
317                  level.set(u, l);
318                  if ( excess[u] > 0 ) {
319                    next.set(u,active[l]);
320                    active[l]=u;
321                  }
322                }
323              }
324            }
325            b=n-2;
326            }
327           
328        }
329         
330         
331        if ( !G.valid(active[b]) ) --b;
332        else {
333          end=false; 
334
335          Node w=active[b];
336          active[b]=next[w];
337          int lev=level[w];
338          T exc=excess[w];
339          int newlevel=n;       //bound on the next level of w
340         
341          OutEdgeIt e;
342          for(G.first(e,w); G.valid(e); G.next(e)) {
343           
344            if ( flow[e] == capacity[e] ) continue;
345            Node v=G.head(e);           
346            //e=wv         
347           
348            if( lev > level[v] ) {     
349              /*Push is allowed now*/
350             
351              if ( excess[v]==0 && v!=t && v!=s ) {
352                int lev_v=level[v];
353                next.set(v,active[lev_v]);
354                active[lev_v]=v;
355              }
356             
357              T cap=capacity[e];
358              T flo=flow[e];
359              T remcap=cap-flo;
360             
361              if ( remcap >= exc ) {       
362                /*A nonsaturating push.*/
363               
364                flow.set(e, flo+exc);
365                excess.set(v, excess[v]+exc);
366                exc=0;
367                break;
368               
369              } else {
370                /*A saturating push.*/
371               
372                flow.set(e, cap);
373                excess.set(v, excess[v]+remcap);
374                exc-=remcap;
375              }
376            } else if ( newlevel > level[v] ){
377              newlevel = level[v];
378            }       
379           
380          } //for out edges wv
381       
382       
383        if ( exc > 0 ) {       
384          InEdgeIt e;
385          for(G.first(e,w); G.valid(e); G.next(e)) {
386           
387            if( flow[e] == 0 ) continue;
388            Node v=G.tail(e); 
389            //e=vw
390           
391            if( lev > level[v] ) { 
392              /*Push is allowed now*/
393             
394              if ( excess[v]==0 && v!=t && v!=s ) {
395                int lev_v=level[v];
396                next.set(v,active[lev_v]);
397                active[lev_v]=v;
398              }
399             
400              T flo=flow[e];
401             
402              if ( flo >= exc ) {
403                /*A nonsaturating push.*/
404               
405                flow.set(e, flo-exc);
406                excess.set(v, excess[v]+exc);
407                exc=0;
408                break;
409              } else {                                               
410                /*A saturating push.*/
411               
412                excess.set(v, excess[v]+flo);
413                exc-=flo;
414                flow.set(e,0);
415              } 
416            } else if ( newlevel > level[v] ) {
417              newlevel = level[v];
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
440            Node right_n=right[w];
441            Node left_n=left[w];
442
443            if ( G.valid(right_n) ) {
444              if ( G.valid(left_n) ) {
445                right.set(left_n, right_n);
446                left.set(right_n, left_n);
447              } else {
448                level_list[lev]=right_n;   
449                left.set(right_n, INVALID);
450              }
451            } else {
452              if ( G.valid(left_n) ) {
453                right.set(left_n, INVALID);
454              } else {
455                level_list[lev]=INVALID;   
456              }
457            }
458            //unlacing ends
459               
460            if ( !G.valid(level_list[lev]) ) {
461             
462               //gapping starts
463              for (int i=lev; i!=k ; ) {
464                Node v=level_list[++i];
465                while ( G.valid(v) ) {
466                  level.set(v,n);
467                  v=right[v];
468                }
469                level_list[i]=INVALID;
470                if ( !what_heur ) active[i]=INVALID;
471              }     
472
473              level.set(w,n);
474              b=lev-1;
475              k=b;
476              //gapping ends
477           
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;
486                if ( k < newlevel ) ++k;      //now k=newlevel
487                Node first=level_list[newlevel];
488                if ( G.valid(first) ) left.set(first,w);
489                right.set(w,first);
490                left.set(w,INVALID);
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
520      value = excess[t];
521      /*Max flow value.*/
522     
523    } //void run()
524
525
526
527
528
529    /*
530      Returns the maximum value of a flow.
531     */
532
533    T flowValue() {
534      return value;
535    }
536
537
538    FlowMap Flow() {
539      return flow;
540      }
541
542
543   
544    void Flow(FlowMap& _flow ) {
545      NodeIt v;
546      for(G.first(v) ; G.valid(v); G.next(v))
547        _flow.set(v,flow[v]);
548    }
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   
559      std::queue<Node> queue;
560     
561      M.set(s,true);     
562      queue.push(s);
563
564      while (!queue.empty()) {
565        Node w=queue.front();
566        queue.pop();
567
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] ) {
572            queue.push(v);
573            M.set(v, true);
574          }
575        }
576
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 ) {
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   
598      std::queue<Node> queue;
599     
600      M.set(t,true);       
601      queue.push(t);
602
603      while (!queue.empty()) {
604        Node w=queue.front();
605        queue.pop();
606
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] ) {
612            queue.push(v);
613            M.set(v, true);
614          }
615        }
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 ) {
621            queue.push(v);
622            M.set(v, true);
623          }
624        }
625      }
626
627      NodeIt v;
628      for(G.first(v) ; G.valid(v); G.next(v)) {
629        M.set(v, !M[v]);
630      }
631
632    }
633
634
635
636    template<typename CutMap>
637    void minCut(CutMap& M) {
638      minMinCut(M);
639    }
640
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
655
656  };
657
658} //namespace hugo
659
660#endif //PREFLOW_H
661
662
663
664
Note: See TracBrowser for help on using the repository browser.