COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/lemon/preflow.h @ 1284:b941d044f87b

Last change on this file since 1284:b941d044f87b was 1227:01f668e3e168, checked in by Alpar Juttner, 15 years ago
  • A primitive function type interface for Preflow.
  • A compilation bug fixed
File size: 23.2 KB
Line 
1/* -*- C++ -*-
2 * src/lemon/preflow.h - Part of LEMON, a generic C++ optimization library
3 *
4 * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
5 * (Egervary Combinatorial Optimization Research Group, 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/// 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 setCap and \ref
46  ///setFlow.
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 _G 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 setCap and \ref
144    ///setFlow, 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, though 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, though 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 \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.
292    void phase2()
293    {
294
295      int k=_node_num-2;  //bound on the highest level under n containing a node
296      int b=k;    //bound on the highest level under n of an active node
297
298   
299      VecNode first(_node_num, INVALID);
300      NNMap next(*_g, INVALID);
301      level.set(_source,0);
302      std::queue<Node> bfs_queue;
303      bfs_queue.push(_source);
304
305      while ( !bfs_queue.empty() ) {
306
307        Node v=bfs_queue.front();
308        bfs_queue.pop();
309        int l=level[v]+1;
310
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 ) {
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
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 ) {
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      }
337      b=_node_num-2;
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
364    /// of the target node \c t. This value equals to the value of
365    /// the maximum flow already after running \ref phase1.
366    Num flowValue() const {
367      return excess[_target];
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
376    ///\ref phase1.  \pre M should be a bool-valued node-map. \pre
377    ///If \ref minCut() is called after \ref phase2() then M should
378    ///be initialized to false.
379    template<typename _CutMap>
380    void minCut(_CutMap& M) const {
381      switch ( status ) {
382        case AFTER_PREFLOW_PHASE_1:
383        for(NodeIt v(*_g); v!=INVALID; ++v) {
384          if (level[v] < _node_num) {
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;
410      M.set(_source,true);
411      queue.push(_source);
412     
413      while (!queue.empty()) {
414        Node w=queue.front();
415        queue.pop();
416       
417        for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
418          Node v=_g->target(e);
419          if (!M[v] && (*_flow)[e] < (*_capacity)[e] ) {
420            queue.push(v);
421            M.set(v, true);
422          }
423        }
424       
425        for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
426          Node v=_g->source(e);
427          if (!M[v] && (*_flow)[e] > 0 ) {
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.
440    ///\pre \ref phase2() or run() should already be run.
441    template<typename _CutMap>
442    void maxMinCut(_CutMap& M) const {
443
444      for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true);
445
446      std::queue<Node> queue;
447
448      M.set(_target,false);
449      queue.push(_target);
450
451      while (!queue.empty()) {
452        Node w=queue.front();
453        queue.pop();
454
455        for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
456          Node v=_g->source(e);
457          if (M[v] && (*_flow)[e] < (*_capacity)[e] ) {
458            queue.push(v);
459            M.set(v, false);
460          }
461        }
462
463        for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
464          Node v=_g->target(e);
465          if (M[v] && (*_flow)[e] > 0 ) {
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    ///
477    void source(Node _s) {
478      _source=_s;
479      if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
480      status=AFTER_NOTHING;
481    }
482
483    ///Returns the source node.
484
485    ///Returns the source node.
486    ///
487    Node source() const {
488      return _source;
489    }
490
491    ///Sets the target node to \c _t.
492
493    ///Sets the target node to \c _t.
494    ///
495    void target(Node _t) {
496      _target=_t;
497      if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
498      status=AFTER_NOTHING;
499    }
500
501    ///Returns the target node.
502
503    ///Returns the target node.
504    ///
505    Node target() const {
506      return _target;
507    }
508
509    /// Sets the edge map of the capacities to _cap.
510
511    /// Sets the edge map of the capacities to _cap.
512    ///
513    void capacityMap(const CapacityMap& _cap) {
514      _capacity=&_cap;
515      status=AFTER_NOTHING;
516    }
517    /// Returns a reference to to capacity map.
518
519    /// Returns a reference to to capacity map.
520    ///
521    const CapacityMap &capacityMap() const {
522      return *_capacity;
523    }
524
525    /// Sets the edge map of the flows to _flow.
526
527    /// Sets the edge map of the flows to _flow.
528    ///
529    void flowMap(FlowMap& _f) {
530      _flow=&_f;
531      flow_prop=NO_FLOW;
532      status=AFTER_NOTHING;
533    }
534     
535    /// Returns a reference to to flow map.
536
537    /// Returns a reference to to flow map.
538    ///
539    const FlowMap &flowMap() const {
540      return *_flow;
541    }
542
543  private:
544
545    int push(Node w, NNMap& next, VecNode& first) {
546
547      int lev=level[w];
548      Num exc=excess[w];
549      int newlevel=_node_num;       //bound on the next level of w
550
551      for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
552        if ( (*_flow)[e] >= (*_capacity)[e] ) continue;
553        Node v=_g->target(e);
554
555        if( lev > level[v] ) { //Push is allowed now
556         
557          if ( excess[v]<=0 && v!=_target && v!=_source ) {
558            next.set(v,first[level[v]]);
559            first[level[v]]=v;
560          }
561
562          Num cap=(*_capacity)[e];
563          Num flo=(*_flow)[e];
564          Num remcap=cap-flo;
565         
566          if ( remcap >= exc ) { //A nonsaturating push.
567           
568            _flow->set(e, flo+exc);
569            excess.set(v, excess[v]+exc);
570            exc=0;
571            break;
572
573          } else { //A saturating push.
574            _flow->set(e, cap);
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 ) {
582        for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
583         
584          if( (*_flow)[e] <= 0 ) continue;
585          Node v=_g->source(e);
586
587          if( lev > level[v] ) { //Push is allowed now
588
589            if ( excess[v]<=0 && v!=_target && v!=_source ) {
590              next.set(v,first[level[v]]);
591              first[level[v]]=v;
592            }
593
594            Num flo=(*_flow)[e];
595
596            if ( flo >= exc ) { //A nonsaturating push.
597
598              _flow->set(e, flo-exc);
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;
606              _flow->set(e,0);
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    {
623      for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num);
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.
629        level.set(_target,0);
630        bfs_queue.push(_target);
631       
632        while ( !bfs_queue.empty() ) {
633         
634          Node v=bfs_queue.front();
635          bfs_queue.pop();
636          int l=level[v]+1;
637         
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 ) {
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         
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 ) {
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: 
669        for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0);
670        case ZERO_FLOW:
671        for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
672       
673        //Reverse_bfs from t, to find the starting level.
674        level.set(_target,0);
675        bfs_queue.push(_target);
676       
677        while ( !bfs_queue.empty() ) {
678         
679          Node v=bfs_queue.front();
680          bfs_queue.pop();
681          int l=level[v]+1;
682         
683          for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
684            Node w=_g->source(e);
685            if ( level[w] == _node_num && w != _source ) {
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
697        for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
698          Num c=(*_capacity)[e];
699          if ( c <= 0 ) continue;
700          Node w=_g->target(e);
701          if ( level[w] < _node_num ) {
702            if ( excess[w] <= 0 && w!=_target ) { //putting into the stack
703              next.set(w,first[level[w]]);
704              first[level[w]]=w;
705            }
706            _flow->set(e, c);
707            excess.set(w, excess[w]+c);
708          }
709        }
710        break;
711
712        case GEN_FLOW:
713        for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
714        {
715          Num exc=0;
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);
719        }
720
721        //the starting flow
722        for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e)  {
723          Num rem=(*_capacity)[e]-(*_flow)[e];
724          if ( rem <= 0 ) continue;
725          Node w=_g->target(e);
726          if ( level[w] < _node_num ) {
727            if ( excess[w] <= 0 && w!=_target ) { //putting into the stack
728              next.set(w,first[level[w]]);
729              first[level[w]]=w;
730            }   
731            _flow->set(e, (*_capacity)[e]);
732            excess.set(w, excess[w]+rem);
733          }
734        }
735       
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 ) {
741              next.set(w,first[level[w]]);
742              first[level[w]]=w;
743            } 
744            excess.set(w, excess[w]+(*_flow)[e]);
745            _flow->set(e, 0);
746          }
747        }
748        break;
749
750        case PRE_FLOW: 
751        //the starting flow
752        for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
753          Num rem=(*_capacity)[e]-(*_flow)[e];
754          if ( rem <= 0 ) continue;
755          Node w=_g->target(e);
756          if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]);
757        }
758       
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);
763        }
764       
765        //computing the excess
766        for(NodeIt w(*_g); w!=INVALID; ++w) {
767          Num exc=0;
768          for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e];
769          for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e];
770          excess.set(w,exc);
771         
772          //putting the active nodes into the stack
773          int lev=level[w];
774            if ( exc > 0 && lev < _node_num && Node(w) != _target ) {
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 ) {
818            level.set(v,_node_num);
819            v=right[v];
820          }
821          level_list[i]=INVALID;
822          if ( !what_heur ) first[i]=INVALID;
823        }
824
825        level.set(w,_node_num);
826        b=lev-1;
827        k=b;
828        //gapping ends
829
830      } else {
831
832        if ( newlevel == _node_num ) level.set(w,_node_num);
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  };
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
866} //namespace lemon
867
868#endif //LEMON_PREFLOW_H
Note: See TracBrowser for help on using the repository browser.