COIN-OR::LEMON - Graph Library

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

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

marci -> hugo replacements
resize -> update replacements

File size: 9.7 KB
RevLine 
[101]1// -*- C++ -*-
2/*
3preflow_hl3.h
4by jacint.
5Runs the highest label variant of the preflow push algorithm with
6running time O(n^2\sqrt(m)), with the felszippantos 'empty level'
7and with the two-phase heuristic: if there is no active node of
8level at most n, then we go into phase 1, do a bfs
9from s, and flow the excess back to s.
10
11In phase 1 we shift everything downwards by n.
12
13'A' is a parameter for the empty_level heuristic
14
15Member functions:
16
17void run() : runs the algorithm
18
19 The following functions should be used after run() was already run.
20
21T maxflow() : returns the value of a maximum flow
22
23T flowonedge(EdgeIt e) : for a fixed maximum flow x it returns x(e)
24
25FlowMap allflow() : returns the fixed maximum flow x
26
27void mincut(CutMap& M) : sets M to the characteristic vector of a
28     minimum cut. M should be a map of bools initialized to false.
29
30void min_mincut(CutMap& M) : sets M to the characteristic vector of the
31     minimum min cut. M should be a map of bools initialized to false.
32
33void max_mincut(CutMap& M) : sets M to the characteristic vector of the
34     maximum min cut. M should be a map of bools initialized to false.
35
36*/
37
38#ifndef PREFLOW_HL3_H
39#define PREFLOW_HL3_H
40
41#define A 1
42
43#include <vector>
44#include <stack>
45#include <queue>
46
[105]47namespace hugo {
[101]48
49  template <typename Graph, typename T,
50    typename FlowMap=typename Graph::EdgeMap<T>, typename CapMap=typename Graph::EdgeMap<T>,
51    typename IntMap=typename Graph::NodeMap<int>, typename TMap=typename Graph::NodeMap<T> >
52  class preflow_hl3 {
53   
54    typedef typename Graph::NodeIt NodeIt;
55    typedef typename Graph::EdgeIt EdgeIt;
56    typedef typename Graph::EachNodeIt EachNodeIt;
57    typedef typename Graph::OutEdgeIt OutEdgeIt;
58    typedef typename Graph::InEdgeIt InEdgeIt;
59   
60    Graph& G;
61    NodeIt s;
62    NodeIt t;
63    FlowMap flow;
64    CapMap& capacity; 
65    T value;
66   
67  public:
68
69    preflow_hl3(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
70      G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity) { }
71
72
73    void run() {
74 
75      bool phase=0;
76      int n=G.nodeNum();
77      int b=n-2;
78      /*
79        b is a bound on the highest level of the stack.
80        In the beginning it is at most n-2.
81      */
82
83      IntMap level(G,n);     
84      TMap excess(G);
85     
86      std::vector<int> numb(n);   
87      /*
88        The number of nodes on level i < n. It is
89        initialized to n+1, because of the reverse_bfs-part.
90        Needed only in phase 0.
91      */
92
93      std::vector<std::stack<NodeIt> > stack(n);   
94      //Stack of the active nodes in level i < n.
95      //We use it in both phases.
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 ) {
112            bfs_queue.push(w);
113            ++numb[l];
114            level.set(w, l);
115          }
116        }
117      }
118     
119      level.set(s,n);
120
121
122
123      /* Starting flow. It is everywhere 0 at the moment. */     
124      for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e)
125        {
126          T c=capacity.get(e);
127          if ( c == 0 ) continue;
128          NodeIt w=G.head(e);
129          if ( level.get(w) < n ) {       
130            if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w);
131            flow.set(e, c);
132            excess.set(w, excess.get(w)+c);
133          }
134        }
135
136      /*
137         End of preprocessing
138      */
139
140
141
142      /*
143        Push/relabel on the highest level active nodes.
144      */       
145      /*While there exists an active node.*/
146      while ( true ) {
147
148        if ( b == 0 ) {
149          if ( phase ) break;
150          phase=1;
151
152          level.set(s,0);
153
154          std::queue<NodeIt> bfs_queue;
155          bfs_queue.push(s);
156         
157          while (!bfs_queue.empty()) {
158           
159            NodeIt v=bfs_queue.front();
160            bfs_queue.pop();
161            int l=level.get(v)+1;
162
163            for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
164              if ( capacity.get(e) == flow.get(e) ) continue;
165              NodeIt u=G.tail(e);
166              if ( level.get(u) == n ) {
167                bfs_queue.push(u);
168                level.set(u, l);
169                if ( excess.get(u) > 0 ) stack[l].push(u);
170              }
171            }
172
173            for(OutEdgeIt e=G.template first<OutEdgeIt>(v); e.valid(); ++e) {
174              if ( 0 == flow.get(e) ) continue;
175              NodeIt u=G.head(e);
176              if ( level.get(u) == n ) {
177                bfs_queue.push(u);
178                level.set(u, l);
179                if ( excess.get(u) > 0 ) stack[l].push(u);
180              }
181            }
182          }
183
184          b=n-2;
185        }
186
187        if ( stack[b].empty() ) --b;
188        else {
189         
190          NodeIt w=stack[b].top();        //w is a highest label active node.
191          stack[b].pop();           
192          int lev=level.get(w);
193          int exc=excess.get(w);
194          int newlevel=n;                 //In newlevel we bound the next level of w.
195         
196          for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
197           
198            if ( flow.get(e) == capacity.get(e) ) continue;
199            NodeIt v=G.head(e);           
200            //e=wv         
201           
202            if( lev > level.get(v) ) {     
203              /*Push is allowed now*/
204             
205              if ( excess.get(v)==0 && v !=t && v!=s )
206                stack[level.get(v)].push(v);
207              /*v becomes active.*/
208             
209              int cap=capacity.get(e);
210              int flo=flow.get(e);
211              int remcap=cap-flo;
212             
213              if ( remcap >= exc ) {       
214                /*A nonsaturating push.*/
215               
216                flow.set(e, flo+exc);
217                excess.set(v, excess.get(v)+exc);
218                exc=0;
219                break;
220               
221              } else {
222                /*A saturating push.*/
223               
224                flow.set(e, cap );
225                excess.set(v, excess.get(v)+remcap);
226                exc-=remcap;
227              }
228            } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
229           
230          } //for out edges wv
231       
232       
233        if ( exc > 0 ) {       
234          for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
235           
236            if( flow.get(e) == 0 ) continue;
237            NodeIt v=G.tail(e); 
238            //e=vw
239           
240            if( lev > level.get(v) ) { 
241              /*Push is allowed now*/
242             
243              if ( excess.get(v)==0 && v !=t && v!=s )
244                stack[level.get(v)].push(v);
245              /*v becomes active.*/
246             
247              int flo=flow.get(e);
248             
249              if ( flo >= exc ) {
250                /*A nonsaturating push.*/
251               
252                flow.set(e, flo-exc);
253                excess.set(v, excess.get(v)+exc);
254                exc=0;
255                break;
256              } else {                                               
257                /*A saturating push.*/
258               
259                excess.set(v, excess.get(v)+flo);
260                exc-=flo;
261                flow.set(e,0);
262              } 
263            } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
264           
265          } //for in edges vw
266         
267        } // if w still has excess after the out edge for cycle
268       
269        excess.set(w, exc);
270         
271
272        /*
273          Relabel
274        */
275       
276        if ( exc > 0 ) {
277          //now 'lev' is the old level of w
278       
279          if ( phase ) {
280            level.set(w,++newlevel);
281            stack[newlevel].push(w);
282            b=newlevel;
283          } else {
284
285            if ( newlevel >= n-2 || --numb[lev] == 0 ) {
286             
287              level.set(w,n);
288             
289              if ( newlevel < n ) {
290               
291                std::queue<NodeIt> bfs_queue;
292                bfs_queue.push(w);
293
294                while (!bfs_queue.empty()) {
295
296                  NodeIt v=bfs_queue.front();   
297                  bfs_queue.pop();
298
299                  for(OutEdgeIt e=G.template first<OutEdgeIt>(v); e.valid(); ++e) {
300                    if ( capacity.get(e) == flow.get(e) ) continue;
301                    NodeIt u=G.head(e);
302                    if ( level.get(u) < n ) {
303                      bfs_queue.push(u);
304                      --numb[level.get(u)];
305                      level.set(u, n);
306                    }
307                  }
308
309                  for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
310                    if ( 0 == flow.get(e) ) continue;
311                    NodeIt u=G.tail(e);
312                    if ( level.get(u) < n ) {
313                      bfs_queue.push(u);
314                      --numb[level.get(u)];
315                      level.set(u, n);
316                    }
317                  }
318                }
319              }
320              b=n-1;
321
322            } else {
323              level.set(w,++newlevel);
324              stack[newlevel].push(w);
325              ++numb[newlevel];
326              b=newlevel;
327            }
328          }
329        }
330
331 
332       
333        } // if stack[b] is nonempty
334       
335      } // while(true)
336
337
338      value = excess.get(t);
339      /*Max flow value.*/
340
341
342    } //void run()
343
344
345
346
347
348    /*
349      Returns the maximum value of a flow.
350     */
351
352    T maxflow() {
353      return value;
354    }
355
356
357
358    /*
359      For the maximum flow x found by the algorithm, it returns the flow value on Edge e, i.e. x(e).
360    */
361
362    T flowonedge(EdgeIt e) {
363      return flow.get(e);
364    }
365
366
367
368    /*
369      Returns the maximum flow x found by the algorithm.
370    */
371
372    FlowMap allflow() {
373      return flow;
374    }
375
376
377
378
379    /*
380      Returns the minimum min cut, by a bfs from s in the residual graph.
381    */
382   
383    template<typename CutMap>
384    void mincut(CutMap& M) {
385   
386      std::queue<NodeIt> queue;
387     
388      M.set(s,true);     
389      queue.push(s);
390
391      while (!queue.empty()) {
392        NodeIt w=queue.front();
393        queue.pop();
394
395        for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
396          NodeIt v=G.head(e);
397          if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
398            queue.push(v);
399            M.set(v, true);
400          }
401        }
402
403        for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
404          NodeIt v=G.tail(e);
405          if (!M.get(v) && flow.get(e) > 0 ) {
406            queue.push(v);
407            M.set(v, true);
408          }
409        }
410
411      }
412
413    }
414
415
416
417    /*
418      Returns the maximum min cut, by a reverse bfs
419      from t in the residual graph.
420    */
421   
422    template<typename CutMap>
423    void max_mincut(CutMap& M) {
424   
425      std::queue<NodeIt> queue;
426     
427      M.set(t,true);       
428      queue.push(t);
429
430      while (!queue.empty()) {
431        NodeIt w=queue.front();
432        queue.pop();
433
434        for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
435          NodeIt v=G.tail(e);
436          if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
437            queue.push(v);
438            M.set(v, true);
439          }
440        }
441
442        for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
443          NodeIt v=G.head(e);
444          if (!M.get(v) && flow.get(e) > 0 ) {
445            queue.push(v);
446            M.set(v, true);
447          }
448        }
449      }
450
451      for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v) {
452        M.set(v, !M.get(v));
453      }
454
455    }
456
457
458
459    template<typename CutMap>
460    void min_mincut(CutMap& M) {
461      mincut(M);
462    }
463
464
465
466  };
[105]467}//namespace hugo
[101]468#endif
469
470
471
472
Note: See TracBrowser for help on using the repository browser.