COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/preflow_hl4.h @ 109:fc5982b39e10

Last change on this file since 109:fc5982b39e10 was 109:fc5982b39e10, checked in by jacint, 20 years ago

Flows with test files. The best is preflow.h

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