COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/preflow.h @ 137:6364b07f8cd4

Last change on this file since 137:6364b07f8cd4 was 113:cf7b01232d86, checked in by jacint, 17 years ago

* empty log message *

File size: 11.0 KB
Line 
1// -*- C++ -*-
2/*
3preflow.h
4by jacint.
5Heuristics:
6 2 phase
7 gap
8 list 'level_list' on the nodes on level i implemented by hand
9 stack 'active' on the active nodes on level i implemented by hand
10 runs heuristic 'highest label' for H1*n relabels
11 runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
12 
13Parameters H0 and H1 are initialized to 20 and 10.
14
15The best preflow I could ever write.
16
17The constructor runs the algorithm.
18
19Members:
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 Flow() : returns the fixed maximum flow x
26
27void minMinCut(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 maxMinCut(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
33void minCut(CutMap& M) : sets M to the characteristic vector of
34     a min cut. M should be a map of bools initialized to false.
35
36*/
37
38#ifndef PREFLOW_H
39#define PREFLOW_H
40
41#define H0 20
42#define H1 1
43
44#include <vector>
45#include <queue>
46
47#include <time_measure.h>
48
49namespace hugo {
50
51  template <typename Graph, typename T,
52    typename FlowMap=typename Graph::EdgeMap<T>,
53    typename CapMap=typename Graph::EdgeMap<T> >
54  class preflow {
55   
56    typedef typename Graph::NodeIt NodeIt;
57    typedef typename Graph::EdgeIt EdgeIt;
58    typedef typename Graph::EachNodeIt EachNodeIt;
59    typedef typename Graph::OutEdgeIt OutEdgeIt;
60    typedef typename Graph::InEdgeIt InEdgeIt;
61   
62    Graph& G;
63    NodeIt s;
64    NodeIt t;
65    FlowMap flow;
66    CapMap& capacity; 
67    T value;
68
69  public:
70    double time;
71    preflow(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity ) :
72      G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity)
73    {
74
75      bool phase=0;        //phase 0 is the 1st phase, phase 1 is the 2nd
76      int n=G.nodeNum();
77      int heur0=(int)(H0*n);  //time while running 'bound decrease'
78      int heur1=(int)(H1*n);  //time while running 'highest label'
79      int heur=heur1;         //starting time interval (#of relabels)
80      bool what_heur=1;       
81      /*
82        what_heur is 0 in case 'bound decrease'
83        and 1 in case 'highest label'
84      */
85      bool end=false;     
86      /*
87        Needed for 'bound decrease', 'true'
88        means no active nodes are above bound b.
89      */
90      int relabel=0;
91      int k=n-2;  //bound on the highest level under n containing a node
92      int b=k;    //bound on the highest level under n of an active node
93     
94      typename Graph::NodeMap<int> level(G,n);     
95      typename Graph::NodeMap<T> excess(G);
96
97      std::vector<NodeIt> active(n);
98      typename Graph::NodeMap<NodeIt> next(G);
99      //Stack of the active nodes in level i < n.
100      //We use it in both phases.
101
102      typename Graph::NodeMap<NodeIt> left(G);
103      typename Graph::NodeMap<NodeIt> right(G);
104      std::vector<NodeIt> level_list(n);
105      /*
106        List of the nodes in level i<n.
107      */
108
109      /*Reverse_bfs from t, to find the starting level.*/
110      level.set(t,0);
111      std::queue<NodeIt> bfs_queue;
112      bfs_queue.push(t);
113
114      while (!bfs_queue.empty()) {
115
116        NodeIt v=bfs_queue.front();     
117        bfs_queue.pop();
118        int l=level.get(v)+1;
119
120        for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
121          NodeIt w=G.tail(e);
122          if ( level.get(w) == n && w != s ) {
123            bfs_queue.push(w);
124            NodeIt first=level_list[l];
125            if ( first != 0 ) left.set(first,w);
126            right.set(w,first);
127            level_list[l]=w;
128            level.set(w, l);
129          }
130        }
131      }
132     
133      level.set(s,n);
134     
135
136      /* Starting flow. It is everywhere 0 at the moment. */     
137      for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e)
138        {
139          T c=capacity.get(e);
140          if ( c == 0 ) continue;
141          NodeIt w=G.head(e);
142          if ( level.get(w) < n ) {       
143            if ( excess.get(w) == 0 && w!=t ) {
144              next.set(w,active[level.get(w)]);
145              active[level.get(w)]=w;
146            }
147            flow.set(e, c);
148            excess.set(w, excess.get(w)+c);
149          }
150        }
151
152      /*
153         End of preprocessing
154      */
155
156
157
158      /*
159        Push/relabel on the highest level active nodes.
160      */       
161      while ( true ) {
162       
163        if ( b == 0 ) {
164          if ( phase ) break;
165         
166          if ( !what_heur && !end && k > 0 ) {
167            b=k;
168            end=true;
169          } else {
170            phase=1;
171            time=currTime();
172            level.set(s,0);
173            std::queue<NodeIt> bfs_queue;
174            bfs_queue.push(s);
175           
176            while (!bfs_queue.empty()) {
177             
178              NodeIt v=bfs_queue.front();       
179              bfs_queue.pop();
180              int l=level.get(v)+1;
181             
182              for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
183                if ( capacity.get(e) == flow.get(e) ) continue;
184                NodeIt u=G.tail(e);
185                if ( level.get(u) >= n ) {
186                  bfs_queue.push(u);
187                  level.set(u, l);
188                  if ( excess.get(u) > 0 ) {
189                    next.set(u,active[l]);
190                    active[l]=u;
191                  }
192                }
193              }
194           
195              for(OutEdgeIt e=G.template first<OutEdgeIt>(v); e.valid(); ++e) {
196                if ( 0 == flow.get(e) ) continue;
197                NodeIt u=G.head(e);
198                if ( level.get(u) >= n ) {
199                  bfs_queue.push(u);
200                  level.set(u, l);
201                  if ( excess.get(u) > 0 ) {
202                    next.set(u,active[l]);
203                    active[l]=u;
204                  }
205                }
206              }
207            }
208            b=n-2;
209            }
210           
211        }
212         
213         
214        if ( active[b] == 0 ) --b;
215        else {
216          end=false; 
217
218          NodeIt w=active[b];
219          active[b]=next.get(w);
220          int lev=level.get(w);
221          T exc=excess.get(w);
222          int newlevel=n;       //bound on the next level of w
223         
224          for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
225           
226            if ( flow.get(e) == capacity.get(e) ) continue;
227            NodeIt v=G.head(e);           
228            //e=wv         
229           
230            if( lev > level.get(v) ) {     
231              /*Push is allowed now*/
232             
233              if ( excess.get(v)==0 && v!=t && v!=s ) {
234                int lev_v=level.get(v);
235                next.set(v,active[lev_v]);
236                active[lev_v]=v;
237              }
238             
239              T cap=capacity.get(e);
240              T flo=flow.get(e);
241              T remcap=cap-flo;
242             
243              if ( remcap >= exc ) {       
244                /*A nonsaturating push.*/
245               
246                flow.set(e, flo+exc);
247                excess.set(v, excess.get(v)+exc);
248                exc=0;
249                break;
250               
251              } else {
252                /*A saturating push.*/
253               
254                flow.set(e, cap);
255                excess.set(v, excess.get(v)+remcap);
256                exc-=remcap;
257              }
258            } else if ( newlevel > level.get(v) ){
259              newlevel = level.get(v);
260            }       
261           
262          } //for out edges wv
263       
264       
265        if ( exc > 0 ) {       
266          for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
267           
268            if( flow.get(e) == 0 ) continue;
269            NodeIt v=G.tail(e); 
270            //e=vw
271           
272            if( lev > level.get(v) ) { 
273              /*Push is allowed now*/
274             
275              if ( excess.get(v)==0 && v!=t && v!=s ) {
276                int lev_v=level.get(v);
277                next.set(v,active[lev_v]);
278                active[lev_v]=v;
279              }
280             
281              T flo=flow.get(e);
282             
283              if ( flo >= exc ) {
284                /*A nonsaturating push.*/
285               
286                flow.set(e, flo-exc);
287                excess.set(v, excess.get(v)+exc);
288                exc=0;
289                break;
290              } else {                                               
291                /*A saturating push.*/
292               
293                excess.set(v, excess.get(v)+flo);
294                exc-=flo;
295                flow.set(e,0);
296              } 
297            } else if ( newlevel > level.get(v) ) {
298              newlevel = level.get(v);
299            }       
300          } //for in edges vw
301         
302        } // if w still has excess after the out edge for cycle
303       
304        excess.set(w, exc);
305         
306        /*
307          Relabel
308        */
309       
310
311        if ( exc > 0 ) {
312          //now 'lev' is the old level of w
313       
314          if ( phase ) {
315            level.set(w,++newlevel);
316            next.set(w,active[newlevel]);
317            active[newlevel]=w;
318            b=newlevel;
319          } else {
320            //unlacing starts
321            NodeIt right_n=right.get(w);
322            NodeIt left_n=left.get(w);
323
324            if ( right_n != 0 ) {
325              if ( left_n != 0 ) {
326                right.set(left_n, right_n);
327                left.set(right_n, left_n);
328              } else {
329                level_list[lev]=right_n;   
330                left.set(right_n, 0);
331              }
332            } else {
333              if ( left_n != 0 ) {
334                right.set(left_n, 0);
335              } else {
336                level_list[lev]=0;   
337
338              }
339            }
340            //unlacing ends
341               
342            //gapping starts
343            if ( level_list[lev]==0 ) {
344             
345              for (int i=lev; i!=k ; ) {
346                NodeIt v=level_list[++i];
347                while ( v != 0 ) {
348                  level.set(v,n);
349                  v=right.get(v);
350                }
351                level_list[i]=0;
352                if ( !what_heur ) active[i]=0;
353              }     
354
355              level.set(w,n);
356              b=lev-1;
357              k=b;
358              //gapping ends
359            } else {
360             
361              if ( newlevel == n ) level.set(w,n);
362              else {
363                level.set(w,++newlevel);
364                next.set(w,active[newlevel]);
365                active[newlevel]=w;
366                if ( what_heur ) b=newlevel;
367                if ( k < newlevel ) ++k;
368                NodeIt first=level_list[newlevel];
369                if ( first != 0 ) left.set(first,w);
370                right.set(w,first);
371                left.set(w,0);
372                level_list[newlevel]=w;
373              }
374            }
375
376
377            ++relabel;
378            if ( relabel >= heur ) {
379              relabel=0;
380              if ( what_heur ) {
381                what_heur=0;
382                heur=heur0;
383                end=false;
384              } else {
385                what_heur=1;
386                heur=heur1;
387                b=k;
388              }
389            }
390          } //phase 0
391         
392         
393        } // if ( exc > 0 )
394         
395       
396        }  // if stack[b] is nonempty
397       
398      } // while(true)
399
400
401      value = excess.get(t);
402      /*Max flow value.*/
403     
404    } //void run()
405
406
407
408
409
410    /*
411      Returns the maximum value of a flow.
412     */
413
414    T maxFlow() {
415      return value;
416    }
417
418
419
420    /*
421      For the maximum flow x found by the algorithm,
422      it returns the flow value on edge e, i.e. x(e).
423    */
424   
425    T flowOnEdge(EdgeIt e) {
426      return flow.get(e);
427    }
428
429
430
431    FlowMap Flow() {
432      return flow;
433      }
434
435
436   
437    void Flow(FlowMap& _flow ) {
438      for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v)
439        _flow.set(v,flow.get(v));
440        }
441
442
443
444    /*
445      Returns the minimum min cut, by a bfs from s in the residual graph.
446    */
447   
448    template<typename _CutMap>
449    void minMinCut(_CutMap& M) {
450   
451      std::queue<NodeIt> queue;
452     
453      M.set(s,true);     
454      queue.push(s);
455
456      while (!queue.empty()) {
457        NodeIt w=queue.front();
458        queue.pop();
459
460        for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
461          NodeIt v=G.head(e);
462          if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
463            queue.push(v);
464            M.set(v, true);
465          }
466        }
467
468        for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
469          NodeIt v=G.tail(e);
470          if (!M.get(v) && flow.get(e) > 0 ) {
471            queue.push(v);
472            M.set(v, true);
473          }
474        }
475      }
476    }
477
478
479 
480    /*
481      Returns the maximum min cut, by a reverse bfs
482      from t in the residual graph.
483    */
484   
485    template<typename _CutMap>
486    void maxMinCut(_CutMap& M) {
487   
488      std::queue<NodeIt> queue;
489     
490      M.set(t,true);       
491      queue.push(t);
492
493      while (!queue.empty()) {
494        NodeIt w=queue.front();
495        queue.pop();
496
497        for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
498          NodeIt v=G.tail(e);
499          if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
500            queue.push(v);
501            M.set(v, true);
502          }
503        }
504
505        for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
506          NodeIt v=G.head(e);
507          if (!M.get(v) && flow.get(e) > 0 ) {
508            queue.push(v);
509            M.set(v, true);
510          }
511        }
512      }
513
514      for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v) {
515        M.set(v, !M.get(v));
516      }
517
518    }
519
520
521
522    template<typename CutMap>
523    void minCut(CutMap& M) {
524      minMinCut(M);
525    }
526
527
528  };
529}//namespace marci
530#endif
531
532
533
534
Note: See TracBrowser for help on using the repository browser.