COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/preflow.h @ 1975:64db671eda28

Last change on this file since 1975:64db671eda28 was 1956:a055123339d5, checked in by Alpar Juttner, 18 years ago

Unified copyright notices

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