COIN-OR::LEMON - Graph Library

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