COIN-OR::LEMON - Graph Library

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