COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/preflow_excess.h @ 493:bbd1db03f0fe

Last change on this file since 493:bbd1db03f0fe was 437:9853b743d830, checked in by jacint, 20 years ago

Testing preprocess.

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