COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/preflow.h @ 451:6b36be4cffa4

Last change on this file since 451:6b36be4cffa4 was 451:6b36be4cffa4, checked in by jacint, 16 years ago

Changes in the interface and new test program added.

File size: 14.0 KB
Line 
1// -*- C++ -*-
2
3/*
4Heuristics:
5 2 phase
6 gap
7 list 'level_list' on the nodes on level i implemented by hand
8 stack 'active' on the active nodes on level i
9 runs heuristic 'highest label' for H1*n relabels
10 runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
11 
12Parameters H0 and H1 are initialized to 20 and 1.
13
14Constructors:
15
16Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if
17     FlowMap is not constant zero, and should be true if it is
18
19Members:
20
21void run()
22
23T flowValue() : returns the value of a maximum flow
24
25void minMinCut(CutMap& M) : sets M to the characteristic vector of the
26     minimum min cut. M should be a map of bools initialized to false. ??Is it OK?
27
28void maxMinCut(CutMap& M) : sets M to the characteristic vector of the
29     maximum min cut. M should be a map of bools initialized to false.
30
31void minCut(CutMap& M) : sets M to the characteristic vector of
32     a min cut. M should be a map of bools initialized to false.
33
34*/
35
36#ifndef HUGO_PREFLOW_H
37#define HUGO_PREFLOW_H
38
39#define H0 20
40#define H1 1
41
42#include <vector>
43#include <queue>
44#include <stack>
45
46namespace hugo {
47
48  template <typename Graph, typename T,
49            typename CapMap=typename Graph::template EdgeMap<T>,
50            typename FlowMap=typename Graph::template EdgeMap<T> >
51  class Preflow {
52   
53    typedef typename Graph::Node Node;
54    typedef typename Graph::NodeIt NodeIt;
55    typedef typename Graph::OutEdgeIt OutEdgeIt;
56    typedef typename Graph::InEdgeIt InEdgeIt;
57
58    typedef typename std::vector<std::stack<Node> > VecStack;
59    typedef typename Graph::template NodeMap<Node> NNMap;
60    typedef typename std::vector<Node> VecNode;
61
62    const Graph& G;
63    Node s;
64    Node t;
65    CapMap* capacity; 
66    FlowMap* flow;
67    int n;      //the number of nodes of G
68    typename Graph::template NodeMap<int> level;     
69    typename Graph::template NodeMap<T> excess;
70
71
72  public:
73 
74    enum flowEnum{
75      ZERO_FLOW=0,
76      GEN_FLOW=1,
77      PREFLOW=2
78    };
79
80    Preflow(Graph& _G, Node _s, Node _t, CapMap& _capacity,
81            FlowMap& _flow) :
82      G(_G), s(_s), t(_t), capacity(&_capacity),
83      flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {}
84
85    void run() {
86      preflow( ZERO_FLOW );
87    }
88   
89    void preflow( flowEnum fe ) {
90      preflowPhase0(fe);
91      preflowPhase1();
92    }
93
94    void preflowPhase0( flowEnum fe ) {
95     
96      int heur0=(int)(H0*n);  //time while running 'bound decrease'
97      int heur1=(int)(H1*n);  //time while running 'highest label'
98      int heur=heur1;         //starting time interval (#of relabels)
99      int numrelabel=0;
100     
101      bool what_heur=1;       
102      //It is 0 in case 'bound decrease' and 1 in case 'highest label'
103
104      bool end=false;     
105      //Needed for 'bound decrease', true means no active nodes are above bound b.
106
107      int k=n-2;  //bound on the highest level under n containing a node
108      int b=k;    //bound on the highest level under n of an active node
109     
110      VecStack active(n);
111     
112      NNMap left(G,INVALID);
113      NNMap right(G,INVALID);
114      VecNode level_list(n,INVALID);
115      //List of the nodes in level i<n, set to n.
116
117      NodeIt v;
118      for(G.first(v); G.valid(v); G.next(v)) level.set(v,n);
119      //setting each node to level n
120     
121      switch ( fe ) {
122      case PREFLOW:
123        {
124          //counting the excess
125          NodeIt v;
126          for(G.first(v); G.valid(v); G.next(v)) {
127            T exc=0;
128         
129            InEdgeIt e;
130            for(G.first(e,v); G.valid(e); G.next(e)) exc+=(*flow)[e];
131            OutEdgeIt f;
132            for(G.first(f,v); G.valid(f); G.next(f)) exc-=(*flow)[f];
133           
134            excess.set(v,exc);   
135           
136            //putting the active nodes into the stack
137            int lev=level[v];
138            if ( exc > 0 && lev < n && v != t ) active[lev].push(v);
139          }
140          break;
141        }
142      case GEN_FLOW:
143        {
144          //Counting the excess of t
145          T exc=0;
146         
147          InEdgeIt e;
148          for(G.first(e,t); G.valid(e); G.next(e)) exc+=(*flow)[e];
149          OutEdgeIt f;
150          for(G.first(f,t); G.valid(f); G.next(f)) exc-=(*flow)[f];
151         
152          excess.set(t,exc);   
153         
154          break;
155        }
156      }
157     
158      preflowPreproc( fe, active, level_list, left, right );
159      //End of preprocessing
160     
161     
162      //Push/relabel on the highest level active nodes.
163      while ( true ) {
164        if ( b == 0 ) {
165          if ( !what_heur && !end && k > 0 ) {
166            b=k;
167            end=true;
168          } else break;
169        }
170       
171        if ( active[b].empty() ) --b;
172        else {
173          end=false; 
174          Node w=active[b].top();
175          active[b].pop();
176          int newlevel=push(w,active);
177          if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list,
178                                       left, right, b, k, what_heur);
179         
180          ++numrelabel;
181          if ( numrelabel >= heur ) {
182            numrelabel=0;
183            if ( what_heur ) {
184              what_heur=0;
185              heur=heur0;
186              end=false;
187            } else {
188              what_heur=1;
189              heur=heur1;
190              b=k;
191            }
192          }
193        }
194      }
195    }
196
197
198    void preflowPhase1() {
199     
200      int k=n-2;  //bound on the highest level under n containing a node
201      int b=k;    //bound on the highest level under n of an active node
202     
203      VecStack active(n);
204      level.set(s,0);
205      std::queue<Node> bfs_queue;
206      bfs_queue.push(s);
207           
208      while (!bfs_queue.empty()) {
209       
210        Node v=bfs_queue.front();       
211        bfs_queue.pop();
212        int l=level[v]+1;
213             
214        InEdgeIt e;
215        for(G.first(e,v); G.valid(e); G.next(e)) {
216          if ( (*capacity)[e] == (*flow)[e] ) continue;
217          Node u=G.tail(e);
218          if ( level[u] >= n ) {
219            bfs_queue.push(u);
220            level.set(u, l);
221            if ( excess[u] > 0 ) active[l].push(u);
222          }
223        }
224       
225        OutEdgeIt f;
226        for(G.first(f,v); G.valid(f); G.next(f)) {
227          if ( 0 == (*flow)[f] ) continue;
228          Node u=G.head(f);
229          if ( level[u] >= n ) {
230            bfs_queue.push(u);
231            level.set(u, l);
232            if ( excess[u] > 0 ) active[l].push(u);
233          }
234        }
235      }
236      b=n-2;
237
238      while ( true ) {
239       
240        if ( b == 0 ) break;
241
242        if ( active[b].empty() ) --b;
243        else {
244          Node w=active[b].top();
245          active[b].pop();
246          int newlevel=push(w,active);   
247
248          //relabel
249          if ( excess[w] > 0 ) {
250            level.set(w,++newlevel);
251            active[newlevel].push(w);
252            b=newlevel;
253          }
254        }  // if stack[b] is nonempty
255      } // while(true)
256    }
257
258
259    //Returns the maximum value of a flow.
260    T flowValue() {
261      return excess[t];
262    }
263
264    //should be used only between preflowPhase0 and preflowPhase1
265    template<typename _CutMap>
266    void actMinCut(_CutMap& M) {
267      NodeIt v;
268      for(G.first(v); G.valid(v); G.next(v))
269        if ( level[v] < n ) M.set(v,false);
270        else M.set(v,true);
271    }
272
273
274
275    /*
276      Returns the minimum min cut, by a bfs from s in the residual graph.
277    */
278    template<typename _CutMap>
279    void minMinCut(_CutMap& M) {
280   
281      std::queue<Node> queue;
282     
283      M.set(s,true);     
284      queue.push(s);
285
286      while (!queue.empty()) {
287        Node w=queue.front();
288        queue.pop();
289
290        OutEdgeIt e;
291        for(G.first(e,w) ; G.valid(e); G.next(e)) {
292          Node v=G.head(e);
293          if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
294            queue.push(v);
295            M.set(v, true);
296          }
297        }
298
299        InEdgeIt f;
300        for(G.first(f,w) ; G.valid(f); G.next(f)) {
301          Node v=G.tail(f);
302          if (!M[v] && (*flow)[f] > 0 ) {
303            queue.push(v);
304            M.set(v, true);
305          }
306        }
307      }
308    }
309
310
311 
312    /*
313      Returns the maximum min cut, by a reverse bfs
314      from t in the residual graph.
315    */
316   
317    template<typename _CutMap>
318    void maxMinCut(_CutMap& M) {
319
320      NodeIt v;
321      for(G.first(v) ; G.valid(v); G.next(v)) {
322        M.set(v, true);
323      }
324
325      std::queue<Node> queue;
326     
327      M.set(t,false);       
328      queue.push(t);
329
330      while (!queue.empty()) {
331        Node w=queue.front();
332        queue.pop();
333
334
335        InEdgeIt e;
336        for(G.first(e,w) ; G.valid(e); G.next(e)) {
337          Node v=G.tail(e);
338          if (M[v] && (*flow)[e] < (*capacity)[e] ) {
339            queue.push(v);
340            M.set(v, false);
341          }
342        }
343       
344        OutEdgeIt f;
345        for(G.first(f,w) ; G.valid(f); G.next(f)) {
346          Node v=G.head(f);
347          if (M[v] && (*flow)[f] > 0 ) {
348            queue.push(v);
349            M.set(v, false);
350          }
351        }
352      }
353    }
354
355
356    template<typename CutMap>
357    void minCut(CutMap& M) {
358      minMinCut(M);
359    }
360
361   
362    void resetTarget (const Node _t) {t=_t;}
363
364    void resetSource (const Node _s) {s=_s;}
365   
366    void resetCap (const CapMap& _cap) {
367      capacity=&_cap;
368    }
369   
370    void resetFlow (FlowMap& _flow) {
371      flow=&_flow;
372    }
373
374
375  private:
376
377    int push(const Node w, VecStack& active) {
378     
379      int lev=level[w];
380      T exc=excess[w];
381      int newlevel=n;       //bound on the next level of w
382         
383      OutEdgeIt e;
384      for(G.first(e,w); G.valid(e); G.next(e)) {
385           
386        if ( (*flow)[e] == (*capacity)[e] ) continue;
387        Node v=G.head(e);           
388           
389        if( lev > level[v] ) { //Push is allowed now
390         
391          if ( excess[v]==0 && v!=t && v!=s ) {
392            int lev_v=level[v];
393            active[lev_v].push(v);
394          }
395         
396          T cap=(*capacity)[e];
397          T flo=(*flow)[e];
398          T remcap=cap-flo;
399         
400          if ( remcap >= exc ) { //A nonsaturating push.
401           
402            flow->set(e, flo+exc);
403            excess.set(v, excess[v]+exc);
404            exc=0;
405            break;
406           
407          } else { //A saturating push.
408            flow->set(e, cap);
409            excess.set(v, excess[v]+remcap);
410            exc-=remcap;
411          }
412        } else if ( newlevel > level[v] ) newlevel = level[v];
413      } //for out edges wv
414     
415      if ( exc > 0 ) { 
416        InEdgeIt e;
417        for(G.first(e,w); G.valid(e); G.next(e)) {
418         
419          if( (*flow)[e] == 0 ) continue;
420          Node v=G.tail(e);
421         
422          if( lev > level[v] ) { //Push is allowed now
423           
424            if ( excess[v]==0 && v!=t && v!=s ) {
425              int lev_v=level[v];
426              active[lev_v].push(v);
427            }
428           
429            T flo=(*flow)[e];
430           
431            if ( flo >= exc ) { //A nonsaturating push.
432             
433              flow->set(e, flo-exc);
434              excess.set(v, excess[v]+exc);
435              exc=0;
436              break;
437            } else {  //A saturating push.
438             
439              excess.set(v, excess[v]+flo);
440              exc-=flo;
441              flow->set(e,0);
442            } 
443          } else if ( newlevel > level[v] ) newlevel = level[v];
444        } //for in edges vw
445       
446      } // if w still has excess after the out edge for cycle
447     
448      excess.set(w, exc);
449     
450      return newlevel;
451     }
452
453
454    void preflowPreproc ( flowEnum fe, VecStack& active,
455                          VecNode& level_list, NNMap& left, NNMap& right ) {
456
457      std::queue<Node> bfs_queue;
458     
459      switch ( fe ) {
460      case ZERO_FLOW:
461        {
462          //Reverse_bfs from t, to find the starting level.
463          level.set(t,0);
464          bfs_queue.push(t);
465       
466          while (!bfs_queue.empty()) {
467           
468            Node v=bfs_queue.front();   
469            bfs_queue.pop();
470            int l=level[v]+1;
471           
472            InEdgeIt e;
473            for(G.first(e,v); G.valid(e); G.next(e)) {
474              Node w=G.tail(e);
475              if ( level[w] == n && w != s ) {
476                bfs_queue.push(w);
477                Node first=level_list[l];
478                if ( G.valid(first) ) left.set(first,w);
479                right.set(w,first);
480                level_list[l]=w;
481                level.set(w, l);
482              }
483            }
484          }
485         
486          //the starting flow
487          OutEdgeIt e;
488          for(G.first(e,s); G.valid(e); G.next(e))
489            {
490              T c=(*capacity)[e];
491              if ( c == 0 ) continue;
492              Node w=G.head(e);
493              if ( level[w] < n ) {       
494                if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);
495                flow->set(e, c);
496                excess.set(w, excess[w]+c);
497              }
498            }
499          break;
500        }
501       
502      case GEN_FLOW:
503      case PREFLOW:
504        {
505          //Reverse_bfs from t in the residual graph,
506          //to find the starting level.
507          level.set(t,0);
508          bfs_queue.push(t);
509         
510          while (!bfs_queue.empty()) {
511           
512            Node v=bfs_queue.front();   
513            bfs_queue.pop();
514            int l=level[v]+1;
515           
516            InEdgeIt e;
517            for(G.first(e,v); G.valid(e); G.next(e)) {
518              if ( (*capacity)[e] == (*flow)[e] ) continue;
519              Node w=G.tail(e);
520              if ( level[w] == n && w != s ) {
521                bfs_queue.push(w);
522                Node first=level_list[l];
523                if ( G.valid(first) ) left.set(first,w);
524                right.set(w,first);
525                level_list[l]=w;
526                level.set(w, l);
527              }
528            }
529           
530            OutEdgeIt f;
531            for(G.first(f,v); G.valid(f); G.next(f)) {
532              if ( 0 == (*flow)[f] ) continue;
533              Node w=G.head(f);
534              if ( level[w] == n && w != s ) {
535                bfs_queue.push(w);
536                Node first=level_list[l];
537                if ( G.valid(first) ) left.set(first,w);
538                right.set(w,first);
539                level_list[l]=w;
540                level.set(w, l);
541              }
542            }
543          }
544         
545         
546          //the starting flow
547          OutEdgeIt e;
548          for(G.first(e,s); G.valid(e); G.next(e))
549            {
550              T rem=(*capacity)[e]-(*flow)[e];
551              if ( rem == 0 ) continue;
552              Node w=G.head(e);
553              if ( level[w] < n ) {       
554                if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);
555                flow->set(e, (*capacity)[e]);
556                excess.set(w, excess[w]+rem);
557              }
558            }
559         
560          InEdgeIt f;
561          for(G.first(f,s); G.valid(f); G.next(f))
562            {
563              if ( (*flow)[f] == 0 ) continue;
564              Node w=G.tail(f);
565              if ( level[w] < n ) {       
566                if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);
567                excess.set(w, excess[w]+(*flow)[f]);
568                flow->set(f, 0);
569              }
570            } 
571          break;
572        } //case PREFLOW
573      }
574    } //preflowPreproc
575
576
577
578    void relabel( const Node w, int newlevel, VecStack& active, 
579                  VecNode& level_list, NNMap& left,
580                  NNMap& right, int& b, int& k, const bool what_heur ) {
581
582      T lev=level[w];   
583     
584      Node right_n=right[w];
585      Node left_n=left[w];
586     
587      //unlacing starts
588      if ( G.valid(right_n) ) {
589        if ( G.valid(left_n) ) {
590          right.set(left_n, right_n);
591          left.set(right_n, left_n);
592        } else {
593          level_list[lev]=right_n;   
594          left.set(right_n, INVALID);
595        }
596      } else {
597        if ( G.valid(left_n) ) {
598          right.set(left_n, INVALID);
599        } else {
600          level_list[lev]=INVALID;   
601        }
602      }
603      //unlacing ends
604               
605      if ( !G.valid(level_list[lev]) ) {
606             
607        //gapping starts
608        for (int i=lev; i!=k ; ) {
609          Node v=level_list[++i];
610          while ( G.valid(v) ) {
611            level.set(v,n);
612            v=right[v];
613          }
614          level_list[i]=INVALID;
615          if ( !what_heur ) {
616            while ( !active[i].empty() ) {
617              active[i].pop();    //FIXME: ezt szebben kene
618            }
619          }         
620        }
621       
622        level.set(w,n);
623        b=lev-1;
624        k=b;
625        //gapping ends
626       
627      } else {
628       
629        if ( newlevel == n ) level.set(w,n);
630        else {
631          level.set(w,++newlevel);
632          active[newlevel].push(w);
633          if ( what_heur ) b=newlevel;
634          if ( k < newlevel ) ++k;      //now k=newlevel
635          Node first=level_list[newlevel];
636          if ( G.valid(first) ) left.set(first,w);
637          right.set(w,first);
638          left.set(w,INVALID);
639          level_list[newlevel]=w;
640        }
641      }
642     
643    } //relabel
644   
645
646  };
647
648} //namespace hugo
649
650#endif //HUGO_PREFLOW_H
651
652
653
654
Note: See TracBrowser for help on using the repository browser.