COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/preflow.h @ 278:c11f84e3da21

Last change on this file since 278:c11f84e3da21 was 278:c11f84e3da21, checked in by marci, 16 years ago

const Graph&, const CapMap?& in preflow constructor

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