COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/preflowproba.h @ 374:0fc9cd9b854a

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