COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/preflow_hl2.h @ 105:a3c73e9b9b2e

Last change on this file since 105:a3c73e9b9b2e was 105:a3c73e9b9b2e, checked in by Alpar Juttner, 16 years ago

marci -> hugo replacements
resize -> update replacements

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