COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/hugo/max_flow.h @ 780:e06d0d16595f

Last change on this file since 780:e06d0d16595f was 774:4297098d9677, checked in by Alpar Juttner, 20 years ago

Merge back the whole branches/hugo++ to trunk.

File size: 23.6 KB
Line 
1// -*- C++ -*-
2#ifndef HUGO_MAX_FLOW_H
3#define HUGO_MAX_FLOW_H
4
5#include <vector>
6#include <queue>
7
8//#include <hugo/graph_wrapper.h>
9#include <hugo/invalid.h>
10#include <hugo/maps.h>
11
12/// \file
13/// \ingroup flowalgs
14
15namespace hugo {
16
17  /// \addtogroup flowalgs
18  /// @{                                                   
19
20  ///Maximum flow algorithms class.
21
22  ///This class provides various algorithms for finding a flow of
23  ///maximum value in a directed graph. The \e source node, the \e
24  ///target node, the \e capacity of the edges and the \e starting \e
25  ///flow value of the edges should be passed to the algorithm through the
26  ///constructor. It is possible to change these quantities using the
27  ///functions \ref setSource, \ref setTarget, \ref setCap and
28  ///\ref setFlow. Before any subsequent runs of any algorithm of
29  ///the class \ref setFlow should be called.
30  ///
31  ///After running an algorithm of the class, the actual flow value
32  ///can be obtained by calling \ref flowValue(). The minimum
33  ///value cut can be written into a \c node map of \c bools by
34  ///calling \ref minCut. (\ref minMinCut and \ref maxMinCut writes
35  ///the inclusionwise minimum and maximum of the minimum value
36  ///cuts, resp.)
37  ///
38  ///\param Graph The directed graph type the algorithm runs on.
39  ///\param Num The number type of the capacities and the flow values.
40  ///\param CapMap The capacity map type.
41  ///\param FlowMap The flow map type.
42  ///
43  ///\author Marton Makai, Jacint Szabo
44  template <typename Graph, typename Num,
45            typename CapMap=typename Graph::template EdgeMap<Num>,
46            typename FlowMap=typename Graph::template EdgeMap<Num> >
47  class MaxFlow {
48  protected:
49    typedef typename Graph::Node Node;
50    typedef typename Graph::NodeIt NodeIt;
51    typedef typename Graph::EdgeIt EdgeIt;
52    typedef typename Graph::OutEdgeIt OutEdgeIt;
53    typedef typename Graph::InEdgeIt InEdgeIt;
54
55    typedef typename std::vector<Node> VecFirst;
56    typedef typename Graph::template NodeMap<Node> NNMap;
57    typedef typename std::vector<Node> VecNode;
58
59    const Graph* g;
60    Node s;
61    Node t;
62    const CapMap* capacity;
63    FlowMap* flow;
64    int n;      //the number of nodes of G
65    //    typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;   
66    //typedef ExpResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
67    //    typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
68    //    typedef typename ResGW::Edge ResGWEdge;
69    typedef typename Graph::template NodeMap<int> ReachedMap;
70
71
72    //level works as a bool map in augmenting path algorithms and is
73    //used by bfs for storing reached information.  In preflow, it
74    //shows the levels of nodes.     
75    ReachedMap level;
76
77    //excess is needed only in preflow
78    typename Graph::template NodeMap<Num> excess;
79
80    // constants used for heuristics
81    static const int H0=20;
82    static const int H1=1;
83
84  public:
85
86    ///Indicates the property of the starting flow.
87
88    ///Indicates the property of the starting flow. The meanings are as follows:
89    ///- \c ZERO_FLOW: constant zero flow
90    ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to
91    ///the sum of the out-flows in every node except the \e source and
92    ///the \e target.
93    ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at
94    ///least the sum of the out-flows in every node except the \e source.
95    ///- \c NO_FLOW: indicates an unspecified edge map. \ref flow will be
96    ///set to the constant zero flow in the beginning of the algorithm in this case.
97    enum FlowEnum{
98      ZERO_FLOW,
99      GEN_FLOW,
100      PRE_FLOW,
101      NO_FLOW
102    };
103
104    enum StatusEnum {
105      AFTER_NOTHING,
106      AFTER_AUGMENTING,
107      AFTER_FAST_AUGMENTING,
108      AFTER_PRE_FLOW_PHASE_1,     
109      AFTER_PRE_FLOW_PHASE_2
110    };
111
112    /// Do not needle this flag only if necessary.
113    StatusEnum status;
114
115    //     int number_of_augmentations;
116
117
118    //     template<typename IntMap>
119    //     class TrickyReachedMap {
120    //     protected:
121    //       IntMap* map;
122    //       int* number_of_augmentations;
123    //     public:
124    //       TrickyReachedMap(IntMap& _map, int& _number_of_augmentations) :
125    //  map(&_map), number_of_augmentations(&_number_of_augmentations) { }
126    //       void set(const Node& n, bool b) {
127    //  if (b)
128    //    map->set(n, *number_of_augmentations);
129    //  else
130    //    map->set(n, *number_of_augmentations-1);
131    //       }
132    //       bool operator[](const Node& n) const {
133    //  return (*map)[n]==*number_of_augmentations;
134    //       }
135    //     };
136   
137    ///Constructor
138
139    ///\todo Document, please.
140    ///
141    MaxFlow(const Graph& _G, Node _s, Node _t,
142            const CapMap& _capacity, FlowMap& _flow) :
143      g(&_G), s(_s), t(_t), capacity(&_capacity),
144      flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0),
145      status(AFTER_NOTHING) { }
146
147    ///Runs a maximum flow algorithm.
148
149    ///Runs a preflow algorithm, which is the fastest maximum flow
150    ///algorithm up-to-date. The default for \c fe is ZERO_FLOW.
151    ///\pre The starting flow must be
152    /// - a constant zero flow if \c fe is \c ZERO_FLOW,
153    /// - an arbitary flow if \c fe is \c GEN_FLOW,
154    /// - an arbitary preflow if \c fe is \c PRE_FLOW,
155    /// - any map if \c fe is NO_FLOW.
156    void run(FlowEnum fe=ZERO_FLOW) {
157      preflow(fe);
158    }
159
160                                                                             
161    ///Runs a preflow algorithm. 
162
163    ///Runs a preflow algorithm. The preflow algorithms provide the
164    ///fastest way to compute a maximum flow in a directed graph.
165    ///\pre The starting flow must be
166    /// - a constant zero flow if \c fe is \c ZERO_FLOW,
167    /// - an arbitary flow if \c fe is \c GEN_FLOW,
168    /// - an arbitary preflow if \c fe is \c PRE_FLOW,
169    /// - any map if \c fe is NO_FLOW.
170    ///
171    ///\todo NO_FLOW should be the default flow.
172    void preflow(FlowEnum fe) {
173      preflowPhase1(fe);
174      preflowPhase2();
175    }
176    // Heuristics:
177    //   2 phase
178    //   gap
179    //   list 'level_list' on the nodes on level i implemented by hand
180    //   stack 'active' on the active nodes on level i                                                                                   
181    //   runs heuristic 'highest label' for H1*n relabels
182    //   runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
183    //   Parameters H0 and H1 are initialized to 20 and 1.
184
185    ///Runs the first phase of the preflow algorithm.
186
187    ///The preflow algorithm consists of two phases, this method runs the
188    ///first phase. After the first phase the maximum flow value and a
189    ///minimum value cut can already be computed, though a maximum flow
190    ///is not yet obtained. So after calling this method \ref flowValue
191    ///and \ref actMinCut gives proper results.
192    ///\warning: \ref minCut, \ref minMinCut and \ref maxMinCut do not
193    ///give minimum value cuts unless calling \ref preflowPhase2.
194    ///\pre The starting flow must be
195    /// - a constant zero flow if \c fe is \c ZERO_FLOW,
196    /// - an arbitary flow if \c fe is \c GEN_FLOW,
197    /// - an arbitary preflow if \c fe is \c PRE_FLOW,
198    /// - any map if \c fe is NO_FLOW.
199    void preflowPhase1(FlowEnum fe)
200    {
201
202      int heur0=(int)(H0*n);  //time while running 'bound decrease'
203      int heur1=(int)(H1*n);  //time while running 'highest label'
204      int heur=heur1;         //starting time interval (#of relabels)
205      int numrelabel=0;
206
207      bool what_heur=1;
208      //It is 0 in case 'bound decrease' and 1 in case 'highest label'
209
210      bool end=false;
211      //Needed for 'bound decrease', true means no active nodes are above bound
212      //b.
213
214      int k=n-2;  //bound on the highest level under n containing a node
215      int b=k;    //bound on the highest level under n of an active node
216
217      VecFirst first(n, INVALID);
218      NNMap next(*g, INVALID); //maybe INVALID is not needed
219
220      NNMap left(*g, INVALID);
221      NNMap right(*g, INVALID);
222      VecNode level_list(n,INVALID);
223      //List of the nodes in level i<n, set to n.
224
225      preflowPreproc(fe, next, first, level_list, left, right);
226      //End of preprocessing
227
228      //Push/relabel on the highest level active nodes.
229      while ( true ) {
230        if ( b == 0 ) {
231          if ( !what_heur && !end && k > 0 ) {
232            b=k;
233            end=true;
234          } else break;
235        }
236
237        if ( first[b]==INVALID ) --b;
238        else {
239          end=false;
240          Node w=first[b];
241          first[b]=next[w];
242          int newlevel=push(w, next, first);
243          if ( excess[w] > 0 ) relabel(w, newlevel, next, first, level_list,
244                                       left, right, b, k, what_heur);
245
246          ++numrelabel;
247          if ( numrelabel >= heur ) {
248            numrelabel=0;
249            if ( what_heur ) {
250              what_heur=0;
251              heur=heur0;
252              end=false;
253            } else {
254              what_heur=1;
255              heur=heur1;
256              b=k;
257            }
258          }
259        }
260      }
261
262      status=AFTER_PRE_FLOW_PHASE_1;
263    }
264
265
266    ///Runs the second phase of the preflow algorithm.
267
268    ///The preflow algorithm consists of two phases, this method runs
269    ///the second phase. After calling \ref preflowPhase1 and then
270    ///\ref preflowPhase2 the methods \ref flowValue, \ref minCut,
271    ///\ref minMinCut and \ref maxMinCut give proper results.
272    ///\pre \ref preflowPhase1 must be called before.
273    void preflowPhase2()
274    {
275
276      int k=n-2;  //bound on the highest level under n containing a node
277      int b=k;    //bound on the highest level under n of an active node
278
279   
280      VecFirst first(n, INVALID);
281      NNMap next(*g, INVALID); //maybe INVALID is not needed
282      level.set(s,0);
283      std::queue<Node> bfs_queue;
284      bfs_queue.push(s);
285
286      while (!bfs_queue.empty()) {
287
288        Node v=bfs_queue.front();
289        bfs_queue.pop();
290        int l=level[v]+1;
291
292        for(InEdgeIt e(*g,v); e!=INVALID; ++e) {
293          if ( (*capacity)[e] <= (*flow)[e] ) continue;
294          Node u=g->tail(e);
295          if ( level[u] >= n ) {
296            bfs_queue.push(u);
297            level.set(u, l);
298            if ( excess[u] > 0 ) {
299              next.set(u,first[l]);
300              first[l]=u;
301            }
302          }
303        }
304
305        for(OutEdgeIt e(*g,v); e!=INVALID; ++e) {
306          if ( 0 >= (*flow)[e] ) continue;
307          Node u=g->head(e);
308          if ( level[u] >= n ) {
309            bfs_queue.push(u);
310            level.set(u, l);
311            if ( excess[u] > 0 ) {
312              next.set(u,first[l]);
313              first[l]=u;
314            }
315          }
316        }
317      }
318      b=n-2;
319
320      while ( true ) {
321
322        if ( b == 0 ) break;
323
324        if ( first[b]==INVALID ) --b;
325        else {
326
327          Node w=first[b];
328          first[b]=next[w];
329          int newlevel=push(w,next, first/*active*/);
330
331          //relabel
332          if ( excess[w] > 0 ) {
333            level.set(w,++newlevel);
334            next.set(w,first[newlevel]);
335            first[newlevel]=w;
336            b=newlevel;
337          }
338        }
339      } // while(true)
340
341      status=AFTER_PRE_FLOW_PHASE_2;
342    }
343
344
345    /// Returns the value of the maximum flow.
346
347    /// Returns the excess of the target node \ref t.
348    /// After running \ref preflowPhase1, this is the value of
349    /// the maximum flow.
350    /// It can be called already after running \ref preflowPhase1.
351    Num flowValue() const {
352      //       Num a=0;
353      //       for(InEdgeIt e(*g,t);g->valid(e);g->next(e)) a+=(*flow)[e];
354      //       for(OutEdgeIt e(*g,t);g->valid(e);g->next(e)) a-=(*flow)[e];
355      //       return a;
356      return excess[t];
357      //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan   
358    }
359
360
361    ///Returns a minimum value cut after calling \ref preflowPhase1.
362
363    ///After the first phase of the preflow algorithm the maximum flow
364    ///value and a minimum value cut can already be computed. This
365    ///method can be called after running \ref preflowPhase1 for
366    ///obtaining a minimum value cut.
367    /// \warning Gives proper result only right after calling \ref
368    /// preflowPhase1.
369    /// \todo We have to make some status variable which shows the
370    /// actual state
371    /// of the class. This enables us to determine which methods are valid
372    /// for MinCut computation
373    template<typename _CutMap>
374    void actMinCut(_CutMap& M) const {
375      switch (status) {
376        case AFTER_PRE_FLOW_PHASE_1:
377        for(NodeIt v(*g); v!=INVALID; ++v) {
378          if (level[v] < n) {
379            M.set(v, false);
380          } else {
381            M.set(v, true);
382          }
383        }
384        break;
385        case AFTER_PRE_FLOW_PHASE_2:
386        case AFTER_NOTHING:
387        case AFTER_AUGMENTING:
388        case AFTER_FAST_AUGMENTING:
389        minMinCut(M);
390        break;
391      }
392    }
393
394    ///Returns the inclusionwise minimum of the minimum value cuts.
395
396    ///Sets \c M to the characteristic vector of the minimum value cut
397    ///which is inclusionwise minimum. It is computed by processing
398    ///a bfs from the source node \c s in the residual graph.
399    ///\pre M should be a node map of bools initialized to false.
400    ///\pre \c flow must be a maximum flow.
401    template<typename _CutMap>
402    void minMinCut(_CutMap& M) const {
403      std::queue<Node> queue;
404
405      M.set(s,true);
406      queue.push(s);
407
408      while (!queue.empty()) {
409        Node w=queue.front();
410        queue.pop();
411
412        for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
413          Node v=g->head(e);
414          if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
415            queue.push(v);
416            M.set(v, true);
417          }
418        }
419
420        for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
421          Node v=g->tail(e);
422          if (!M[v] && (*flow)[e] > 0 ) {
423            queue.push(v);
424            M.set(v, true);
425          }
426        }
427      }
428    }
429
430    ///Returns the inclusionwise maximum of the minimum value cuts.
431
432    ///Sets \c M to the characteristic vector of the minimum value cut
433    ///which is inclusionwise maximum. It is computed by processing a
434    ///backward bfs from the target node \c t in the residual graph.
435    ///\pre M should be a node map of bools initialized to false.
436    ///\pre \c flow must be a maximum flow.
437    template<typename _CutMap>
438    void maxMinCut(_CutMap& M) const {
439
440      for(NodeIt v(*g) ; v!=INVALID; ++v) M.set(v, true);
441
442      std::queue<Node> queue;
443
444      M.set(t,false);
445      queue.push(t);
446
447      while (!queue.empty()) {
448        Node w=queue.front();
449        queue.pop();
450
451        for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
452          Node v=g->tail(e);
453          if (M[v] && (*flow)[e] < (*capacity)[e] ) {
454            queue.push(v);
455            M.set(v, false);
456          }
457        }
458
459        for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
460          Node v=g->head(e);
461          if (M[v] && (*flow)[e] > 0 ) {
462            queue.push(v);
463            M.set(v, false);
464          }
465        }
466      }
467    }
468
469    ///Returns a minimum value cut.
470
471    ///Sets \c M to the characteristic vector of a minimum value cut.
472    ///\pre M should be a node map of bools initialized to false.
473    ///\pre \c flow must be a maximum flow.   
474    template<typename CutMap>
475    void minCut(CutMap& M) const { minMinCut(M); }
476
477    ///Sets the source node to \c _s.
478
479    ///Sets the source node to \c _s.
480    ///
481    void setSource(Node _s) { s=_s; status=AFTER_NOTHING; }
482
483    ///Sets the target node to \c _t.
484
485    ///Sets the target node to \c _t.
486    ///
487    void setTarget(Node _t) { t=_t; status=AFTER_NOTHING; }
488
489    /// Sets the edge map of the capacities to _cap.
490
491    /// Sets the edge map of the capacities to _cap.
492    ///
493    void setCap(const CapMap& _cap)
494    { capacity=&_cap; status=AFTER_NOTHING; }
495
496    /// Sets the edge map of the flows to _flow.
497
498    /// Sets the edge map of the flows to _flow.
499    ///
500    void setFlow(FlowMap& _flow) { flow=&_flow; status=AFTER_NOTHING; }
501
502
503  private:
504
505    int push(Node w, NNMap& next, VecFirst& first) {
506
507      int lev=level[w];
508      Num exc=excess[w];
509      int newlevel=n;       //bound on the next level of w
510
511      for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
512        if ( (*flow)[e] >= (*capacity)[e] ) continue;
513        Node v=g->head(e);
514
515        if( lev > level[v] ) { //Push is allowed now
516         
517          if ( excess[v]<=0 && v!=t && v!=s ) {
518            next.set(v,first[level[v]]);
519            first[level[v]]=v;
520          }
521
522          Num cap=(*capacity)[e];
523          Num flo=(*flow)[e];
524          Num remcap=cap-flo;
525         
526          if ( remcap >= exc ) { //A nonsaturating push.
527           
528            flow->set(e, flo+exc);
529            excess.set(v, excess[v]+exc);
530            exc=0;
531            break;
532
533          } else { //A saturating push.
534            flow->set(e, cap);
535            excess.set(v, excess[v]+remcap);
536            exc-=remcap;
537          }
538        } else if ( newlevel > level[v] ) newlevel = level[v];
539      } //for out edges wv
540
541      if ( exc > 0 ) {
542        for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
543         
544          if( (*flow)[e] <= 0 ) continue;
545          Node v=g->tail(e);
546
547          if( lev > level[v] ) { //Push is allowed now
548
549            if ( excess[v]<=0 && v!=t && v!=s ) {
550              next.set(v,first[level[v]]);
551              first[level[v]]=v;
552            }
553
554            Num flo=(*flow)[e];
555
556            if ( flo >= exc ) { //A nonsaturating push.
557
558              flow->set(e, flo-exc);
559              excess.set(v, excess[v]+exc);
560              exc=0;
561              break;
562            } else {  //A saturating push.
563
564              excess.set(v, excess[v]+flo);
565              exc-=flo;
566              flow->set(e,0);
567            }
568          } else if ( newlevel > level[v] ) newlevel = level[v];
569        } //for in edges vw
570
571      } // if w still has excess after the out edge for cycle
572
573      excess.set(w, exc);
574     
575      return newlevel;
576    }
577   
578   
579   
580    void preflowPreproc(FlowEnum fe, NNMap& next, VecFirst& first,
581                        VecNode& level_list, NNMap& left, NNMap& right)
582    {
583      switch (fe) {  //setting excess
584        case NO_FLOW:
585        for(EdgeIt e(*g); e!=INVALID; ++e) flow->set(e,0);
586        for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0);
587        break;
588        case ZERO_FLOW:
589        for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0);
590        break;
591        case GEN_FLOW:
592        for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0);
593        {
594          Num exc=0;
595          for(InEdgeIt e(*g,t) ; e!=INVALID; ++e) exc+=(*flow)[e];
596          for(OutEdgeIt e(*g,t) ; e!=INVALID; ++e) exc-=(*flow)[e];
597          excess.set(t,exc);
598        }
599        break;
600        default:
601        break;
602      }
603     
604      for(NodeIt v(*g); v!=INVALID; ++v) level.set(v,n);
605      //setting each node to level n
606     
607      std::queue<Node> bfs_queue;
608
609
610      switch (fe) {
611      case NO_FLOW:   //flow is already set to const zero
612      case ZERO_FLOW:
613        //Reverse_bfs from t, to find the starting level.
614        level.set(t,0);
615        bfs_queue.push(t);
616       
617        while (!bfs_queue.empty()) {
618         
619          Node v=bfs_queue.front();
620          bfs_queue.pop();
621          int l=level[v]+1;
622         
623          for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) {
624            Node w=g->tail(e);
625            if ( level[w] == n && w != s ) {
626              bfs_queue.push(w);
627              Node z=level_list[l];
628              if ( z!=INVALID ) left.set(z,w);
629              right.set(w,z);
630              level_list[l]=w;
631              level.set(w, l);
632            }
633          }
634        }
635       
636        //the starting flow
637        for(OutEdgeIt e(*g,s) ; e!=INVALID; ++e)
638          {
639            Num c=(*capacity)[e];
640            if ( c <= 0 ) continue;
641            Node w=g->head(e);
642            if ( level[w] < n ) {
643              if ( excess[w] <= 0 && w!=t ) //putting into the stack
644                {
645                  next.set(w,first[level[w]]);
646                  first[level[w]]=w;
647                }
648              flow->set(e, c);
649              excess.set(w, excess[w]+c);
650            }
651          }
652        break;
653      case GEN_FLOW:
654        //Reverse_bfs from t in the residual graph,
655        //to find the starting level.
656        level.set(t,0);
657        bfs_queue.push(t);
658       
659        while (!bfs_queue.empty()) {
660         
661          Node v=bfs_queue.front();
662          bfs_queue.pop();
663          int l=level[v]+1;
664         
665          for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) {
666            if ( (*capacity)[e] <= (*flow)[e] ) continue;
667            Node w=g->tail(e);
668            if ( level[w] == n && w != s ) {
669              bfs_queue.push(w);
670              Node z=level_list[l];
671              if ( z!=INVALID ) left.set(z,w);
672              right.set(w,z);
673              level_list[l]=w;
674              level.set(w, l);
675            }
676          }
677         
678          for(OutEdgeIt e(*g,v) ; e!=INVALID; ++e) {
679            if ( 0 >= (*flow)[e] ) continue;
680            Node w=g->head(e);
681            if ( level[w] == n && w != s ) {
682              bfs_queue.push(w);
683              Node z=level_list[l];
684              if ( z!=INVALID ) left.set(z,w);
685              right.set(w,z);
686              level_list[l]=w;
687              level.set(w, l);
688            }
689          }
690        }
691       
692        //the starting flow
693        for(OutEdgeIt e(*g,s); e!=INVALID; ++e)
694          {
695            Num rem=(*capacity)[e]-(*flow)[e];
696            if ( rem <= 0 ) continue;
697            Node w=g->head(e);
698            if ( level[w] < n ) {
699              if ( excess[w] <= 0 && w!=t ) //putting into the stack
700                {
701                  next.set(w,first[level[w]]);
702                  first[level[w]]=w;
703                }   
704              flow->set(e, (*capacity)[e]);
705              excess.set(w, excess[w]+rem);
706            }
707          }
708       
709        for(InEdgeIt e(*g,s); e!=INVALID; ++e)
710          {
711            if ( (*flow)[e] <= 0 ) continue;
712            Node w=g->tail(e);
713            if ( level[w] < n ) {
714              if ( excess[w] <= 0 && w!=t )
715                {
716                  next.set(w,first[level[w]]);
717                  first[level[w]]=w;
718                } 
719              excess.set(w, excess[w]+(*flow)[e]);
720              flow->set(e, 0);
721            }
722          }
723        break;
724      case PRE_FLOW:
725        //Reverse_bfs from t in the residual graph,
726        //to find the starting level.
727        level.set(t,0);
728        bfs_queue.push(t);
729       
730        while (!bfs_queue.empty()) {
731         
732          Node v=bfs_queue.front();
733          bfs_queue.pop();
734          int l=level[v]+1;
735         
736          for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) {
737            if ( (*capacity)[e] <= (*flow)[e] ) continue;
738            Node w=g->tail(e);
739            if ( level[w] == n && w != s ) {
740              bfs_queue.push(w);
741              Node z=level_list[l];
742              if ( z!=INVALID ) left.set(z,w);
743              right.set(w,z);
744              level_list[l]=w;
745              level.set(w, l);
746            }
747          }
748         
749          for(OutEdgeIt e(*g,v) ; e!=INVALID; ++e) {
750            if ( 0 >= (*flow)[e] ) continue;
751            Node w=g->head(e);
752            if ( level[w] == n && w != s ) {
753              bfs_queue.push(w);
754              Node z=level_list[l];
755              if ( z!=INVALID ) left.set(z,w);
756              right.set(w,z);
757              level_list[l]=w;
758              level.set(w, l);
759            }
760          }
761        }
762       
763       
764        //the starting flow
765        for(OutEdgeIt e(*g,s) ; e!=INVALID; ++e) {
766          Num rem=(*capacity)[e]-(*flow)[e];
767          if ( rem <= 0 ) continue;
768          Node w=g->head(e);
769          if ( level[w] < n ) {
770            flow->set(e, (*capacity)[e]);
771            excess.set(w, excess[w]+rem);
772          }
773        }
774       
775        for(InEdgeIt e(*g,s) ; e!=INVALID; ++e) {
776          if ( (*flow)[e] <= 0 ) continue;
777          Node w=g->tail(e);
778          if ( level[w] < n ) {
779            excess.set(w, excess[w]+(*flow)[e]);
780            flow->set(e, 0);
781          }
782        }
783       
784        //computing the excess
785        for(NodeIt w(*g); w!=INVALID; ++w) {
786          Num exc=0;
787         
788          for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) exc+=(*flow)[e];
789          for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) exc-=(*flow)[e];
790         
791          excess.set(w,exc);
792         
793          //putting the active nodes into the stack
794          int lev=level[w];
795            if ( exc > 0 && lev < n && Node(w) != t )
796              ///\bug       if ( exc > 0 && lev < n && w != t ) temporarily for working with wrappers.
797            {
798              next.set(w,first[lev]);
799              first[lev]=w;
800            }
801        }
802        break;
803      } //switch
804    } //preflowPreproc
805
806
807    void relabel(Node w, int newlevel, NNMap& next, VecFirst& first,
808                 VecNode& level_list, NNMap& left,
809                 NNMap& right, int& b, int& k, bool what_heur )
810    {
811
812      int lev=level[w];
813
814      Node right_n=right[w];
815      Node left_n=left[w];
816
817      //unlacing starts
818      if ( right_n!=INVALID ) {
819        if ( left_n!=INVALID ) {
820          right.set(left_n, right_n);
821          left.set(right_n, left_n);
822        } else {
823          level_list[lev]=right_n;
824          left.set(right_n, INVALID);
825        }
826      } else {
827        if ( left_n!=INVALID ) {
828          right.set(left_n, INVALID);
829        } else {
830          level_list[lev]=INVALID;
831        }
832      }
833      //unlacing ends
834
835      if ( level_list[lev]==INVALID ) {
836
837        //gapping starts
838        for (int i=lev; i!=k ; ) {
839          Node v=level_list[++i];
840          while ( v!=INVALID ) {
841            level.set(v,n);
842            v=right[v];
843          }
844          level_list[i]=INVALID;
845          if ( !what_heur ) first[i]=INVALID;
846        }
847
848        level.set(w,n);
849        b=lev-1;
850        k=b;
851        //gapping ends
852
853      } else {
854
855        if ( newlevel == n ) level.set(w,n);
856        else {
857          level.set(w,++newlevel);
858          next.set(w,first[newlevel]);
859          first[newlevel]=w;
860          if ( what_heur ) b=newlevel;
861          if ( k < newlevel ) ++k;      //now k=newlevel
862          Node z=level_list[newlevel];
863          if ( z!=INVALID ) left.set(z,w);
864          right.set(w,z);
865          left.set(w,INVALID);
866          level_list[newlevel]=w;
867        }
868      }
869    } //relabel
870
871    void printexcess() {////
872      std::cout << "Excesses:" <<std::endl;
873
874      for(NodeIt v(*g); v!=INVALID ; ++v) {
875        std::cout << 1+(g->id(v)) << ":" << excess[v]<<std::endl;
876      }
877    }
878
879    void printlevel() {////
880      std::cout << "Levels:" <<std::endl;
881
882      for(NodeIt v(*g); v!=INVALID ; ++v) {
883        std::cout << 1+(g->id(v)) << ":" << level[v]<<std::endl;
884      }
885    }
886
887    void printactive() {////
888      std::cout << "Levels:" <<std::endl;
889
890      for(NodeIt v(*g); v!=INVALID ; ++v) {
891        std::cout << 1+(g->id(v)) << ":" << level[v]<<std::endl;
892      }
893    }
894
895
896  };  //class MaxFlow
897} //namespace hugo
898
899#endif //HUGO_MAX_FLOW_H
900
901
902
903
Note: See TracBrowser for help on using the repository browser.