COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/lemon/preflow.h @ 1222:a3fb216a267d

Last change on this file since 1222:a3fb216a267d was 1222:a3fb216a267d, checked in by Alpar Juttner, 15 years ago

The first step toward function type interface to Preflow alg:

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