COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/preflow.h @ 1449:ac7e995e47e2

Last change on this file since 1449:ac7e995e47e2 was 1435:8e85e6bbefdf, checked in by Akos Ladanyi, 19 years ago

trunk/src/* move to trunk/

File size: 23.2 KB
RevLine 
[906]1/* -*- C++ -*-
[1435]2 * lemon/preflow.h - Part of LEMON, a generic C++ optimization library
[906]3 *
[1164]4 * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
[1359]5 * (Egervary Research Group on Combinatorial Optimization, EGRES).
[906]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
[921]17#ifndef LEMON_PREFLOW_H
18#define LEMON_PREFLOW_H
[836]19
20#include <vector>
21#include <queue>
22
[921]23#include <lemon/invalid.h>
24#include <lemon/maps.h>
[977]25#include <lemon/graph_utils.h>
[836]26
27/// \file
28/// \ingroup flowalgs
[874]29/// Implementation of the preflow algorithm.
[836]30
[921]31namespace lemon {
[836]32
33  /// \addtogroup flowalgs
34  /// @{                                                   
35
[851]36  ///%Preflow algorithms class.
[836]37
38  ///This class provides an implementation of the \e preflow \e
39  ///algorithm producing a flow of maximum value in a directed
[1222]40  ///graph. The preflow algorithms are the fastest known max flow algorithms
[851]41  ///up to now. The \e source node, the \e target node, the \e
[836]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
[1285]45  ///functions \ref source, \ref target, \ref capacityMap and \ref
46  ///flowMap.
[836]47  ///
[921]48  ///After running \ref lemon::Preflow::phase1() "phase1()"
49  ///or \ref lemon::Preflow::run() "run()", the maximal flow
[836]50  ///value can be obtained by calling \ref flowValue(). The minimum
[851]51  ///value cut can be written into a <tt>bool</tt> node map by
52  ///calling \ref minCut(). (\ref minMinCut() and \ref maxMinCut() writes
[836]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.
[1222]58  ///\param CapacityMap The capacity map type.
[836]59  ///\param FlowMap The flow map type.
60  ///
61  ///\author Jacint Szabo
[1227]62  ///\todo Second template parameter is superfluous
[836]63  template <typename Graph, typename Num,
[1222]64            typename CapacityMap=typename Graph::template EdgeMap<Num>,
[836]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
[1222]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
[836]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
[1222]95    ///Indicates the property of the starting flow map.
96    ///The meanings are as follows:
[836]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.
[911]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.
[836]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
[1222]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.
[836]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.
[1285]137    ///\param _gr The directed graph the algorithm runs on.
[836]138    ///\param _s The source node.
139    ///\param _t The target node.
[1222]140    ///\param _cap The capacity of the edges.
141    ///\param _f The flow of the edges.
[836]142    ///Except the graph, all of these parameters can be reset by
[1285]143    ///calling \ref source, \ref target, \ref capacityMap and \ref
144    ///flowMap, resp.
[1222]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),
[836]149        flow_prop(NO_FLOW), status(AFTER_NOTHING) { }
150
151
152                                                                             
153    ///Runs the preflow algorithm. 
154
[851]155    ///Runs the preflow algorithm.
156    ///
[836]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
[920]179    ///The preflow algorithm consists of two phases, this method runs
180    ///the first phase. After the first phase the maximum flow value
[1285]181    ///and a minimum value cut can already be computed, although a
[920]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.
[836]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
[920]201    ///The preflow algorithm consists of two phases, this method runs
202    ///the first phase. After the first phase the maximum flow value
[1285]203    ///and a minimum value cut can already be computed, although a
[920]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.
[911]207    ///\warning \ref minCut(), \ref minMinCut() and \ref maxMinCut() do not
208    ///give minimum value cuts unless calling \ref phase2().
[836]209    void phase1()
210    {
[1222]211      int heur0=(int)(H0*_node_num);  //time while running 'bound decrease'
212      int heur1=(int)(H1*_node_num);  //time while running 'highest label'
[836]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
[1222]223      int k=_node_num-2;  //bound on the highest level under n containing a node
[836]224      int b=k;    //bound on the highest level under n of an active node
225
[1222]226      VecNode first(_node_num, INVALID);
227      NNMap next(*_g, INVALID);
[836]228
[1222]229      NNMap left(*_g, INVALID);
230      NNMap right(*_g, INVALID);
231      VecNode level_list(_node_num,INVALID);
[836]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
[1222]278    //   runs heuristic 'bound decrease' for H0*n relabels,
279    //        starts with 'highest label'
[836]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
[920]286    ///the second phase. After calling \ref phase1 and then \ref
287    ///phase2, \ref flow contains a maximum flow, \ref flowValue
288    ///returns the value of a maximum flow, \ref minCut returns a
289    ///minimum cut, while the methods \ref minMinCut and \ref
290    ///maxMinCut return the inclusionwise minimum and maximum cuts of
291    ///minimum value, resp.  \pre \ref phase1 must be called before.
[836]292    void phase2()
293    {
294
[1222]295      int k=_node_num-2;  //bound on the highest level under n containing a node
[836]296      int b=k;    //bound on the highest level under n of an active node
297
298   
[1222]299      VecNode first(_node_num, INVALID);
300      NNMap next(*_g, INVALID);
301      level.set(_source,0);
[836]302      std::queue<Node> bfs_queue;
[1222]303      bfs_queue.push(_source);
[836]304
305      while ( !bfs_queue.empty() ) {
306
307        Node v=bfs_queue.front();
308        bfs_queue.pop();
309        int l=level[v]+1;
310
[1222]311        for(InEdgeIt e(*_g,v); e!=INVALID; ++e) {
312          if ( (*_capacity)[e] <= (*_flow)[e] ) continue;
313          Node u=_g->source(e);
314          if ( level[u] >= _node_num ) {
[836]315            bfs_queue.push(u);
316            level.set(u, l);
317            if ( excess[u] > 0 ) {
318              next.set(u,first[l]);
319              first[l]=u;
320            }
321          }
322        }
323
[1222]324        for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) {
325          if ( 0 >= (*_flow)[e] ) continue;
326          Node u=_g->target(e);
327          if ( level[u] >= _node_num ) {
[836]328            bfs_queue.push(u);
329            level.set(u, l);
330            if ( excess[u] > 0 ) {
331              next.set(u,first[l]);
332              first[l]=u;
333            }
334          }
335        }
336      }
[1222]337      b=_node_num-2;
[836]338
339      while ( true ) {
340
341        if ( b == 0 ) break;
342        if ( first[b]==INVALID ) --b;
343        else {
344          Node w=first[b];
345          first[b]=next[w];
346          int newlevel=push(w,next, first);
347         
348          //relabel
349          if ( excess[w] > 0 ) {
350            level.set(w,++newlevel);
351            next.set(w,first[newlevel]);
352            first[newlevel]=w;
353            b=newlevel;
354          }
355        }
356      } // while(true)
357      flow_prop=GEN_FLOW;
358      status=AFTER_PREFLOW_PHASE_2;
359    }
360
361    /// Returns the value of the maximum flow.
362
363    /// Returns the value of the maximum flow by returning the excess
[911]364    /// of the target node \c t. This value equals to the value of
[836]365    /// the maximum flow already after running \ref phase1.
366    Num flowValue() const {
[1222]367      return excess[_target];
[836]368    }
369
370
371    ///Returns a minimum value cut.
372
373    ///Sets \c M to the characteristic vector of a minimum value
374    ///cut. This method can be called both after running \ref
375    ///phase1 and \ref phase2. It is much faster after
[849]376    ///\ref phase1.  \pre M should be a bool-valued node-map. \pre
[911]377    ///If \ref minCut() is called after \ref phase2() then M should
[836]378    ///be initialized to false.
379    template<typename _CutMap>
380    void minCut(_CutMap& M) const {
381      switch ( status ) {
382        case AFTER_PREFLOW_PHASE_1:
[1222]383        for(NodeIt v(*_g); v!=INVALID; ++v) {
384          if (level[v] < _node_num) {
[836]385            M.set(v, false);
386          } else {
387            M.set(v, true);
388          }
389        }
390        break;
391        case AFTER_PREFLOW_PHASE_2:
392        minMinCut(M);
393        break;
394        case AFTER_NOTHING:
395        break;
396      }
397    }
398
399    ///Returns the inclusionwise minimum of the minimum value cuts.
400
401    ///Sets \c M to the characteristic vector of the minimum value cut
402    ///which is inclusionwise minimum. It is computed by processing a
403    ///bfs from the source node \c s in the residual graph.  \pre M
404    ///should be a node map of bools initialized to false.  \pre \ref
405    ///phase2 should already be run.
406    template<typename _CutMap>
407    void minMinCut(_CutMap& M) const {
408
409      std::queue<Node> queue;
[1222]410      M.set(_source,true);
[1227]411      queue.push(_source);
[836]412     
413      while (!queue.empty()) {
414        Node w=queue.front();
415        queue.pop();
416       
[1222]417        for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
418          Node v=_g->target(e);
419          if (!M[v] && (*_flow)[e] < (*_capacity)[e] ) {
[836]420            queue.push(v);
421            M.set(v, true);
422          }
423        }
424       
[1222]425        for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
426          Node v=_g->source(e);
427          if (!M[v] && (*_flow)[e] > 0 ) {
[836]428            queue.push(v);
429            M.set(v, true);
430          }
431        }
432      }
433    }
434   
435    ///Returns the inclusionwise maximum of the minimum value cuts.
436
437    ///Sets \c M to the characteristic vector of the minimum value cut
438    ///which is inclusionwise maximum. It is computed by processing a
439    ///backward bfs from the target node \c t in the residual graph.
[911]440    ///\pre \ref phase2() or run() should already be run.
[836]441    template<typename _CutMap>
442    void maxMinCut(_CutMap& M) const {
443
[1222]444      for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true);
[836]445
446      std::queue<Node> queue;
447
[1222]448      M.set(_target,false);
449      queue.push(_target);
[836]450
451      while (!queue.empty()) {
452        Node w=queue.front();
453        queue.pop();
454
[1222]455        for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
456          Node v=_g->source(e);
457          if (M[v] && (*_flow)[e] < (*_capacity)[e] ) {
[836]458            queue.push(v);
459            M.set(v, false);
460          }
461        }
462
[1222]463        for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
464          Node v=_g->target(e);
465          if (M[v] && (*_flow)[e] > 0 ) {
[836]466            queue.push(v);
467            M.set(v, false);
468          }
469        }
470      }
471    }
472
473    ///Sets the source node to \c _s.
474
475    ///Sets the source node to \c _s.
476    ///
[1222]477    void source(Node _s) {
478      _source=_s;
[836]479      if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
480      status=AFTER_NOTHING;
481    }
482
[1222]483    ///Returns the source node.
484
485    ///Returns the source node.
486    ///
487    Node source() const {
488      return _source;
489    }
490
[836]491    ///Sets the target node to \c _t.
492
493    ///Sets the target node to \c _t.
494    ///
[1222]495    void target(Node _t) {
496      _target=_t;
[836]497      if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
498      status=AFTER_NOTHING;
499    }
500
[1222]501    ///Returns the target node.
502
503    ///Returns the target node.
504    ///
505    Node target() const {
506      return _target;
507    }
508
[836]509    /// Sets the edge map of the capacities to _cap.
510
511    /// Sets the edge map of the capacities to _cap.
512    ///
[1222]513    void capacityMap(const CapacityMap& _cap) {
514      _capacity=&_cap;
[836]515      status=AFTER_NOTHING;
516    }
[1285]517    /// Returns a reference to capacity map.
[1222]518
[1285]519    /// Returns a reference to capacity map.
[1222]520    ///
521    const CapacityMap &capacityMap() const {
522      return *_capacity;
523    }
[836]524
525    /// Sets the edge map of the flows to _flow.
526
527    /// Sets the edge map of the flows to _flow.
528    ///
[1222]529    void flowMap(FlowMap& _f) {
530      _flow=&_f;
[836]531      flow_prop=NO_FLOW;
532      status=AFTER_NOTHING;
533    }
[1222]534     
[1285]535    /// Returns a reference to flow map.
[836]536
[1285]537    /// Returns a reference to flow map.
[1222]538    ///
539    const FlowMap &flowMap() const {
540      return *_flow;
541    }
[836]542
543  private:
544
545    int push(Node w, NNMap& next, VecNode& first) {
546
547      int lev=level[w];
548      Num exc=excess[w];
[1222]549      int newlevel=_node_num;       //bound on the next level of w
[836]550
[1222]551      for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
552        if ( (*_flow)[e] >= (*_capacity)[e] ) continue;
553        Node v=_g->target(e);
[836]554
555        if( lev > level[v] ) { //Push is allowed now
556         
[1222]557          if ( excess[v]<=0 && v!=_target && v!=_source ) {
[836]558            next.set(v,first[level[v]]);
559            first[level[v]]=v;
560          }
561
[1222]562          Num cap=(*_capacity)[e];
563          Num flo=(*_flow)[e];
[836]564          Num remcap=cap-flo;
565         
566          if ( remcap >= exc ) { //A nonsaturating push.
567           
[1222]568            _flow->set(e, flo+exc);
[836]569            excess.set(v, excess[v]+exc);
570            exc=0;
571            break;
572
573          } else { //A saturating push.
[1222]574            _flow->set(e, cap);
[836]575            excess.set(v, excess[v]+remcap);
576            exc-=remcap;
577          }
578        } else if ( newlevel > level[v] ) newlevel = level[v];
579      } //for out edges wv
580
581      if ( exc > 0 ) {
[1222]582        for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
[836]583         
[1222]584          if( (*_flow)[e] <= 0 ) continue;
585          Node v=_g->source(e);
[836]586
587          if( lev > level[v] ) { //Push is allowed now
588
[1222]589            if ( excess[v]<=0 && v!=_target && v!=_source ) {
[836]590              next.set(v,first[level[v]]);
591              first[level[v]]=v;
592            }
593
[1222]594            Num flo=(*_flow)[e];
[836]595
596            if ( flo >= exc ) { //A nonsaturating push.
597
[1222]598              _flow->set(e, flo-exc);
[836]599              excess.set(v, excess[v]+exc);
600              exc=0;
601              break;
602            } else {  //A saturating push.
603
604              excess.set(v, excess[v]+flo);
605              exc-=flo;
[1222]606              _flow->set(e,0);
[836]607            }
608          } else if ( newlevel > level[v] ) newlevel = level[v];
609        } //for in edges vw
610
611      } // if w still has excess after the out edge for cycle
612
613      excess.set(w, exc);
614     
615      return newlevel;
616    }
617   
618   
619   
620    void preflowPreproc(VecNode& first, NNMap& next,
621                        VecNode& level_list, NNMap& left, NNMap& right)
622    {
[1222]623      for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num);
[836]624      std::queue<Node> bfs_queue;
625     
626      if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
627        //Reverse_bfs from t in the residual graph,
628        //to find the starting level.
[1222]629        level.set(_target,0);
630        bfs_queue.push(_target);
[836]631       
632        while ( !bfs_queue.empty() ) {
633         
634          Node v=bfs_queue.front();
635          bfs_queue.pop();
636          int l=level[v]+1;
637         
[1222]638          for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
639            if ( (*_capacity)[e] <= (*_flow)[e] ) continue;
640            Node w=_g->source(e);
641            if ( level[w] == _node_num && w != _source ) {
[836]642              bfs_queue.push(w);
643              Node z=level_list[l];
644              if ( z!=INVALID ) left.set(z,w);
645              right.set(w,z);
646              level_list[l]=w;
647              level.set(w, l);
648            }
649          }
650         
[1222]651          for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
652            if ( 0 >= (*_flow)[e] ) continue;
653            Node w=_g->target(e);
654            if ( level[w] == _node_num && w != _source ) {
[836]655              bfs_queue.push(w);
656              Node z=level_list[l];
657              if ( z!=INVALID ) left.set(z,w);
658              right.set(w,z);
659              level_list[l]=w;
660              level.set(w, l);
661            }
662          }
663        } //while
664      } //if
665
666
667      switch (flow_prop) {
668        case NO_FLOW: 
[1222]669        for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0);
[836]670        case ZERO_FLOW:
[1222]671        for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
[836]672       
673        //Reverse_bfs from t, to find the starting level.
[1222]674        level.set(_target,0);
675        bfs_queue.push(_target);
[836]676       
677        while ( !bfs_queue.empty() ) {
678         
679          Node v=bfs_queue.front();
680          bfs_queue.pop();
681          int l=level[v]+1;
682         
[1222]683          for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
684            Node w=_g->source(e);
685            if ( level[w] == _node_num && w != _source ) {
[836]686              bfs_queue.push(w);
687              Node z=level_list[l];
688              if ( z!=INVALID ) left.set(z,w);
689              right.set(w,z);
690              level_list[l]=w;
691              level.set(w, l);
692            }
693          }
694        }
695       
696        //the starting flow
[1222]697        for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
698          Num c=(*_capacity)[e];
[836]699          if ( c <= 0 ) continue;
[1222]700          Node w=_g->target(e);
701          if ( level[w] < _node_num ) {
702            if ( excess[w] <= 0 && w!=_target ) { //putting into the stack
[836]703              next.set(w,first[level[w]]);
704              first[level[w]]=w;
705            }
[1222]706            _flow->set(e, c);
[836]707            excess.set(w, excess[w]+c);
708          }
709        }
710        break;
711
712        case GEN_FLOW:
[1222]713        for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
[836]714        {
715          Num exc=0;
[1222]716          for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e];
717          for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e];
718          excess.set(_target,exc);
[836]719        }
720
721        //the starting flow
[1222]722        for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e)  {
723          Num rem=(*_capacity)[e]-(*_flow)[e];
[836]724          if ( rem <= 0 ) continue;
[1222]725          Node w=_g->target(e);
726          if ( level[w] < _node_num ) {
727            if ( excess[w] <= 0 && w!=_target ) { //putting into the stack
[836]728              next.set(w,first[level[w]]);
729              first[level[w]]=w;
730            }   
[1222]731            _flow->set(e, (*_capacity)[e]);
[836]732            excess.set(w, excess[w]+rem);
733          }
734        }
735       
[1222]736        for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) {
737          if ( (*_flow)[e] <= 0 ) continue;
738          Node w=_g->source(e);
739          if ( level[w] < _node_num ) {
740            if ( excess[w] <= 0 && w!=_target ) {
[836]741              next.set(w,first[level[w]]);
742              first[level[w]]=w;
743            } 
[1222]744            excess.set(w, excess[w]+(*_flow)[e]);
745            _flow->set(e, 0);
[836]746          }
747        }
748        break;
749
750        case PRE_FLOW: 
751        //the starting flow
[1222]752        for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
753          Num rem=(*_capacity)[e]-(*_flow)[e];
[836]754          if ( rem <= 0 ) continue;
[1222]755          Node w=_g->target(e);
756          if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]);
[836]757        }
758       
[1222]759        for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
760          if ( (*_flow)[e] <= 0 ) continue;
761          Node w=_g->source(e);
762          if ( level[w] < _node_num ) _flow->set(e, 0);
[836]763        }
764       
765        //computing the excess
[1222]766        for(NodeIt w(*_g); w!=INVALID; ++w) {
[836]767          Num exc=0;
[1222]768          for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e];
769          for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e];
[836]770          excess.set(w,exc);
771         
772          //putting the active nodes into the stack
773          int lev=level[w];
[1222]774            if ( exc > 0 && lev < _node_num && Node(w) != _target ) {
[836]775              next.set(w,first[lev]);
776              first[lev]=w;
777            }
778        }
779        break;
780      } //switch
781    } //preflowPreproc
782
783
784    void relabel(Node w, int newlevel, VecNode& first, NNMap& next,
785                 VecNode& level_list, NNMap& left,
786                 NNMap& right, int& b, int& k, bool what_heur )
787    {
788
789      int lev=level[w];
790
791      Node right_n=right[w];
792      Node left_n=left[w];
793
794      //unlacing starts
795      if ( right_n!=INVALID ) {
796        if ( left_n!=INVALID ) {
797          right.set(left_n, right_n);
798          left.set(right_n, left_n);
799        } else {
800          level_list[lev]=right_n;
801          left.set(right_n, INVALID);
802        }
803      } else {
804        if ( left_n!=INVALID ) {
805          right.set(left_n, INVALID);
806        } else {
807          level_list[lev]=INVALID;
808        }
809      }
810      //unlacing ends
811
812      if ( level_list[lev]==INVALID ) {
813
814        //gapping starts
815        for (int i=lev; i!=k ; ) {
816          Node v=level_list[++i];
817          while ( v!=INVALID ) {
[1222]818            level.set(v,_node_num);
[836]819            v=right[v];
820          }
821          level_list[i]=INVALID;
822          if ( !what_heur ) first[i]=INVALID;
823        }
824
[1222]825        level.set(w,_node_num);
[836]826        b=lev-1;
827        k=b;
828        //gapping ends
829
830      } else {
831
[1222]832        if ( newlevel == _node_num ) level.set(w,_node_num);
[836]833        else {
834          level.set(w,++newlevel);
835          next.set(w,first[newlevel]);
836          first[newlevel]=w;
837          if ( what_heur ) b=newlevel;
838          if ( k < newlevel ) ++k;      //now k=newlevel
839          Node z=level_list[newlevel];
840          if ( z!=INVALID ) left.set(z,w);
841          right.set(w,z);
842          left.set(w,INVALID);
843          level_list[newlevel]=w;
844        }
845      }
846    } //relabel
847
848  };
[1227]849
850  ///Function type interface for Preflow algorithm.
851
852  /// \ingroup flowalgs
853  ///Function type interface for Preflow algorithm.
854  ///\sa Preflow
855  template<class GR, class CM, class FM>
856  Preflow<GR,typename CM::Value,CM,FM> preflow(const GR &g,
857                            typename GR::Node source,
858                            typename GR::Node target,
859                            const CM &cap,
860                            FM &flow
861                            )
862  {
863    return Preflow<GR,typename CM::Value,CM,FM>(g,source,target,cap,flow);
864  }
865
[921]866} //namespace lemon
[836]867
[921]868#endif //LEMON_PREFLOW_H
Note: See TracBrowser for help on using the repository browser.