COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/preflow_hl2.h @ 109:fc5982b39e10

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

Flows with test files. The best is preflow.h

File size: 8.0 KB
Line 
1// -*- C++ -*-
2/*
3preflow_hl2.h
4by jacint.
5Runs the highest label variant of the preflow push algorithm with
6running time O(n^2\sqrt(m)).
7
8Heuristics:
9
10  gap: we iterate through the nodes for finding the nodes above
11       the gap and under level n. So it is quite slow.
12  numb: we maintain the number of nodes in level i.
13  highest label
14 
15'A' is a parameter for the gap, we only upgrade the nodes to level n,
16  if the gap is under A*n.
17
18The constructor runs the algorithm.
19
20Members:
21
22T maxFlow() : returns the value of a maximum flow
23
24T flowOnEdge(EdgeIt e) : for a fixed maximum flow x it returns x(e)
25
26FlowMap Flow() : returns the fixed maximum flow x
27
28void minCut(CutMap& M) : sets M to the characteristic vector of a
29     minimum cut. M should be a map of bools initialized to false.
30
31void minMinCut(CutMap& M) : sets M to the characteristic vector of the
32     minimum min cut. M should be a map of bools initialized to false.
33
34void maxMinCut(CutMap& M) : sets M to the characteristic vector of the
35     maximum min cut. M should be a map of bools initialized to false.
36
37*/
38
39#ifndef PREFLOW_HL2_H
40#define PREFLOW_HL2_H
41
42#define A .9
43
44#include <vector>
45#include <stack>
46#include <queue>
47
48namespace hugo {
49
50  template <typename Graph, typename T,
51    typename FlowMap=typename Graph::EdgeMap<T>,
52    typename CapMap=typename Graph::EdgeMap<T> >
53  class preflow_hl2 {
54   
55    typedef typename Graph::NodeIt NodeIt;
56    typedef typename Graph::EdgeIt EdgeIt;
57    typedef typename Graph::EachNodeIt EachNodeIt;
58    typedef typename Graph::OutEdgeIt OutEdgeIt;
59    typedef typename Graph::InEdgeIt InEdgeIt;
60   
61    Graph& G;
62    NodeIt s;
63    NodeIt t;
64    FlowMap flow;
65    CapMap& capacity; 
66    T value;
67   
68  public:
69
70    preflow_hl2(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
71      G(_G), s(_s), t(_t), flow(_G), capacity(_capacity) {
72
73      int n=G.nodeNum();
74      int b=n-2;
75      /*
76        b is a bound on the highest level of an active node.
77      */
78
79      typename Graph::NodeMap<int> level(G,n);     
80      typename Graph::NodeMap<T> excess(G);
81
82      std::vector<int> numb(n);   
83      /*
84        The number of nodes on level i < n. It is
85        initialized to n+1, because of the reverse_bfs-part.
86      */
87
88      std::vector<std::stack<NodeIt> > stack(2*n-1);   
89      //Stack of the active nodes in level i.
90
91
92      /*Reverse_bfs from t, to find the starting level.*/
93      level.set(t,0);
94      std::queue<NodeIt> bfs_queue;
95      bfs_queue.push(t);
96
97      while (!bfs_queue.empty()) {
98
99        NodeIt v=bfs_queue.front();     
100        bfs_queue.pop();
101        int l=level.get(v)+1;
102
103        for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
104          NodeIt w=G.tail(e);
105          if ( level.get(w) == n ) {
106            bfs_queue.push(w);
107            ++numb[l];
108            level.set(w, l);
109          }
110        }
111      }
112       
113      level.set(s,n);
114
115
116
117      /* Starting flow. It is everywhere 0 at the moment. */     
118      for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e)
119        {
120          T c=capacity.get(e);
121          if ( c == 0 ) continue;
122          NodeIt w=G.head(e);
123          if ( w!=s ) {   
124            if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w);
125            flow.set(e, c);
126            excess.set(w, excess.get(w)+c);
127          }
128        }
129     
130      /*
131         End of preprocessing
132      */
133
134
135      /*
136        Push/relabel on the highest level active nodes.
137      */
138      /*While there exists an active node.*/
139      while (b) {
140        if ( stack[b].empty() ) {
141          --b;
142          continue;
143        }
144       
145        NodeIt w=stack[b].top();   
146        stack[b].pop();           
147        int lev=level.get(w);
148        T exc=excess.get(w);
149        int newlevel=2*n;      //In newlevel we bound the next level of w.
150       
151          for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
152           
153            if ( flow.get(e) == capacity.get(e) ) continue;
154            NodeIt v=G.head(e);           
155            //e=wv         
156           
157            if( lev > level.get(v) ) {     
158              /*Push is allowed now*/
159             
160              if ( excess.get(v)==0 && v != s && v !=t )
161                stack[level.get(v)].push(v);
162              /*v becomes active.*/
163             
164              T cap=capacity.get(e);
165              T flo=flow.get(e);
166              T remcap=cap-flo;
167             
168              if ( remcap >= exc ) {       
169                /*A nonsaturating push.*/
170               
171                flow.set(e, flo+exc);
172                excess.set(v, excess.get(v)+exc);
173                exc=0;
174                break;
175               
176              } else {
177                /*A saturating push.*/
178               
179                flow.set(e, cap );
180                excess.set(v, excess.get(v)+remcap);
181                exc-=remcap;
182              }
183            } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
184           
185          } //for out edges wv
186       
187       
188        if ( exc > 0 ) {       
189          for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
190           
191            if( flow.get(e) == 0 ) continue;
192            NodeIt v=G.tail(e); 
193            //e=vw
194           
195            if( lev > level.get(v) ) { 
196              /*Push is allowed now*/
197             
198              if ( excess.get(v)==0 && v != s && v !=t)
199                stack[level.get(v)].push(v);
200              /*v becomes active.*/
201             
202              T flo=flow.get(e);
203             
204              if ( flo >= exc ) {
205                /*A nonsaturating push.*/
206               
207                flow.set(e, flo-exc);
208                excess.set(v, excess.get(v)+exc);
209                exc=0;
210                break;
211              } else {                                               
212                /*A saturating push.*/
213               
214                excess.set(v, excess.get(v)+flo);
215                exc-=flo;
216                flow.set(e,0);
217              } 
218            } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
219           
220          } //for in edges vw
221         
222        } // if w still has excess after the out edge for cycle
223       
224        excess.set(w, exc);
225       
226
227       
228
229        /*
230          Relabel
231        */
232       
233        if ( exc > 0 ) {
234          //now 'lev' is the old level of w
235          level.set(w,++newlevel);
236         
237          if ( lev < n ) {
238            --numb[lev];
239           
240            if ( !numb[lev] && lev < A*n ) {  //If the level of w gets empty.
241             
242              for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) {
243                if (level.get(v) > lev && level.get(v) < n ) level.set(v,n); 
244              }
245              for (int i=lev+1 ; i!=n ; ++i) numb[i]=0;
246              if ( newlevel < n ) newlevel=n;
247            } else {
248              if ( newlevel < n ) ++numb[newlevel];
249            }
250          }
251         
252          stack[newlevel].push(w);
253          b=newlevel;
254         
255        }
256       
257      } // while(b)
258     
259     
260      value = excess.get(t);
261      /*Max flow value.*/
262
263
264    } //void run()
265
266
267
268
269
270    /*
271      Returns the maximum value of a flow.
272     */
273
274    T maxFlow() {
275      return value;
276    }
277
278
279
280    /*
281      For the maximum flow x found by the algorithm,
282      it returns the flow value on edge e, i.e. x(e).
283    */
284
285    T flowOnEdge(const EdgeIt e) {
286      return flow.get(e);
287    }
288
289
290
291    /*
292      Returns the maximum flow x found by the algorithm.
293    */
294
295    FlowMap Flow() {
296      return flow;
297    }
298
299
300
301
302    /*
303      Returns the minimum min cut, by a bfs from s in the residual graph.
304    */
305   
306    template<typename CutMap>
307    void minCut(CutMap& M) {
308   
309      std::queue<NodeIt> queue;
310     
311      M.set(s,true);     
312      queue.push(s);
313
314      while (!queue.empty()) {
315        NodeIt w=queue.front();
316        queue.pop();
317       
318        for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
319          NodeIt v=G.head(e);
320          if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
321            queue.push(v);
322            M.set(v, true);
323          }
324        }
325
326        for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
327          NodeIt v=G.tail(e);
328          if (!M.get(v) && flow.get(e) > 0 ) {
329            queue.push(v);
330            M.set(v, true);
331          }
332        }
333
334      }
335    }
336
337
338
339    /*
340      Returns the maximum min cut, by a reverse bfs
341      from t in the residual graph.
342    */
343   
344    template<typename CutMap>
345    void maxMinCut(CutMap& M) {
346   
347      std::queue<NodeIt> queue;
348     
349      M.set(t,true);       
350      queue.push(t);
351
352      while (!queue.empty()) {
353        NodeIt w=queue.front();
354        queue.pop();
355
356        for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
357          NodeIt v=G.tail(e);
358          if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
359            queue.push(v);
360            M.set(v, true);
361          }
362        }
363
364        for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
365          NodeIt v=G.head(e);
366          if (!M.get(v) && flow.get(e) > 0 ) {
367            queue.push(v);
368            M.set(v, true);
369          }
370        }
371      }
372
373      for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v) {
374        M.set(v, !M.get(v));
375      }
376
377    }
378
379
380
381    template<typename CutMap>
382    void minMinCut(CutMap& M) {
383      minCut(M);
384    }
385
386
387
388  };
389}//namespace marci
390#endif
391
392
393
394
Note: See TracBrowser for help on using the repository browser.