COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/preflow.h @ 1762:3915867b6975

Last change on this file since 1762:3915867b6975 was 1762:3915867b6975, checked in by jacint, 18 years ago

throwing an exception if s=t

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