COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/preflow_hl4.h @ 108:0351b00fd283

Last change on this file since 108:0351b00fd283 was 105:a3c73e9b9b2e, checked in by Alpar Juttner, 21 years ago

marci -> hugo replacements
resize -> update replacements

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