COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/max_save.h @ 854:baf0b6e40211

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