COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/preflow.h @ 1742:2637b9420d0a

Last change on this file since 1742:2637b9420d0a was 1742:2637b9420d0a, checked in by Balazs Dezso, 14 years ago

Show description of the file

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