COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/preflow_hl0.h @ 137:6364b07f8cd4

Last change on this file since 137:6364b07f8cd4 was 109:fc5982b39e10, checked in by jacint, 17 years ago

Flows with test files. The best is preflow.h

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