COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/preflow_res.h @ 478:8c74de352f80

Last change on this file since 478:8c74de352f80 was 444:618c5d6f36b9, checked in by jacint, 20 years ago

debug

File size: 9.7 KB
Line 
1// -*- C++ -*-
2//The same as preflow.h, using ResGraphWrapper
3#ifndef HUGO_PREFLOW_RES_H
4#define HUGO_PREFLOW_RES_H
5
6#define H0 20
7#define H1 1
8
9#include <vector>
10#include <queue>
11#include <graph_wrapper.h>
12
13#include<iostream>
14
15namespace hugo {
16
17  template <typename Graph, typename T,
18            typename CapMap=typename Graph::template EdgeMap<T>,
19            typename FlowMap=typename Graph::template EdgeMap<T> >
20  class PreflowRes {
21   
22    typedef typename Graph::Node Node;
23    typedef typename Graph::Edge Edge;
24    typedef typename Graph::NodeIt NodeIt;
25    typedef typename Graph::OutEdgeIt OutEdgeIt;
26    typedef typename Graph::InEdgeIt InEdgeIt;
27   
28    const Graph& G;
29    Node s;
30    Node t;
31    const CapMap& capacity; 
32    FlowMap& flow;
33    T value;
34    bool constzero;
35
36    typedef ResGraphWrapper<const Graph, T, CapMap, FlowMap> ResGW;
37    typedef typename ResGW::OutEdgeIt ResOutEdgeIt;
38    typedef typename ResGW::InEdgeIt ResInEdgeIt;
39    typedef typename ResGW::Edge ResEdge;
40 
41  public:
42    PreflowRes(Graph& _G, Node _s, Node _t, CapMap& _capacity,
43            FlowMap& _flow, bool _constzero ) :
44      G(_G), s(_s), t(_t), capacity(_capacity), flow(_flow), constzero(_constzero) {}
45   
46   
47    void run() {
48
49      ResGW res_graph(G, capacity, flow);
50
51      value=0;                //for the subsequent runs
52
53      bool phase=0;        //phase 0 is the 1st phase, phase 1 is the 2nd
54      int n=G.nodeNum();
55      int heur0=(int)(H0*n);  //time while running 'bound decrease'
56      int heur1=(int)(H1*n);  //time while running 'highest label'
57      int heur=heur1;         //starting time interval (#of relabels)
58      bool what_heur=1;       
59      /*
60        what_heur is 0 in case 'bound decrease'
61        and 1 in case 'highest label'
62      */
63      bool end=false;     
64      /*
65        Needed for 'bound decrease', 'true'
66        means no active nodes are above bound b.
67      */
68      int relabel=0;
69      int k=n-2;  //bound on the highest level under n containing a node
70      int b=k;    //bound on the highest level under n of an active node
71     
72      typename Graph::template NodeMap<int> level(G,n);     
73      typename Graph::template NodeMap<T> excess(G);
74
75      std::vector<Node> active(n-1,INVALID);
76      typename Graph::template NodeMap<Node> next(G,INVALID);
77      //Stack of the active nodes in level i < n.
78      //We use it in both phases.
79
80      typename Graph::template NodeMap<Node> left(G,INVALID);
81      typename Graph::template NodeMap<Node> right(G,INVALID);
82      std::vector<Node> level_list(n,INVALID);
83      /*
84        List of the nodes in level i<n.
85      */
86
87
88      /*
89        Reverse_bfs from t in the residual graph,
90        to find the starting level.
91      */
92      level.set(t,0);
93      std::queue<Node> bfs_queue;
94      bfs_queue.push(t);
95     
96      while (!bfs_queue.empty()) {
97       
98        Node v=bfs_queue.front();       
99        bfs_queue.pop();
100        int l=level[v]+1;
101       
102        ResInEdgeIt e;
103        for(res_graph.first(e,v); res_graph.valid(e);
104            res_graph.next(e)) {
105          Node w=res_graph.tail(e);
106          if ( level[w] == n && w != s ) {
107            bfs_queue.push(w);
108            Node first=level_list[l];
109            if ( G.valid(first) ) left.set(first,w);
110            right.set(w,first);
111            level_list[l]=w;
112            level.set(w, l);
113          }
114        }
115      }
116     
117       
118      if ( !constzero ) {
119        /*
120          Counting the excess
121        */
122        NodeIt v;
123        for(G.first(v); G.valid(v); G.next(v)) {
124          T exc=0;
125
126          InEdgeIt e;
127          for(G.first(e,v); G.valid(e); G.next(e)) exc+=flow[e];
128          OutEdgeIt f;
129          for(G.first(f,v); G.valid(f); G.next(f)) exc-=flow[f];
130
131          excess.set(v,exc);     
132
133          //putting the active nodes into the stack
134          int lev=level[v];
135          if ( exc > 0 && lev < n ) {
136            next.set(v,active[lev]);
137            active[lev]=v;
138          }
139        }
140      }
141     
142
143
144      //the starting flow
145      ResOutEdgeIt e;
146      for(res_graph.first(e,s); res_graph.valid(e);
147          res_graph.next(e)) {
148          Node w=res_graph.head(e);
149          if ( level[w] < n ) {   
150            if ( excess[w] == 0 && w!=t ) {
151              next.set(w,active[level[w]]);
152              active[level[w]]=w;
153            }
154            T rem=res_graph.resCap(e);
155            excess.set(w, excess[w]+rem);
156            res_graph.augment(e, rem );
157          }
158      }
159       
160
161      /*
162         End of preprocessing
163      */
164
165
166
167      /*
168        Push/relabel on the highest level active nodes.
169      */       
170      while ( true ) {
171       
172        if ( b == 0 ) {
173          if ( phase ) break;
174         
175          if ( !what_heur && !end && k > 0 ) {
176            b=k;
177            end=true;
178          } else {
179            phase=1;
180            level.set(s,0);
181            std::queue<Node> bfs_queue;
182            bfs_queue.push(s);
183           
184            while (!bfs_queue.empty()) {
185             
186              Node v=bfs_queue.front();
187              bfs_queue.pop();
188              int l=level[v]+1;
189             
190              ResInEdgeIt e;
191              for(res_graph.first(e,v);
192                  res_graph.valid(e); res_graph.next(e)) {
193                Node u=res_graph.tail(e);
194                if ( level[u] >= n ) {
195                  bfs_queue.push(u);
196                  level.set(u, l);
197                  if ( excess[u] > 0 ) {
198                    next.set(u,active[l]);
199                    active[l]=u;
200                  }
201                }
202              }
203           
204            }
205            b=n-2;
206          }
207           
208        }
209         
210         
211        if ( !G.valid(active[b]) ) --b;
212        else {
213          end=false; 
214
215          Node w=active[b];
216          active[b]=next[w];
217          int lev=level[w];
218          T exc=excess[w];
219          int newlevel=n;       //bound on the next level of w
220         
221          ResOutEdgeIt e;
222          for(res_graph.first(e,w); res_graph.valid(e); res_graph.next(e)) {
223           
224            Node v=res_graph.head(e);           
225            if( lev > level[v] ) {     
226              /*Push is allowed now*/
227             
228              if ( excess[v]==0 && v!=t && v!=s ) {
229                int lev_v=level[v];
230                next.set(v,active[lev_v]);
231                active[lev_v]=v;
232              }
233             
234              T remcap=res_graph.resCap(e);
235             
236              if ( remcap >= exc ) {       
237                /*A nonsaturating push.*/
238                res_graph.augment(e, exc);
239                excess.set(v, excess[v]+exc);
240                exc=0;
241                break;
242               
243              } else {
244                /*A saturating push.*/
245               
246                res_graph.augment(e, remcap);
247                excess.set(v, excess[v]+remcap);
248                exc-=remcap;
249              }
250            } else if ( newlevel > level[v] ){
251              newlevel = level[v];
252            }       
253           
254          }
255
256        excess.set(w, exc);
257         
258        /*
259          Relabel
260        */
261       
262
263        if ( exc > 0 ) {
264          //now 'lev' is the old level of w
265       
266          if ( phase ) {
267            level.set(w,++newlevel);
268            next.set(w,active[newlevel]);
269            active[newlevel]=w;
270            b=newlevel;
271          } else {
272            //unlacing starts
273            Node right_n=right[w];
274            Node left_n=left[w];
275
276            if ( G.valid(right_n) ) {
277              if ( G.valid(left_n) ) {
278                right.set(left_n, right_n);
279                left.set(right_n, left_n);
280              } else {
281                level_list[lev]=right_n;   
282                left.set(right_n, INVALID);
283              }
284            } else {
285              if ( G.valid(left_n) ) {
286                right.set(left_n, INVALID);
287              } else {
288                level_list[lev]=INVALID;   
289              }
290            }
291            //unlacing ends
292               
293            if ( !G.valid(level_list[lev]) ) {
294             
295               //gapping starts
296              for (int i=lev; i!=k ; ) {
297                Node v=level_list[++i];
298                while ( G.valid(v) ) {
299                  level.set(v,n);
300                  v=right[v];
301                }
302                level_list[i]=INVALID;
303                if ( !what_heur ) active[i]=INVALID;
304              }     
305
306              level.set(w,n);
307              b=lev-1;
308              k=b;
309              //gapping ends
310           
311            } else {
312             
313              if ( newlevel == n ) level.set(w,n);
314              else {
315                level.set(w,++newlevel);
316                next.set(w,active[newlevel]);
317                active[newlevel]=w;
318                if ( what_heur ) b=newlevel;
319                if ( k < newlevel ) ++k;      //now k=newlevel
320                Node first=level_list[newlevel];
321                if ( G.valid(first) ) left.set(first,w);
322                right.set(w,first);
323                left.set(w,INVALID);
324                level_list[newlevel]=w;
325              }
326            }
327
328
329            ++relabel;
330            if ( relabel >= heur ) {
331              relabel=0;
332              if ( what_heur ) {
333                what_heur=0;
334                heur=heur0;
335                end=false;
336              } else {
337                what_heur=1;
338                heur=heur1;
339                b=k;
340              }
341            }
342          } //phase 0
343         
344         
345        } // if ( exc > 0 )
346         
347       
348        }  // if stack[b] is nonempty
349       
350      } // while(true)
351
352
353      value = excess[t];
354      /*Max flow value.*/
355     
356    } //void run()
357
358
359
360
361
362    /*
363      Returns the maximum value of a flow.
364     */
365
366    T flowValue() {
367      return value;
368    }
369
370
371    FlowMap Flow() {
372      return flow;
373      }
374
375
376   
377    void Flow(FlowMap& _flow ) {
378      NodeIt v;
379      for(G.first(v) ; G.valid(v); G.next(v))
380        _flow.set(v,flow[v]);
381    }
382
383
384
385    /*
386      Returns the minimum min cut, by a bfs from s in the residual graph.
387    */
388   
389    template<typename _CutMap>
390    void minMinCut(_CutMap& M) {
391   
392      std::queue<Node> queue;
393     
394      M.set(s,true);     
395      queue.push(s);
396
397      while (!queue.empty()) {
398        Node w=queue.front();
399        queue.pop();
400
401        OutEdgeIt e;
402        for(G.first(e,w) ; G.valid(e); G.next(e)) {
403          Node v=G.head(e);
404          if (!M[v] && flow[e] < capacity[e] ) {
405            queue.push(v);
406            M.set(v, true);
407          }
408        }
409
410        InEdgeIt f;
411        for(G.first(f,w) ; G.valid(f); G.next(f)) {
412          Node v=G.tail(f);
413          if (!M[v] && flow[f] > 0 ) {
414            queue.push(v);
415            M.set(v, true);
416          }
417        }
418      }
419    }
420
421
422 
423    /*
424      Returns the maximum min cut, by a reverse bfs
425      from t in the residual graph.
426    */
427   
428    template<typename _CutMap>
429    void maxMinCut(_CutMap& M) {
430   
431      std::queue<Node> queue;
432     
433      M.set(t,true);       
434      queue.push(t);
435
436      while (!queue.empty()) {
437        Node w=queue.front();
438        queue.pop();
439
440
441        InEdgeIt e;
442        for(G.first(e,w) ; G.valid(e); G.next(e)) {
443          Node v=G.tail(e);
444          if (!M[v] && flow[e] < capacity[e] ) {
445            queue.push(v);
446            M.set(v, true);
447          }
448        }
449       
450        OutEdgeIt f;
451        for(G.first(f,w) ; G.valid(f); G.next(f)) {
452          Node v=G.head(f);
453          if (!M[v] && flow[f] > 0 ) {
454            queue.push(v);
455            M.set(v, true);
456          }
457        }
458      }
459
460      NodeIt v;
461      for(G.first(v) ; G.valid(v); G.next(v)) {
462        M.set(v, !M[v]);
463      }
464
465    }
466
467
468
469    template<typename CutMap>
470    void minCut(CutMap& M) {
471      minMinCut(M);
472    }
473
474   
475   
476    void resetTarget (Node _t) {t=_t;}
477    void resetSource (Node _s) {s=_s;}
478   
479    void resetCap (CapMap _cap) {capacity=_cap;}
480
481    void resetFlow (FlowMap _flow, bool _constzero) {
482      flow=_flow;
483      constzero=_constzero;
484    }
485
486
487  };
488
489} //namespace hugo
490
491#endif //HUGO_PREFLOW_RES_H
492
493
494
495
Note: See TracBrowser for help on using the repository browser.