COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/preflow.h @ 350:3a9a767b841e

Last change on this file since 350:3a9a767b841e was 330:7ac0d4e8a31c, checked in by marci, 20 years ago

In the resgraphwrapper interface, and in the constructor,
the order of FlowMap? and CapacityMap? is changed.

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