COIN-OR::LEMON - Graph Library

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