COIN-OR::LEMON - Graph Library

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

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