COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/preflow_aug.h @ 140:ca164520d31a

Last change on this file since 140:ca164520d31a was 140:ca164520d31a, checked in by jacint, 16 years ago

* empty log message *

File size: 12.3 KB
Line 
1// -*- C++ -*-
2
3//hogy kell azt megcsinalni hogy a konstruktorban a FlowMap-nek legyen default erteke (a 0 map)
4
5/*
6preflow_aug.h
7by jacint.
8
9The same as preflow.h, the only difference is that here
10we start from a given flow.
11
12
13The constructor runs the algorithm:
14
15template <typename Graph, typename T,
16    typename Graph_EdgeMap_T_1=typename Graph::EdgeMap<T>,
17    typename Graph_EdgeMap_T_2=typename Graph::EdgeMap<T> >
18PreflowAug(Graph G, G_NodeIt s, G_NodeIt t, G_EdgeMap_T_1 capacity,
19   G_EdgeMap_T_2 flow, bool flow_type)
20
21'capacity' must be non-negative, if flow_type=0 then
22'flow' must be a flow, otherwise it must be a preflow.
23
24
25Members:
26
27T maxFlow() : returns the value of a maximum flow
28
29T flow(G_EdgeIt e) : for a fixed maximum flow x it returns x(e)
30
31G_EdgeMap_T_2 flow() : returns the fixed maximum flow x
32
33void minMinCut(G_NodeMap_bool& M) :
34     sets M to the characteristic vector of the
35     minimum min cut. M must be initialized to false.
36
37void maxMinCut(G_NodeMap_bool& M) :
38     sets M to the characteristic vector of the
39     maximum min cut. M must be initialized to false.
40
41void minCut(G_NodeMap_bool& M) :
42     sets M to the characteristic vector of a 
43     min cut. M must be initialized to false.
44*/
45
46#ifndef PREFLOW_H
47#define PREFLOW_H
48
49#define H0 20
50#define H1 1
51
52#include <vector>
53#include <queue>
54
55#include <iostream> //for error handling
56
57#include <time_measure.h>
58
59namespace hugo {
60
61  template <typename Graph, typename T,
62    typename CapMap=typename Graph::EdgeMap<T>,
63    typename FlowMap=typename Graph::EdgeMap<T> >
64  class PreflowAug {
65   
66    typedef typename Graph::NodeIt NodeIt;
67    typedef typename Graph::EdgeIt EdgeIt;
68    typedef typename Graph::EachNodeIt EachNodeIt;
69    typedef typename Graph::EachEdgeIt EachEdgeIt;
70    typedef typename Graph::OutEdgeIt OutEdgeIt;
71    typedef typename Graph::InEdgeIt InEdgeIt;
72
73    Graph& G;
74    NodeIt s;
75    NodeIt t;
76    CapMap& capacity; 
77    FlowMap& _flow;
78    bool flow_type;
79    T value;
80
81  public:
82    double time;
83    PreflowAug(Graph& _G, NodeIt _s, NodeIt _t,
84               CapMap& _capacity, FlowMap& __flow, bool _flow_type ) :
85      G(_G), s(_s), t(_t), capacity(_capacity), _flow(__flow),
86      flow_type(_flow_type)
87    {
88
89
90      bool phase=0;        //phase 0 is the 1st phase, phase 1 is the 2nd
91      int n=G.nodeNum();
92      int heur0=(int)(H0*n);  //time while running 'bound decrease'
93      int heur1=(int)(H1*n);  //time while running 'highest label'
94      int heur=heur1;         //starting time interval (#of relabels)
95      bool what_heur=1;       
96      /*
97        what_heur is 0 in case 'bound decrease'
98        and 1 in case 'highest label'
99      */
100      bool end=false;     
101      /*
102        Needed for 'bound decrease', 'true'
103        means no active nodes are above bound b.
104      */
105      int relabel=0;
106      int k=n-2;  //bound on the highest level under n containing a node
107      int b=k;    //bound on the highest level under n of an active node
108     
109      typename Graph::NodeMap<int> level(G,n);     
110      typename Graph::NodeMap<T> excess(G);
111
112      std::vector<NodeIt> active(n);
113      typename Graph::NodeMap<NodeIt> next(G);
114      //Stack of the active nodes in level i < n.
115      //We use it in both phases.
116
117      typename Graph::NodeMap<NodeIt> left(G);
118      typename Graph::NodeMap<NodeIt> right(G);
119      std::vector<NodeIt> level_list(n);
120      /*
121        List of the nodes in level i<n.
122      */
123
124      /*Reverse_bfs from t, to find the starting level.*/
125      level.set(t,0);
126      std::queue<NodeIt> bfs_queue;
127      bfs_queue.push(t);
128
129
130      while (!bfs_queue.empty()) {
131       
132        NodeIt v=bfs_queue.front();     
133        bfs_queue.pop();
134        int l=level.get(v)+1;
135             
136        for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
137          if ( capacity.get(e) == _flow.get(e) ) continue;
138          NodeIt w=G.tail(e);
139          if ( level.get(w) == n && w != s ) {
140            bfs_queue.push(w);
141            NodeIt first=level_list[l];
142            if ( first.valid() ) left.set(first,w);
143            right.set(w,first);
144            level_list[l]=w;
145            level.set(w, l);
146          }
147        }
148           
149        for(OutEdgeIt e=G.template first<OutEdgeIt>(v); e.valid(); ++e) {
150          if ( 0 == _flow.get(e) ) continue;
151          NodeIt w=G.head(e);
152          if ( level.get(w) == n && w != s ) {
153            bfs_queue.push(w);
154            NodeIt first=level_list[l];
155            if ( first != 0 ) left.set(first,w);
156            right.set(w,first);
157            level_list[l]=w;
158            level.set(w, l);
159          }
160        }
161      }
162
163      level.set(s,n);
164
165
166
167      /*The starting excess.*/
168      if ( flow_type ) {
169        for(EachEdgeIt e=G.template first<EachEdgeIt>(); e.valid(); ++e) {
170          if ( _flow.get(e) > 0 ) {
171            T flo=_flow.get(e);
172            NodeIt u=G.tail(e);
173            NodeIt v=G.head(e);
174            excess.set(u, excess.get(u)-flo);
175            excess.set(v, excess.get(v)+flo);
176          }
177        }
178
179        for(EachNodeIt v=G.template first<EachNodeIt>(); v.valid(); ++v) {
180          if ( excess.get(v) < 0 ) {
181            std::cerr<<"It is not a pre_flow."<<std::endl;
182            exit(1);
183          } else {
184            next.set(v,active[level.get(v)]);
185            active[level.get(v)]=v;}
186        }
187      }
188
189
190      /* Starting flow. It is everywhere 0 at the moment. */     
191      for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e)
192        {
193          T c=capacity.get(e)-_flow.get(e);
194          if ( c == 0 ) continue;
195          NodeIt w=G.head(e);
196          if ( level.get(w) < n ) {       
197            if ( excess.get(w) == 0 && w!=t ) {
198              next.set(w,active[level.get(w)]);
199              active[level.get(w)]=w;
200            }
201            _flow.set(e, capacity.get(e));
202            excess.set(w, excess.get(w)+c);
203          }
204        }
205
206      /*
207         End of preprocessing
208      */
209
210
211
212      /*
213        Push/relabel on the highest level active nodes.
214      */       
215      while ( true ) {
216       
217        if ( b == 0 ) {
218          if ( phase ) break;
219         
220          if ( !what_heur && !end && k > 0 ) {
221            b=k;
222            end=true;
223          } else {
224            phase=1;
225            time=currTime();
226            level.set(s,0);
227            std::queue<NodeIt> bfs_queue;
228            bfs_queue.push(s);
229           
230            while (!bfs_queue.empty()) {
231             
232              NodeIt v=bfs_queue.front();       
233              bfs_queue.pop();
234              int l=level.get(v)+1;
235             
236              for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
237                if ( capacity.get(e) == _flow.get(e) ) continue;
238                NodeIt u=G.tail(e);
239                if ( level.get(u) >= n ) {
240                  bfs_queue.push(u);
241                  level.set(u, l);
242                  if ( excess.get(u) > 0 ) {
243                    next.set(u,active[l]);
244                    active[l]=u;
245                  }
246                }
247              }
248           
249              for(OutEdgeIt e=G.template first<OutEdgeIt>(v); e.valid(); ++e) {
250                if ( 0 == _flow.get(e) ) continue;
251                NodeIt u=G.head(e);
252                if ( level.get(u) >= n ) {
253                  bfs_queue.push(u);
254                  level.set(u, l);
255                  if ( excess.get(u) > 0 ) {
256                    next.set(u,active[l]);
257                    active[l]=u;
258                  }
259                }
260              }
261            }
262            b=n-2;
263            }
264           
265        }
266         
267         
268        if ( active[b] == 0 ) --b;
269        else {
270          end=false; 
271
272          NodeIt w=active[b];
273          active[b]=next.get(w);
274          int lev=level.get(w);
275          T exc=excess.get(w);
276          int newlevel=n;       //bound on the next level of w
277         
278          for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
279           
280            if ( _flow.get(e) == capacity.get(e) ) continue;
281            NodeIt v=G.head(e);           
282            //e=wv         
283           
284            if( lev > level.get(v) ) {     
285              /*Push is allowed now*/
286             
287              if ( excess.get(v)==0 && v!=t && v!=s ) {
288                int lev_v=level.get(v);
289                next.set(v,active[lev_v]);
290                active[lev_v]=v;
291              }
292             
293              T cap=capacity.get(e);
294              T flo=_flow.get(e);
295              T remcap=cap-flo;
296             
297              if ( remcap >= exc ) {       
298                /*A nonsaturating push.*/
299               
300                _flow.set(e, flo+exc);
301                excess.set(v, excess.get(v)+exc);
302                exc=0;
303                break;
304               
305              } else {
306                /*A saturating push.*/
307               
308                _flow.set(e, cap);
309                excess.set(v, excess.get(v)+remcap);
310                exc-=remcap;
311              }
312            } else if ( newlevel > level.get(v) ){
313              newlevel = level.get(v);
314            }       
315           
316          } //for out edges wv
317       
318       
319        if ( exc > 0 ) {       
320          for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
321           
322            if( _flow.get(e) == 0 ) continue;
323            NodeIt v=G.tail(e); 
324            //e=vw
325           
326            if( lev > level.get(v) ) { 
327              /*Push is allowed now*/
328             
329              if ( excess.get(v)==0 && v!=t && v!=s ) {
330                int lev_v=level.get(v);
331                next.set(v,active[lev_v]);
332                active[lev_v]=v;
333              }
334             
335              T flo=_flow.get(e);
336             
337              if ( flo >= exc ) {
338                /*A nonsaturating push.*/
339               
340                _flow.set(e, flo-exc);
341                excess.set(v, excess.get(v)+exc);
342                exc=0;
343                break;
344              } else {                                               
345                /*A saturating push.*/
346               
347                excess.set(v, excess.get(v)+flo);
348                exc-=flo;
349                _flow.set(e,0);
350              } 
351            } else if ( newlevel > level.get(v) ) {
352              newlevel = level.get(v);
353            }       
354          } //for in edges vw
355         
356        } // if w still has excess after the out edge for cycle
357       
358        excess.set(w, exc);
359         
360        /*
361          Relabel
362        */
363       
364
365        if ( exc > 0 ) {
366          //now 'lev' is the old level of w
367       
368          if ( phase ) {
369            level.set(w,++newlevel);
370            next.set(w,active[newlevel]);
371            active[newlevel]=w;
372            b=newlevel;
373          } else {
374            //unlacing starts
375            NodeIt right_n=right.get(w);
376            NodeIt left_n=left.get(w);
377
378            if ( right_n != 0 ) {
379              if ( left_n != 0 ) {
380                right.set(left_n, right_n);
381                left.set(right_n, left_n);
382              } else {
383                level_list[lev]=right_n;   
384                left.set(right_n, 0);
385              }
386            } else {
387              if ( left_n != 0 ) {
388                right.set(left_n, 0);
389              } else {
390                level_list[lev]=0;   
391
392              }
393            }
394            //unlacing ends
395               
396            //gapping starts
397            if ( level_list[lev]==0 ) {
398             
399              for (int i=lev; i!=k ; ) {
400                NodeIt v=level_list[++i];
401                while ( v != 0 ) {
402                  level.set(v,n);
403                  v=right.get(v);
404                }
405                level_list[i]=0;
406                if ( !what_heur ) active[i]=0;
407              }     
408
409              level.set(w,n);
410              b=lev-1;
411              k=b;
412              //gapping ends
413            } else {
414             
415              if ( newlevel == n ) level.set(w,n);
416              else {
417                level.set(w,++newlevel);
418                next.set(w,active[newlevel]);
419                active[newlevel]=w;
420                if ( what_heur ) b=newlevel;
421                if ( k < newlevel ) ++k;
422                NodeIt first=level_list[newlevel];
423                if ( first != 0 ) left.set(first,w);
424                right.set(w,first);
425                left.set(w,0);
426                level_list[newlevel]=w;
427              }
428            }
429
430
431            ++relabel;
432            if ( relabel >= heur ) {
433              relabel=0;
434              if ( what_heur ) {
435                what_heur=0;
436                heur=heur0;
437                end=false;
438              } else {
439                what_heur=1;
440                heur=heur1;
441                b=k;
442              }
443            }
444          } //phase 0
445         
446         
447        } // if ( exc > 0 )
448         
449       
450        }  // if stack[b] is nonempty
451       
452      } // while(true)
453
454
455      value = excess.get(t);
456      /*Max flow value.*/
457     
458    } //void run()
459
460
461
462
463
464    /*
465      Returns the maximum value of a flow.
466     */
467
468    T maxFlow() {
469      return value;
470    }
471
472
473
474    /*
475      For the maximum flow x found by the algorithm,
476      it returns the flow value on edge e, i.e. x(e).
477    */
478   
479    T flow(EdgeIt e) {
480      return _flow.get(e);
481    }
482
483
484
485    FlowMap flow() {
486      return _flow;
487      }
488
489
490   
491    void flow(FlowMap& __flow ) {
492      for(EachNodeIt v=G.template first<EachNodeIt>(); v.valid(); ++v)
493        __flow.set(v,_flow.get(v));
494        }
495
496
497
498    /*
499      Returns the minimum min cut, by a bfs from s in the residual graph.
500    */
501   
502    template<typename _CutMap>
503    void minMinCut(_CutMap& M) {
504   
505      std::queue<NodeIt> queue;
506     
507      M.set(s,true);     
508      queue.push(s);
509
510      while (!queue.empty()) {
511        NodeIt w=queue.front();
512        queue.pop();
513
514        for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
515          NodeIt v=G.head(e);
516          if (!M.get(v) && _flow.get(e) < capacity.get(e) ) {
517            queue.push(v);
518            M.set(v, true);
519          }
520        }
521
522        for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
523          NodeIt v=G.tail(e);
524          if (!M.get(v) && _flow.get(e) > 0 ) {
525            queue.push(v);
526            M.set(v, true);
527          }
528        }
529      }
530    }
531
532
533 
534    /*
535      Returns the maximum min cut, by a reverse bfs
536      from t in the residual graph.
537    */
538   
539    template<typename _CutMap>
540    void maxMinCut(_CutMap& M) {
541   
542      std::queue<NodeIt> queue;
543     
544      M.set(t,true);       
545      queue.push(t);
546
547      while (!queue.empty()) {
548        NodeIt w=queue.front();
549        queue.pop();
550
551        for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
552          NodeIt v=G.tail(e);
553          if (!M.get(v) && _flow.get(e) < capacity.get(e) ) {
554            queue.push(v);
555            M.set(v, true);
556          }
557        }
558
559        for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
560          NodeIt v=G.head(e);
561          if (!M.get(v) && _flow.get(e) > 0 ) {
562            queue.push(v);
563            M.set(v, true);
564          }
565        }
566      }
567
568      for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v) {
569        M.set(v, !M.get(v));
570      }
571
572    }
573
574
575
576    template<typename CutMap>
577    void minCut(CutMap& M) {
578      minMinCut(M);
579    }
580
581
582  };
583}//namespace marci
584#endif
585
586
587
588
Note: See TracBrowser for help on using the repository browser.