COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/preflow_jgraph.h @ 140:ca164520d31a

Last change on this file since 140:ca164520d31a was 140:ca164520d31a, checked in by jacint, 16 years ago

* empty log message *

File size: 10.8 KB
Line 
1// -*- C++ -*-
2/*
3preflow.h with 'j_graph interface'
4by jacint.
5Heuristics:
6 2 phase
7 gap
8 list 'level_list' on the nodes on level i implemented by hand
9 stack 'active' on the active nodes on level i implemented by hand
10 runs heuristic 'highest label' for H1*n relabels
11 runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
12 
13Parameters H0 and H1 are initialized to 20 and 10.
14
15The best preflow I could ever write.
16
17The constructor runs the algorithm.
18
19Members:
20
21T maxFlow() : returns the value of a maximum flow
22
23T flowOnEdge(EdgeIt e) : for a fixed maximum flow x it returns x(e)
24
25FlowMap Flow() : returns the fixed maximum flow x
26
27void minMinCut(CutMap& M) : sets M to the characteristic vector of the
28     minimum min cut. M should be a map of bools initialized to false.
29
30void maxMinCut(CutMap& M) : sets M to the characteristic vector of the
31     maximum min cut. M should be a map of bools initialized to false.
32
33void minCut(CutMap& M) : sets M to the characteristic vector of
34     a min cut. M should be a map of bools initialized to false.
35
36*/
37
38#ifndef PREFLOW_H
39#define PREFLOW_H
40
41#define H0 20
42#define H1 1
43
44#include <vector>
45#include <queue>
46
47#include<iostream>
48
49#include <time_measure.h>
50
51namespace hugo {
52
53  template <typename Graph, typename T,
54    typename FlowMap=typename Graph::EdgeMap<T>,
55    typename CapMap=typename Graph::EdgeMap<T> >
56  class preflow {
57   
58    typedef typename Graph::TrivNodeIt NodeIt;
59    typedef typename Graph::TrivEdgeIt EdgeIt;
60    typedef typename Graph::NodeIt EachNodeIt;
61    typedef typename Graph::OutEdgeIt OutEdgeIt;
62    typedef typename Graph::InEdgeIt InEdgeIt;
63   
64    Graph& G;
65    NodeIt s;
66    NodeIt t;
67    FlowMap flow;
68    CapMap& capacity; 
69    T value;
70
71  public:
72    double time;
73    preflow(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity ) :
74      G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity)
75    {
76
77      bool phase=0;        //phase 0 is the 1st phase, phase 1 is the 2nd
78      int n=G.numNodes();
79      int heur0=(int)(H0*n);  //time while running 'bound decrease'
80      int heur1=(int)(H1*n);  //time while running 'highest label'
81      int heur=heur1;         //starting time interval (#of relabels)
82      bool what_heur=1;       
83      /*
84        what_heur is 0 in case 'bound decrease'
85        and 1 in case 'highest label'
86      */
87      bool end=false;     
88      /*
89        Needed for 'bound decrease', 'true'
90        means no active nodes are above bound b.
91      */
92      int relabel=0;
93      int k=n-2;  //bound on the highest level under n containing a node
94      int b=k;    //bound on the highest level under n of an active node
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        List of the nodes in level i<n.
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.firstIn(v); e; G.next(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 ) 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.firstOut(s); e; G.next(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            phase=1;
173            time=currTime();
174            level.set(s,0);
175            std::queue<NodeIt> bfs_queue;
176            bfs_queue.push(s);
177           
178            while (!bfs_queue.empty()) {
179             
180              NodeIt v=bfs_queue.front();       
181              bfs_queue.pop();
182              int l=level.get(v)+1;
183             
184              for(InEdgeIt e=G.firstIn(v); e; G.next(e)) {
185                if ( capacity.get(e) == flow.get(e) ) continue;
186                NodeIt u=G.tail(e);
187                if ( level.get(u) >= n ) {
188                  bfs_queue.push(u);
189                  level.set(u, l);
190                  if ( excess.get(u) > 0 ) {
191                    next.set(u,active[l]);
192                    active[l]=u;
193                  }
194                }
195              }
196           
197              for(OutEdgeIt e=G.firstOut(v); e; G.next(e)) {
198                if ( 0 == flow.get(e) ) continue;
199                NodeIt u=G.head(e);
200                if ( level.get(u) >= n ) {
201                  bfs_queue.push(u);
202                  level.set(u, l);
203                  if ( excess.get(u) > 0 ) {
204                    next.set(u,active[l]);
205                    active[l]=u;
206                  }
207                }
208              }
209            }
210            b=n-2;
211            }
212           
213        }
214         
215         
216        if ( !active[b]  ) --b;
217        else {
218          end=false; 
219
220          NodeIt w=active[b];
221          active[b]=next.get(w);
222          int lev=level.get(w);
223          T exc=excess.get(w);
224          int newlevel=n;       //bound on the next level of w
225         
226          for(OutEdgeIt e=G.firstOut(w); e; G.next(e)) {
227           
228            if ( flow.get(e) == capacity.get(e) ) continue;
229            NodeIt v=G.head(e);           
230            //e=wv         
231           
232            if( lev > level.get(v) ) {     
233              /*Push is allowed now*/
234             
235              if ( excess.get(v)==0 && v!=t && v!=s ) {
236                int lev_v=level.get(v);
237                next.set(v,active[lev_v]);
238                active[lev_v]=v;
239              }
240             
241              T cap=capacity.get(e);
242              T flo=flow.get(e);
243              T remcap=cap-flo;
244             
245              if ( remcap >= exc ) {       
246                /*A nonsaturating push.*/
247               
248                flow.set(e, flo+exc);
249                excess.set(v, excess.get(v)+exc);
250                exc=0;
251                break;
252               
253              } else {
254                /*A saturating push.*/
255               
256                flow.set(e, cap);
257                excess.set(v, excess.get(v)+remcap);
258                exc-=remcap;
259              }
260            } else if ( newlevel > level.get(v) ){
261              newlevel = level.get(v);
262            }       
263           
264          } //for out edges wv
265       
266       
267        if ( exc > 0 ) {       
268          for( InEdgeIt e=G.firstIn(w); e; G.next(e)) {
269           
270            if( flow.get(e) == 0 ) continue;
271            NodeIt v=G.tail(e); 
272            //e=vw
273           
274            if( lev > level.get(v) ) { 
275              /*Push is allowed now*/
276             
277              if ( excess.get(v)==0 && v!=t && v!=s ) {
278                int lev_v=level.get(v);
279                next.set(v,active[lev_v]);
280                active[lev_v]=v;
281              }
282             
283              T flo=flow.get(e);
284             
285              if ( flo >= exc ) {
286                /*A nonsaturating push.*/
287               
288                flow.set(e, flo-exc);
289                excess.set(v, excess.get(v)+exc);
290                exc=0;
291                break;
292              } else {                                               
293/*A saturating push.*/
294               
295                excess.set(v, excess.get(v)+flo);
296                exc-=flo;
297                flow.set(e,0);
298              } 
299            } else if ( newlevel > level.get(v) ) {
300              newlevel = level.get(v);
301            }       
302          } //for in edges vw
303         
304        } // if w still has excess after the out edge for cycle
305       
306        excess.set(w, exc);
307         
308        /*
309          Relabel
310        */
311       
312
313        if ( exc > 0 ) {
314          //now 'lev' is the old level of w
315       
316          if ( phase ) {
317            level.set(w,++newlevel);
318            next.set(w,active[newlevel]);
319            active[newlevel]=w;
320            b=newlevel;
321          } else {
322            //unlacing starts
323            NodeIt right_n=right.get(w);
324            NodeIt left_n=left.get(w);
325
326            if ( right_n ) {
327              if ( left_n ) {
328                right.set(left_n, right_n);
329                left.set(right_n, left_n);
330              } else {
331                level_list[lev]=right_n;   
332                left.set(right_n, NodeIt());
333              }
334            } else {
335              if ( left_n ) {
336                right.set(left_n, NodeIt());
337              } else {
338                level_list[lev]=NodeIt();   
339
340              }
341            }
342            //unlacing ends
343               
344            //gapping starts
345            if ( !level_list[lev] ) {
346             
347              for (int i=lev; i!=k ; ) {
348                NodeIt v=level_list[++i];
349                while ( v ) {
350                  level.set(v,n);
351                  v=right.get(v);
352                }
353                level_list[i]=NodeIt();
354                if ( !what_heur ) active[i]=NodeIt();
355              }     
356
357              level.set(w,n);
358              b=lev-1;
359              k=b;
360              //gapping ends
361            } else {
362             
363              if ( newlevel == n ) level.set(w,n);
364              else {
365                level.set(w,++newlevel);
366                next.set(w,active[newlevel]);
367                active[newlevel]=w;
368                if ( what_heur ) b=newlevel;
369                if ( k < newlevel ) ++k;
370                NodeIt first=level_list[newlevel];
371                if ( first ) left.set(first,w);
372                right.set(w,first);
373                left.set(w,NodeIt());
374                level_list[newlevel]=w;
375              }
376            }
377
378
379            ++relabel;
380            if ( relabel >= heur ) {
381              relabel=0;
382              if ( what_heur ) {
383                what_heur=0;
384                heur=heur0;
385                end=false;
386              } else {
387                what_heur=1;
388                heur=heur1;
389                b=k;
390              }
391            }
392          } //phase 0
393         
394         
395        } // if ( exc > 0 )
396         
397       
398        }  // if stack[b] is nonempty
399       
400      } // while(true)
401
402
403      value = excess.get(t);
404      /*Max flow value.*/
405     
406    } //void run()
407
408
409
410
411
412    /*
413      Returns the maximum value of a flow.
414     */
415
416    T maxFlow() {
417      return value;
418    }
419
420
421
422    /*
423      For the maximum flow x found by the algorithm,
424      it returns the flow value on edge e, i.e. x(e).
425    */
426   
427    T flowOnEdge(EdgeIt e) {
428      return flow.get(e);
429    }
430
431
432
433    FlowMap Flow() {
434      return flow;
435      }
436
437
438   
439    void Flow(FlowMap& _flow ) {
440      for(EachNodeIt v=G.firstNode() ; v; G.next(v))
441        _flow.set(v,flow.get(v));
442      }
443   
444
445
446    /*
447      Returns the minimum min cut, by a bfs from s in the residual graph.
448    */
449   
450    template<typename _CutMap>
451    void minMinCut(_CutMap& M) {
452   
453      std::queue<NodeIt> queue;
454     
455      M.set(s,true);     
456      queue.push(s);
457
458      while (!queue.empty()) {
459        NodeIt w=queue.front();
460        queue.pop();
461
462        for(OutEdgeIt e=G.firstOut(w) ; e; G.next(e)) {
463          NodeIt v=G.head(e);
464          if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
465            queue.push(v);
466            M.set(v, true);
467          }
468        }
469
470        for(InEdgeIt e=G.firstIn(w) ; e; G.next(e)) {
471          NodeIt v=G.tail(e);
472          if (!M.get(v) && flow.get(e) > 0 ) {
473            queue.push(v);
474            M.set(v, true);
475          }
476        }
477      }
478    }
479
480
481 
482    /*
483      Returns the maximum min cut, by a reverse bfs
484      from t in the residual graph.
485    */
486   
487    template<typename _CutMap>
488    void maxMinCut(_CutMap& M) {
489   
490      std::queue<NodeIt> queue;
491     
492      M.set(t,true);       
493      queue.push(t);
494
495      while (!queue.empty()) {
496        NodeIt w=queue.front();
497        queue.pop();
498
499        for(InEdgeIt e=G.firstIn(w) ; e; G.next(e)) {
500          NodeIt v=G.tail(e);
501          if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
502            queue.push(v);
503            M.set(v, true);
504          }
505        }
506
507        for(OutEdgeIt e=G.firstOut(w) ; e; G.next(e)) {
508          NodeIt v=G.head(e);
509          if (!M.get(v) && flow.get(e) > 0 ) {
510            queue.push(v);
511            M.set(v, true);
512          }
513        }
514      }
515
516      for(EachNodeIt v=G.firstNode() ; v; G.next(v)) {
517        M.set(v, !M.get(v));
518      }
519
520    }
521
522
523
524    template<typename CutMap>
525    void minCut(CutMap& M) {
526      minMinCut(M);
527    }
528
529
530  };
531}//namespace marci
532#endif
533
534
535
536
Note: See TracBrowser for help on using the repository browser.