COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/preflow.h @ 1933:a876a3d6a4c7

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