COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/max_flow_no_stack.h @ 719:cb9efd4cc9db

Last change on this file since 719:cb9efd4cc9db was 719:cb9efd4cc9db, checked in by Alpar Juttner, 16 years ago

Indenting

File size: 35.8 KB
Line 
1// -*- C++ -*-
2#ifndef HUGO_MAX_FLOW_NO_STACK_H
3#define HUGO_MAX_FLOW_NO_STACK_H
4
5#include <vector>
6#include <queue>
7//#include <stack>
8
9#include <hugo/graph_wrapper.h>
10#include <bfs_dfs.h>
11#include <hugo/invalid.h>
12#include <hugo/maps.h>
13#include <hugo/for_each_macros.h>
14
15/// \file
16/// \brief The same as max_flow.h, but without using stl stack for the active nodes. Only for test.
17/// \ingroup galgs
18
19namespace hugo {
20
21  /// \addtogroup galgs
22  /// @{                                                                                                                                       
23  ///Maximum flow algorithms class.
24
25  ///This class provides various algorithms for finding a flow of
26  ///maximum value in a directed graph. The \e source node, the \e
27  ///target node, the \e capacity of the edges and the \e starting \e
28  ///flow value of the edges should be passed to the algorithm through the
29  ///constructor. It is possible to change these quantities using the
30  ///functions \ref resetSource, \ref resetTarget, \ref resetCap and
31  ///\ref resetFlow. Before any subsequent runs of any algorithm of
32  ///the class \ref resetFlow should be called.
33
34  ///After running an algorithm of the class, the actual flow value
35  ///can be obtained by calling \ref flowValue(). The minimum
36  ///value cut can be written into a \c node map of \c bools by
37  ///calling \ref minCut. (\ref minMinCut and \ref maxMinCut writes
38  ///the inclusionwise minimum and maximum of the minimum value
39  ///cuts, resp.)                                                                                                                               
40  ///\param Graph The directed graph type the algorithm runs on.
41  ///\param Num The number type of the capacities and the flow values.
42  ///\param CapMap The capacity map type.
43  ///\param FlowMap The flow map type.                                                                                                           
44  ///\author Marton Makai, Jacint Szabo
45  template <typename Graph, typename Num,
46            typename CapMap=typename Graph::template EdgeMap<Num>,
47            typename FlowMap=typename Graph::template EdgeMap<Num> >
48  class MaxFlowNoStack {
49  protected:
50    typedef typename Graph::Node Node;
51    typedef typename Graph::NodeIt NodeIt;
52    typedef typename Graph::EdgeIt EdgeIt;
53    typedef typename Graph::OutEdgeIt OutEdgeIt;
54    typedef typename Graph::InEdgeIt InEdgeIt;
55
56    //    typedef typename std::vector<std::stack<Node> > VecStack;
57    typedef typename std::vector<Node> VecFirst;
58    typedef typename Graph::template NodeMap<Node> NNMap;
59    typedef typename std::vector<Node> VecNode;
60
61    const Graph* g;
62    Node s;
63    Node t;
64    const CapMap* capacity;
65    FlowMap* flow;
66    int n;      //the number of nodes of G
67    typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;   
68    //typedef ExpResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
69    typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
70    typedef typename ResGW::Edge ResGWEdge;
71    //typedef typename ResGW::template NodeMap<bool> ReachedMap;
72    typedef typename Graph::template NodeMap<int> ReachedMap;
73
74
75    //level works as a bool map in augmenting path algorithms and is
76    //used by bfs for storing reached information.  In preflow, it
77    //shows the levels of nodes.     
78    ReachedMap level;
79
80    //excess is needed only in preflow
81    typename Graph::template NodeMap<Num> excess;
82
83    //fixme   
84//   protected:
85    //     MaxFlow() { }
86    //     void set(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,
87    //       FlowMap& _flow)
88    //       {
89    //  g=&_G;
90    //  s=_s;
91    //  t=_t;
92    //  capacity=&_capacity;
93    //  flow=&_flow;
94    //  n=_G.nodeNum;
95    //  level.set (_G); //kellene vmi ilyesmi fv
96    //  excess(_G,0); //itt is
97    //       }
98
99    // constants used for heuristics
100    static const int H0=20;
101    static const int H1=1;
102
103  public:
104
105    ///Indicates the property of the starting flow.
106
107    ///Indicates the property of the starting flow. The meanings are as follows:
108    ///- \c ZERO_FLOW: constant zero flow
109    ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to
110    ///the sum of the out-flows in every node except the \e source and
111    ///the \e target.
112    ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at
113    ///least the sum of the out-flows in every node except the \e source.
114    ///- \c NO_FLOW: indicates an unspecified edge map. \ref flow will be
115    ///set to the constant zero flow in the beginning of the algorithm in this case.
116    enum FlowEnum{
117      ZERO_FLOW,
118      GEN_FLOW,
119      PRE_FLOW,
120      NO_FLOW
121    };
122
123    enum StatusEnum {
124      AFTER_NOTHING,
125      AFTER_AUGMENTING,
126      AFTER_FAST_AUGMENTING,
127      AFTER_PRE_FLOW_PHASE_1,     
128      AFTER_PRE_FLOW_PHASE_2
129    };
130
131    /// Don not needle this flag only if necessary.
132    StatusEnum status;
133    int number_of_augmentations;
134
135
136    template<typename IntMap>
137    class TrickyReachedMap {
138    protected:
139      IntMap* map;
140      int* number_of_augmentations;
141    public:
142      TrickyReachedMap(IntMap& _map, int& _number_of_augmentations) :
143        map(&_map), number_of_augmentations(&_number_of_augmentations) { }
144      void set(const Node& n, bool b) {
145        if (b)
146          map->set(n, *number_of_augmentations);
147        else
148          map->set(n, *number_of_augmentations-1);
149      }
150      bool operator[](const Node& n) const {
151        return (*map)[n]==*number_of_augmentations;
152      }
153    };
154   
155    ///Constructor
156
157    ///\todo Document, please.
158    ///
159    MaxFlowNoStack(const Graph& _G, Node _s, Node _t,
160                   const CapMap& _capacity, FlowMap& _flow) :
161      g(&_G), s(_s), t(_t), capacity(&_capacity),
162      flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0),
163      status(AFTER_NOTHING), number_of_augmentations(0) { }
164
165    ///Runs a maximum flow algorithm.
166
167    ///Runs a preflow algorithm, which is the fastest maximum flow
168    ///algorithm up-to-date. The default for \c fe is ZERO_FLOW.
169    ///\pre The starting flow must be
170    /// - a constant zero flow if \c fe is \c ZERO_FLOW,
171    /// - an arbitary flow if \c fe is \c GEN_FLOW,
172    /// - an arbitary preflow if \c fe is \c PRE_FLOW,
173    /// - any map if \c fe is NO_FLOW.
174    void run(FlowEnum fe=ZERO_FLOW) {
175      preflow(fe);
176    }
177
178                                                                             
179    ///Runs a preflow algorithm. 
180
181    ///Runs a preflow algorithm. The preflow algorithms provide the
182    ///fastest way to compute a maximum flow in a directed graph.
183    ///\pre The starting flow must be
184    /// - a constant zero flow if \c fe is \c ZERO_FLOW,
185    /// - an arbitary flow if \c fe is \c GEN_FLOW,
186    /// - an arbitary preflow if \c fe is \c PRE_FLOW,
187    /// - any map if \c fe is NO_FLOW.
188    ///
189    ///\todo NO_FLOW should be the default flow.
190    void preflow(FlowEnum fe) {
191      preflowPhase1(fe);
192      preflowPhase2();
193    }
194    // Heuristics:
195    //   2 phase
196    //   gap
197    //   list 'level_list' on the nodes on level i implemented by hand
198    //   stack 'active' on the active nodes on level i                                                                                   
199    //   runs heuristic 'highest label' for H1*n relabels
200    //   runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
201    //   Parameters H0 and H1 are initialized to 20 and 1.
202
203    ///Runs the first phase of the preflow algorithm.
204
205    ///The preflow algorithm consists of two phases, this method runs the
206    ///first phase. After the first phase the maximum flow value and a
207    ///minimum value cut can already be computed, though a maximum flow
208    ///is net yet obtained. So after calling this method \ref flowValue
209    ///and \ref actMinCut gives proper results.
210    ///\warning: \ref minCut, \ref minMinCut and \ref maxMinCut do not
211    ///give minimum value cuts unless calling \ref preflowPhase2.
212    ///\pre The starting flow must be
213    /// - a constant zero flow if \c fe is \c ZERO_FLOW,
214    /// - an arbitary flow if \c fe is \c GEN_FLOW,
215    /// - an arbitary preflow if \c fe is \c PRE_FLOW,
216    /// - any map if \c fe is NO_FLOW.
217    void preflowPhase1(FlowEnum fe);
218
219    ///Runs the second phase of the preflow algorithm.
220
221    ///The preflow algorithm consists of two phases, this method runs
222    ///the second phase. After calling \ref preflowPhase1 and then
223    ///\ref preflowPhase2 the methods \ref flowValue, \ref minCut,
224    ///\ref minMinCut and \ref maxMinCut give proper results.
225    ///\pre \ref preflowPhase1 must be called before.
226    void preflowPhase2();
227
228    /// Starting from a flow, this method searches for an augmenting path
229    /// according to the Edmonds-Karp algorithm
230    /// and augments the flow on if any.
231    /// The return value shows if the augmentation was succesful.
232    bool augmentOnShortestPath();
233    bool augmentOnShortestPath2();
234
235    /// Starting from a flow, this method searches for an augmenting blocking
236    /// flow according to Dinits' algorithm and augments the flow on if any.
237    /// The blocking flow is computed in a physically constructed
238    /// residual graph of type \c Mutablegraph.
239    /// The return value show sif the augmentation was succesful.
240    template<typename MutableGraph> bool augmentOnBlockingFlow();
241
242    /// The same as \c augmentOnBlockingFlow<MutableGraph> but the
243    /// residual graph is not constructed physically.
244    /// The return value shows if the augmentation was succesful.
245    bool augmentOnBlockingFlow2();
246
247    /// Returns the maximum value of a flow.
248
249    /// Returns the maximum value of a flow, by counting the
250    /// over-flow of the target node \ref t.
251    /// It can be called already after running \ref preflowPhase1.
252    Num flowValue() const {
253      Num a=0;
254      FOR_EACH_INC_LOC(InEdgeIt, e, *g, t) a+=(*flow)[e];
255      FOR_EACH_INC_LOC(OutEdgeIt, e, *g, t) a-=(*flow)[e];
256      return a;
257      //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan   
258    }
259
260    ///Returns a minimum value cut after calling \ref preflowPhase1.
261
262    ///After the first phase of the preflow algorithm the maximum flow
263    ///value and a minimum value cut can already be computed. This
264    ///method can be called after running \ref preflowPhase1 for
265    ///obtaining a minimum value cut.
266    /// \warning Gives proper result only right after calling \ref
267    /// preflowPhase1.
268    /// \todo We have to make some status variable which shows the
269    /// actual state
270    /// of the class. This enables us to determine which methods are valid
271    /// for MinCut computation
272    template<typename _CutMap>
273    void actMinCut(_CutMap& M) const {
274      NodeIt v;
275      switch (status) {
276      case AFTER_PRE_FLOW_PHASE_1:
277        for(g->first(v); g->valid(v); g->next(v)) {
278          if (level[v] < n) {
279            M.set(v, false);
280          } else {
281            M.set(v, true);
282          }
283        }
284        break;
285      case AFTER_PRE_FLOW_PHASE_2:
286      case AFTER_NOTHING:
287        minMinCut(M);
288        break;
289      case AFTER_AUGMENTING:
290        for(g->first(v); g->valid(v); g->next(v)) {
291          if (level[v]) {
292            M.set(v, true);
293          } else {
294            M.set(v, false);
295          }
296        }
297        break;
298      case AFTER_FAST_AUGMENTING:
299        for(g->first(v); g->valid(v); g->next(v)) {
300          if (level[v]==number_of_augmentations) {
301            M.set(v, true);
302          } else {
303            M.set(v, false);
304          }
305        }
306        break;
307      }
308    }
309
310    ///Returns the inclusionwise minimum of the minimum value cuts.
311
312    ///Sets \c M to the characteristic vector of the minimum value cut
313    ///which is inclusionwise minimum. It is computed by processing
314    ///a bfs from the source node \c s in the residual graph.
315    ///\pre M should be a node map of bools initialized to false.
316    ///\pre \c flow must be a maximum flow.
317    template<typename _CutMap>
318    void minMinCut(_CutMap& M) const {
319      std::queue<Node> queue;
320
321      M.set(s,true);
322      queue.push(s);
323
324      while (!queue.empty()) {
325        Node w=queue.front();
326        queue.pop();
327
328        OutEdgeIt e;
329        for(g->first(e,w) ; g->valid(e); g->next(e)) {
330          Node v=g->head(e);
331          if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
332            queue.push(v);
333            M.set(v, true);
334          }
335        }
336
337        InEdgeIt f;
338        for(g->first(f,w) ; g->valid(f); g->next(f)) {
339          Node v=g->tail(f);
340          if (!M[v] && (*flow)[f] > 0 ) {
341            queue.push(v);
342            M.set(v, true);
343          }
344        }
345      }
346    }
347
348    ///Returns the inclusionwise maximum of the minimum value cuts.
349
350    ///Sets \c M to the characteristic vector of the minimum value cut
351    ///which is inclusionwise maximum. It is computed by processing a
352    ///backward bfs from the target node \c t in the residual graph.
353    ///\pre M should be a node map of bools initialized to false.
354    ///\pre \c flow must be a maximum flow.
355    template<typename _CutMap>
356    void maxMinCut(_CutMap& M) const {
357
358      NodeIt v;
359      for(g->first(v) ; g->valid(v); g->next(v)) {
360        M.set(v, true);
361      }
362
363      std::queue<Node> queue;
364
365      M.set(t,false);
366      queue.push(t);
367
368      while (!queue.empty()) {
369        Node w=queue.front();
370        queue.pop();
371
372        InEdgeIt e;
373        for(g->first(e,w) ; g->valid(e); g->next(e)) {
374          Node v=g->tail(e);
375          if (M[v] && (*flow)[e] < (*capacity)[e] ) {
376            queue.push(v);
377            M.set(v, false);
378          }
379        }
380
381        OutEdgeIt f;
382        for(g->first(f,w) ; g->valid(f); g->next(f)) {
383          Node v=g->head(f);
384          if (M[v] && (*flow)[f] > 0 ) {
385            queue.push(v);
386            M.set(v, false);
387          }
388        }
389      }
390    }
391
392    ///Returns a minimum value cut.
393
394    ///Sets \c M to the characteristic vector of a minimum value cut.
395    ///\pre M should be a node map of bools initialized to false.
396    ///\pre \c flow must be a maximum flow.   
397    template<typename CutMap>
398    void minCut(CutMap& M) const { minMinCut(M); }
399
400    ///Resets the source node to \c _s.
401
402    ///Resets the source node to \c _s.
403    ///
404    void resetSource(Node _s) { s=_s; status=AFTER_NOTHING; }
405
406    ///Resets the target node to \c _t.
407
408    ///Resets the target node to \c _t.
409    ///
410    void resetTarget(Node _t) { t=_t; status=AFTER_NOTHING; }
411
412    /// Resets the edge map of the capacities to _cap.
413
414    /// Resets the edge map of the capacities to _cap.
415    ///
416    void resetCap(const CapMap& _cap)
417    { capacity=&_cap; status=AFTER_NOTHING; }
418
419    /// Resets the edge map of the flows to _flow.
420
421    /// Resets the edge map of the flows to _flow.
422    ///
423    void resetFlow(FlowMap& _flow) { flow=&_flow; status=AFTER_NOTHING; }
424
425
426  private:
427
428    int push(Node w, NNMap& next, VecFirst& first) {
429
430      int lev=level[w];
431      Num exc=excess[w];
432      int newlevel=n;       //bound on the next level of w
433
434      OutEdgeIt e;
435      for(g->first(e,w); g->valid(e); g->next(e)) {
436
437        if ( (*flow)[e] >= (*capacity)[e] ) continue;
438        Node v=g->head(e);
439
440        if( lev > level[v] ) { //Push is allowed now
441
442          if ( excess[v]<=0 && v!=t && v!=s ) {
443            next.set(v,first[level[v]]);
444            first[level[v]]=v;
445            //      int lev_v=level[v];
446            //active[lev_v].push(v);
447          }
448
449          Num cap=(*capacity)[e];
450          Num flo=(*flow)[e];
451          Num remcap=cap-flo;
452
453          if ( remcap >= exc ) { //A nonsaturating push.
454
455            flow->set(e, flo+exc);
456            excess.set(v, excess[v]+exc);
457            exc=0;
458            break;
459
460          } else { //A saturating push.
461            flow->set(e, cap);
462            excess.set(v, excess[v]+remcap);
463            exc-=remcap;
464          }
465        } else if ( newlevel > level[v] ) newlevel = level[v];
466      } //for out edges wv
467
468      if ( exc > 0 ) {
469        InEdgeIt e;
470        for(g->first(e,w); g->valid(e); g->next(e)) {
471
472          if( (*flow)[e] <= 0 ) continue;
473          Node v=g->tail(e);
474
475          if( lev > level[v] ) { //Push is allowed now
476
477            if ( excess[v]<=0 && v!=t && v!=s ) {
478              next.set(v,first[level[v]]);
479              first[level[v]]=v;
480              //int lev_v=level[v];
481              //active[lev_v].push(v);
482            }
483
484            Num flo=(*flow)[e];
485
486            if ( flo >= exc ) { //A nonsaturating push.
487
488              flow->set(e, flo-exc);
489              excess.set(v, excess[v]+exc);
490              exc=0;
491              break;
492            } else {  //A saturating push.
493
494              excess.set(v, excess[v]+flo);
495              exc-=flo;
496              flow->set(e,0);
497            }
498          } else if ( newlevel > level[v] ) newlevel = level[v];
499        } //for in edges vw
500
501      } // if w still has excess after the out edge for cycle
502
503      excess.set(w, exc);
504
505      return newlevel;
506    }
507
508
509    void preflowPreproc(FlowEnum fe, NNMap& next, VecFirst& first,
510                        VecNode& level_list, NNMap& left, NNMap& right)
511    {
512      std::queue<Node> bfs_queue;
513
514      switch (fe) {
515      case NO_FLOW:   //flow is already set to const zero in this case
516      case ZERO_FLOW:
517        {
518          //Reverse_bfs from t, to find the starting level.
519          level.set(t,0);
520          bfs_queue.push(t);
521
522          while (!bfs_queue.empty()) {
523
524            Node v=bfs_queue.front();
525            bfs_queue.pop();
526            int l=level[v]+1;
527
528            InEdgeIt e;
529            for(g->first(e,v); g->valid(e); g->next(e)) {
530              Node w=g->tail(e);
531              if ( level[w] == n && w != s ) {
532                bfs_queue.push(w);
533                Node z=level_list[l];
534                if ( g->valid(z) ) left.set(z,w);
535                right.set(w,z);
536                level_list[l]=w;
537                level.set(w, l);
538              }
539            }
540          }
541
542          //the starting flow
543          OutEdgeIt e;
544          for(g->first(e,s); g->valid(e); g->next(e))
545            {
546              Num c=(*capacity)[e];
547              if ( c <= 0 ) continue;
548              Node w=g->head(e);
549              if ( level[w] < n ) {
550                if ( excess[w] <= 0 && w!=t )
551                  {
552                    next.set(w,first[level[w]]);
553                    first[level[w]]=w;
554                    //active[level[w]].push(w);
555                  }
556                flow->set(e, c);
557                excess.set(w, excess[w]+c);
558              }
559            }
560          break;
561        }
562
563      case GEN_FLOW:
564      case PRE_FLOW:
565        {
566          //Reverse_bfs from t in the residual graph,
567          //to find the starting level.
568          level.set(t,0);
569          bfs_queue.push(t);
570
571          while (!bfs_queue.empty()) {
572
573            Node v=bfs_queue.front();
574            bfs_queue.pop();
575            int l=level[v]+1;
576
577            InEdgeIt e;
578            for(g->first(e,v); g->valid(e); g->next(e)) {
579              if ( (*capacity)[e] <= (*flow)[e] ) continue;
580              Node w=g->tail(e);
581              if ( level[w] == n && w != s ) {
582                bfs_queue.push(w);
583                Node z=level_list[l];
584                if ( g->valid(z) ) left.set(z,w);
585                right.set(w,z);
586                level_list[l]=w;
587                level.set(w, l);
588              }
589            }
590
591            OutEdgeIt f;
592            for(g->first(f,v); g->valid(f); g->next(f)) {
593              if ( 0 >= (*flow)[f] ) continue;
594              Node w=g->head(f);
595              if ( level[w] == n && w != s ) {
596                bfs_queue.push(w);
597                Node z=level_list[l];
598                if ( g->valid(z) ) left.set(z,w);
599                right.set(w,z);
600                level_list[l]=w;
601                level.set(w, l);
602              }
603            }
604          }
605
606
607          //the starting flow
608          OutEdgeIt e;
609          for(g->first(e,s); g->valid(e); g->next(e))
610            {
611              Num rem=(*capacity)[e]-(*flow)[e];
612              if ( rem <= 0 ) continue;
613              Node w=g->head(e);
614              if ( level[w] < n ) {
615                if ( excess[w] <= 0 && w!=t )
616                  {
617                    next.set(w,first[level[w]]);
618                    first[level[w]]=w;
619                    //active[level[w]].push(w);
620                  }   
621                flow->set(e, (*capacity)[e]);
622                excess.set(w, excess[w]+rem);
623              }
624            }
625
626          InEdgeIt f;
627          for(g->first(f,s); g->valid(f); g->next(f))
628            {
629              if ( (*flow)[f] <= 0 ) continue;
630              Node w=g->tail(f);
631              if ( level[w] < n ) {
632                if ( excess[w] <= 0 && w!=t )
633                  {
634                    next.set(w,first[level[w]]);
635                    first[level[w]]=w;
636                    //active[level[w]].push(w);
637                  }   
638                excess.set(w, excess[w]+(*flow)[f]);
639                flow->set(f, 0);
640              }
641            }
642          break;
643        } //case PRE_FLOW
644      }
645    } //preflowPreproc
646
647
648
649    void relabel(Node w, int newlevel, NNMap& next, VecFirst& first,
650                 VecNode& level_list, NNMap& left,
651                 NNMap& right, int& b, int& k, bool what_heur )
652    {
653
654      Num lev=level[w];
655
656      Node right_n=right[w];
657      Node left_n=left[w];
658
659      //unlacing starts
660      if ( g->valid(right_n) ) {
661        if ( g->valid(left_n) ) {
662          right.set(left_n, right_n);
663          left.set(right_n, left_n);
664        } else {
665          level_list[lev]=right_n;
666          left.set(right_n, INVALID);
667        }
668      } else {
669        if ( g->valid(left_n) ) {
670          right.set(left_n, INVALID);
671        } else {
672          level_list[lev]=INVALID;
673        }
674      }
675      //unlacing ends
676
677      if ( !g->valid(level_list[lev]) ) {
678
679        //gapping starts
680        for (int i=lev; i!=k ; ) {
681          Node v=level_list[++i];
682          while ( g->valid(v) ) {
683            level.set(v,n);
684            v=right[v];
685          }
686          level_list[i]=INVALID;
687          if ( !what_heur ) first[i]=INVALID;
688          /*{
689            while ( !active[i].empty() ) {
690            active[i].pop();    //FIXME: ezt szebben kene
691            }
692            }*/
693        }
694
695        level.set(w,n);
696        b=lev-1;
697        k=b;
698        //gapping ends
699
700      } else {
701
702        if ( newlevel == n ) level.set(w,n);
703        else {
704          level.set(w,++newlevel);
705          next.set(w,first[newlevel]);
706          first[newlevel]=w;
707          //      active[newlevel].push(w);
708          if ( what_heur ) b=newlevel;
709          if ( k < newlevel ) ++k;      //now k=newlevel
710          Node z=level_list[newlevel];
711          if ( g->valid(z) ) left.set(z,w);
712          right.set(w,z);
713          left.set(w,INVALID);
714          level_list[newlevel]=w;
715        }
716      }
717
718    } //relabel
719
720
721    template<typename MapGraphWrapper>
722    class DistanceMap {
723    protected:
724      const MapGraphWrapper* g;
725      typename MapGraphWrapper::template NodeMap<int> dist;
726    public:
727      DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { }
728      void set(const typename MapGraphWrapper::Node& n, int a) {
729        dist.set(n, a);
730      }
731      int operator[](const typename MapGraphWrapper::Node& n) const {
732        return dist[n];
733      }
734      //       int get(const typename MapGraphWrapper::Node& n) const {
735      //        return dist[n]; }
736      //       bool get(const typename MapGraphWrapper::Edge& e) const {
737      //        return (dist.get(g->tail(e))<dist.get(g->head(e))); }
738      bool operator[](const typename MapGraphWrapper::Edge& e) const {
739        return (dist[g->tail(e)]<dist[g->head(e)]);
740      }
741    };
742
743  };
744
745
746  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
747  void MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::preflowPhase1(FlowEnum fe)
748  {
749
750    int heur0=(int)(H0*n);  //time while running 'bound decrease'
751    int heur1=(int)(H1*n);  //time while running 'highest label'
752    int heur=heur1;         //starting time interval (#of relabels)
753    int numrelabel=0;
754
755    bool what_heur=1;
756    //It is 0 in case 'bound decrease' and 1 in case 'highest label'
757
758    bool end=false;
759    //Needed for 'bound decrease', true means no active nodes are above bound
760    //b.
761
762    int k=n-2;  //bound on the highest level under n containing a node
763    int b=k;    //bound on the highest level under n of an active node
764
765    VecFirst first(n, INVALID);
766    NNMap next(*g, INVALID); //maybe INVALID is not needed
767    //    VecStack active(n);
768
769    NNMap left(*g, INVALID);
770    NNMap right(*g, INVALID);
771    VecNode level_list(n,INVALID);
772    //List of the nodes in level i<n, set to n.
773
774    NodeIt v;
775    for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
776    //setting each node to level n
777
778    if ( fe == NO_FLOW ) {
779      EdgeIt e;
780      for(g->first(e); g->valid(e); g->next(e)) flow->set(e,0);
781    }
782
783    switch (fe) { //computing the excess
784    case PRE_FLOW:
785      {
786        NodeIt v;
787        for(g->first(v); g->valid(v); g->next(v)) {
788          Num exc=0;
789
790          InEdgeIt e;
791          for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
792          OutEdgeIt f;
793          for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
794
795          excess.set(v,exc);
796
797          //putting the active nodes into the stack
798          int lev=level[v];
799          if ( exc > 0 && lev < n && v != t )
800            {
801              next.set(v,first[lev]);
802              first[lev]=v;
803            }
804          //      active[lev].push(v);
805        }
806        break;
807      }
808    case GEN_FLOW:
809      {
810        NodeIt v;
811        for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
812
813        Num exc=0;
814        InEdgeIt e;
815        for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
816        OutEdgeIt f;
817        for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
818        excess.set(t,exc);
819        break;
820      }
821    case ZERO_FLOW:
822    case NO_FLOW:
823      {
824        NodeIt v;
825        for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
826        break;
827      }
828    }
829
830    preflowPreproc(fe, next, first,/*active*/ level_list, left, right);
831    //End of preprocessing
832
833
834    //Push/relabel on the highest level active nodes.
835    while ( true ) {
836      if ( b == 0 ) {
837        if ( !what_heur && !end && k > 0 ) {
838          b=k;
839          end=true;
840        } else break;
841      }
842
843      if ( !g->valid(first[b])/*active[b].empty()*/ ) --b;
844      else {
845        end=false;
846        Node w=first[b];
847        first[b]=next[w];
848        /*      Node w=active[b].top();
849                active[b].pop();*/
850        int newlevel=push(w,/*active*/next, first);
851        if ( excess[w] > 0 ) relabel(w, newlevel, /*active*/next, first, level_list,
852                                     left, right, b, k, what_heur);
853
854        ++numrelabel;
855        if ( numrelabel >= heur ) {
856          numrelabel=0;
857          if ( what_heur ) {
858            what_heur=0;
859            heur=heur0;
860            end=false;
861          } else {
862            what_heur=1;
863            heur=heur1;
864            b=k;
865          }
866        }
867      }
868    }
869
870    status=AFTER_PRE_FLOW_PHASE_1;
871  }
872
873
874
875  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
876  void MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::preflowPhase2()
877  {
878
879    int k=n-2;  //bound on the highest level under n containing a node
880    int b=k;    //bound on the highest level under n of an active node
881
882   
883    VecFirst first(n, INVALID);
884    NNMap next(*g, INVALID); //maybe INVALID is not needed
885    //    VecStack active(n);
886    level.set(s,0);
887    std::queue<Node> bfs_queue;
888    bfs_queue.push(s);
889
890    while (!bfs_queue.empty()) {
891
892      Node v=bfs_queue.front();
893      bfs_queue.pop();
894      int l=level[v]+1;
895
896      InEdgeIt e;
897      for(g->first(e,v); g->valid(e); g->next(e)) {
898        if ( (*capacity)[e] <= (*flow)[e] ) continue;
899        Node u=g->tail(e);
900        if ( level[u] >= n ) {
901          bfs_queue.push(u);
902          level.set(u, l);
903          if ( excess[u] > 0 ) {
904            next.set(u,first[l]);
905            first[l]=u;
906            //active[l].push(u);
907          }
908        }
909      }
910
911      OutEdgeIt f;
912      for(g->first(f,v); g->valid(f); g->next(f)) {
913        if ( 0 >= (*flow)[f] ) continue;
914        Node u=g->head(f);
915        if ( level[u] >= n ) {
916          bfs_queue.push(u);
917          level.set(u, l);
918          if ( excess[u] > 0 ) {
919            next.set(u,first[l]);
920            first[l]=u;
921            //active[l].push(u);
922          }
923        }
924      }
925    }
926    b=n-2;
927
928    while ( true ) {
929
930      if ( b == 0 ) break;
931
932      if ( !g->valid(first[b])/*active[b].empty()*/ ) --b;
933      else {
934
935        Node w=first[b];
936        first[b]=next[w];
937        /*      Node w=active[b].top();
938                active[b].pop();*/
939        int newlevel=push(w,next, first/*active*/);
940
941        //relabel
942        if ( excess[w] > 0 ) {
943          level.set(w,++newlevel);
944          next.set(w,first[newlevel]);
945          first[newlevel]=w;
946          //active[newlevel].push(w);
947          b=newlevel;
948        }
949      }  // if stack[b] is nonempty
950    } // while(true)
951
952    status=AFTER_PRE_FLOW_PHASE_2;
953  }
954
955
956
957  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
958  bool MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath()
959  {
960    ResGW res_graph(*g, *capacity, *flow);
961    bool _augment=false;
962
963    //ReachedMap level(res_graph);
964    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
965    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
966    bfs.pushAndSetReached(s);
967
968    typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
969    pred.set(s, INVALID);
970
971    typename ResGW::template NodeMap<Num> free(res_graph);
972
973    //searching for augmenting path
974    while ( !bfs.finished() ) {
975      ResGWOutEdgeIt e=bfs;
976      if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
977        Node v=res_graph.tail(e);
978        Node w=res_graph.head(e);
979        pred.set(w, e);
980        if (res_graph.valid(pred[v])) {
981          free.set(w, std::min(free[v], res_graph.resCap(e)));
982        } else {
983          free.set(w, res_graph.resCap(e));
984        }
985        if (res_graph.head(e)==t) { _augment=true; break; }
986      }
987
988      ++bfs;
989    } //end of searching augmenting path
990
991    if (_augment) {
992      Node n=t;
993      Num augment_value=free[t];
994      while (res_graph.valid(pred[n])) {
995        ResGWEdge e=pred[n];
996        res_graph.augment(e, augment_value);
997        n=res_graph.tail(e);
998      }
999    }
1000
1001    status=AFTER_AUGMENTING;
1002    return _augment;
1003  }
1004
1005
1006  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
1007  bool MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath2()
1008  {
1009    ResGW res_graph(*g, *capacity, *flow);
1010    bool _augment=false;
1011
1012    if (status!=AFTER_FAST_AUGMENTING) {
1013      FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
1014      number_of_augmentations=1;
1015    } else {
1016      ++number_of_augmentations;
1017    }
1018    TrickyReachedMap<ReachedMap>
1019      tricky_reached_map(level, number_of_augmentations);
1020    //ReachedMap level(res_graph);
1021//    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
1022    BfsIterator<ResGW, TrickyReachedMap<ReachedMap> >
1023      bfs(res_graph, tricky_reached_map);
1024    bfs.pushAndSetReached(s);
1025
1026    typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
1027    pred.set(s, INVALID);
1028
1029    typename ResGW::template NodeMap<Num> free(res_graph);
1030
1031    //searching for augmenting path
1032    while ( !bfs.finished() ) {
1033      ResGWOutEdgeIt e=bfs;
1034      if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
1035        Node v=res_graph.tail(e);
1036        Node w=res_graph.head(e);
1037        pred.set(w, e);
1038        if (res_graph.valid(pred[v])) {
1039          free.set(w, std::min(free[v], res_graph.resCap(e)));
1040        } else {
1041          free.set(w, res_graph.resCap(e));
1042        }
1043        if (res_graph.head(e)==t) { _augment=true; break; }
1044      }
1045
1046      ++bfs;
1047    } //end of searching augmenting path
1048
1049    if (_augment) {
1050      Node n=t;
1051      Num augment_value=free[t];
1052      while (res_graph.valid(pred[n])) {
1053        ResGWEdge e=pred[n];
1054        res_graph.augment(e, augment_value);
1055        n=res_graph.tail(e);
1056      }
1057    }
1058
1059    status=AFTER_FAST_AUGMENTING;
1060    return _augment;
1061  }
1062
1063
1064  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
1065  template<typename MutableGraph>
1066  bool MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow()
1067  {
1068    typedef MutableGraph MG;
1069    bool _augment=false;
1070
1071    ResGW res_graph(*g, *capacity, *flow);
1072
1073    //bfs for distances on the residual graph
1074    //ReachedMap level(res_graph);
1075    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
1076    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
1077    bfs.pushAndSetReached(s);
1078    typename ResGW::template NodeMap<int>
1079      dist(res_graph); //filled up with 0's
1080
1081    //F will contain the physical copy of the residual graph
1082    //with the set of edges which are on shortest paths
1083    MG F;
1084    typename ResGW::template NodeMap<typename MG::Node>
1085      res_graph_to_F(res_graph);
1086    {
1087      typename ResGW::NodeIt n;
1088      for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
1089        res_graph_to_F.set(n, F.addNode());
1090      }
1091    }
1092
1093    typename MG::Node sF=res_graph_to_F[s];
1094    typename MG::Node tF=res_graph_to_F[t];
1095    typename MG::template EdgeMap<ResGWEdge> original_edge(F);
1096    typename MG::template EdgeMap<Num> residual_capacity(F);
1097
1098    while ( !bfs.finished() ) {
1099      ResGWOutEdgeIt e=bfs;
1100      if (res_graph.valid(e)) {
1101        if (bfs.isBNodeNewlyReached()) {
1102          dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
1103          typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)],
1104                                        res_graph_to_F[res_graph.head(e)]);
1105          original_edge.update();
1106          original_edge.set(f, e);
1107          residual_capacity.update();
1108          residual_capacity.set(f, res_graph.resCap(e));
1109        } else {
1110          if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
1111            typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)],
1112                                          res_graph_to_F[res_graph.head(e)]);
1113            original_edge.update();
1114            original_edge.set(f, e);
1115            residual_capacity.update();
1116            residual_capacity.set(f, res_graph.resCap(e));
1117          }
1118        }
1119      }
1120      ++bfs;
1121    } //computing distances from s in the residual graph
1122
1123    bool __augment=true;
1124
1125    while (__augment) {
1126      __augment=false;
1127      //computing blocking flow with dfs
1128      DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F);
1129      typename MG::template NodeMap<typename MG::Edge> pred(F);
1130      pred.set(sF, INVALID);
1131      //invalid iterators for sources
1132
1133      typename MG::template NodeMap<Num> free(F);
1134
1135      dfs.pushAndSetReached(sF);
1136      while (!dfs.finished()) {
1137        ++dfs;
1138        if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
1139          if (dfs.isBNodeNewlyReached()) {
1140            typename MG::Node v=F.aNode(dfs);
1141            typename MG::Node w=F.bNode(dfs);
1142            pred.set(w, dfs);
1143            if (F.valid(pred[v])) {
1144              free.set(w, std::min(free[v], residual_capacity[dfs]));
1145            } else {
1146              free.set(w, residual_capacity[dfs]);
1147            }
1148            if (w==tF) {
1149              __augment=true;
1150              _augment=true;
1151              break;
1152            }
1153
1154          } else {
1155            F.erase(/*typename MG::OutEdgeIt*/(dfs));
1156          }
1157        }
1158      }
1159
1160      if (__augment) {
1161        typename MG::Node n=tF;
1162        Num augment_value=free[tF];
1163        while (F.valid(pred[n])) {
1164          typename MG::Edge e=pred[n];
1165          res_graph.augment(original_edge[e], augment_value);
1166          n=F.tail(e);
1167          if (residual_capacity[e]==augment_value)
1168            F.erase(e);
1169          else
1170            residual_capacity.set(e, residual_capacity[e]-augment_value);
1171        }
1172      }
1173
1174    }
1175
1176    status=AFTER_AUGMENTING;
1177    return _augment;
1178  }
1179
1180
1181
1182
1183  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
1184  bool MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2()
1185  {
1186    bool _augment=false;
1187
1188    ResGW res_graph(*g, *capacity, *flow);
1189
1190    //ReachedMap level(res_graph);
1191    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
1192    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
1193
1194    bfs.pushAndSetReached(s);
1195    DistanceMap<ResGW> dist(res_graph);
1196    while ( !bfs.finished() ) {
1197      ResGWOutEdgeIt e=bfs;
1198      if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
1199        dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
1200      }
1201      ++bfs;
1202    } //computing distances from s in the residual graph
1203
1204      //Subgraph containing the edges on some shortest paths
1205    ConstMap<typename ResGW::Node, bool> true_map(true);
1206    typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>,
1207      DistanceMap<ResGW> > FilterResGW;
1208    FilterResGW filter_res_graph(res_graph, true_map, dist);
1209
1210    //Subgraph, which is able to delete edges which are already
1211    //met by the dfs
1212    typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt>
1213      first_out_edges(filter_res_graph);
1214    typename FilterResGW::NodeIt v;
1215    for(filter_res_graph.first(v); filter_res_graph.valid(v);
1216        filter_res_graph.next(v))
1217      {
1218        typename FilterResGW::OutEdgeIt e;
1219        filter_res_graph.first(e, v);
1220        first_out_edges.set(v, e);
1221      }
1222    typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
1223      template NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW;
1224    ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);
1225
1226    bool __augment=true;
1227
1228    while (__augment) {
1229
1230      __augment=false;
1231      //computing blocking flow with dfs
1232      DfsIterator< ErasingResGW,
1233        typename ErasingResGW::template NodeMap<bool> >
1234        dfs(erasing_res_graph);
1235      typename ErasingResGW::
1236        template NodeMap<typename ErasingResGW::OutEdgeIt>
1237        pred(erasing_res_graph);
1238      pred.set(s, INVALID);
1239      //invalid iterators for sources
1240
1241      typename ErasingResGW::template NodeMap<Num>
1242        free1(erasing_res_graph);
1243
1244      dfs.pushAndSetReached
1245        ///\bug hugo 0.2
1246        (typename ErasingResGW::Node
1247         (typename FilterResGW::Node
1248          (typename ResGW::Node(s)
1249           )
1250          )
1251         );
1252      while (!dfs.finished()) {
1253        ++dfs;
1254        if (erasing_res_graph.valid(typename ErasingResGW::OutEdgeIt(dfs)))
1255          {
1256            if (dfs.isBNodeNewlyReached()) {
1257
1258              typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs);
1259              typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs);
1260
1261              pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs));
1262              if (erasing_res_graph.valid(pred[v])) {
1263                free1.set
1264                  (w, std::min(free1[v], res_graph.resCap
1265                               (typename ErasingResGW::OutEdgeIt(dfs))));
1266              } else {
1267                free1.set
1268                  (w, res_graph.resCap
1269                   (typename ErasingResGW::OutEdgeIt(dfs)));
1270              }
1271
1272              if (w==t) {
1273                __augment=true;
1274                _augment=true;
1275                break;
1276              }
1277            } else {
1278              erasing_res_graph.erase(dfs);
1279            }
1280          }
1281      }
1282
1283      if (__augment) {
1284        typename ErasingResGW::Node
1285          n=typename FilterResGW::Node(typename ResGW::Node(t));
1286        //        typename ResGW::NodeMap<Num> a(res_graph);
1287        //        typename ResGW::Node b;
1288        //        Num j=a[b];
1289        //        typename FilterResGW::NodeMap<Num> a1(filter_res_graph);
1290        //        typename FilterResGW::Node b1;
1291        //        Num j1=a1[b1];
1292        //        typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);
1293        //        typename ErasingResGW::Node b2;
1294        //        Num j2=a2[b2];
1295        Num augment_value=free1[n];
1296        while (erasing_res_graph.valid(pred[n])) {
1297          typename ErasingResGW::OutEdgeIt e=pred[n];
1298          res_graph.augment(e, augment_value);
1299          n=erasing_res_graph.tail(e);
1300          if (res_graph.resCap(e)==0)
1301            erasing_res_graph.erase(e);
1302        }
1303      }
1304
1305    } //while (__augment)
1306
1307    status=AFTER_AUGMENTING;
1308    return _augment;
1309  }
1310
1311
1312} //namespace hugo
1313
1314#endif //HUGO_MAX_FLOW_H
1315
1316
1317
1318
Note: See TracBrowser for help on using the repository browser.