COIN-OR::LEMON - Graph Library

Changeset 2514:57143c09dc20 in lemon-0.x for lemon


Ignore:
Timestamp:
11/17/07 21:58:11 (12 years ago)
Author:
Balazs Dezso
Branch:
default
Phase:
public
Convert:
svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3379
Message:

Redesign the maximum flow algorithms

Redesigned interface
Preflow changed to use elevator
Edmonds-Karp does not use the ResGraphAdaptor?
Goldberg-Tarjan algorithm (Preflow with Dynamic Trees)
Dinitz-Sleator-Tarjan (Blocking flow with Dynamic Tree)

Location:
lemon
Files:
3 added
3 edited

Legend:

Unmodified
Added
Removed
  • lemon/Makefile.am

    r2482 r2514  
    5353        lemon/dfs.h \
    5454        lemon/dijkstra.h \
     55        lemon/dinitz_sleator_tarjan.h \
    5556        lemon/dist_log.h \
    5657        lemon/dim2.h \
    5758        lemon/dimacs.h \
     59        lemon/dynamic_tree.h \
    5860        lemon/edge_set.h \
    5961        lemon/edmonds_karp.h \
     
    7274        lemon/graph_writer.h \
    7375        lemon/grid_ugraph.h \
     76        lemon/goldberg_tarjan.h \
    7477        lemon/hao_orlin.h \
    7578        lemon/hypercube_graph.h \
  • lemon/edmonds_karp.h

    r2391 r2514  
    2424/// \brief Implementation of the Edmonds-Karp algorithm.
    2525
    26 #include <lemon/graph_adaptor.h>
    2726#include <lemon/tolerance.h>
    28 #include <lemon/bfs.h>
     27#include <vector>
    2928
    3029namespace lemon {
    3130
     31  /// \brief Default traits class of EdmondsKarp class.
     32  ///
     33  /// Default traits class of EdmondsKarp class.
     34  /// \param _Graph Graph type.
     35  /// \param _CapacityMap Type of capacity map.
     36  template <typename _Graph, typename _CapacityMap>
     37  struct EdmondsKarpDefaultTraits {
     38
     39    /// \brief The graph type the algorithm runs on.
     40    typedef _Graph Graph;
     41
     42    /// \brief The type of the map that stores the edge capacities.
     43    ///
     44    /// The type of the map that stores the edge capacities.
     45    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
     46    typedef _CapacityMap CapacityMap;
     47
     48    /// \brief The type of the length of the edges.
     49    typedef typename CapacityMap::Value Value;
     50
     51    /// \brief The map type that stores the flow values.
     52    ///
     53    /// The map type that stores the flow values.
     54    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
     55    typedef typename Graph::template EdgeMap<Value> FlowMap;
     56
     57    /// \brief Instantiates a FlowMap.
     58    ///
     59    /// This function instantiates a \ref FlowMap.
     60    /// \param graph The graph, to which we would like to define the flow map.
     61    static FlowMap* createFlowMap(const Graph& graph) {
     62      return new FlowMap(graph);
     63    }
     64
     65    /// \brief The tolerance used by the algorithm
     66    ///
     67    /// The tolerance used by the algorithm to handle inexact computation.
     68    typedef Tolerance<Value> Tolerance;
     69
     70  };
     71
    3272  /// \ingroup max_flow
     73  ///
    3374  /// \brief Edmonds-Karp algorithms class.
    3475  ///
    3576  /// This class provides an implementation of the \e Edmonds-Karp \e
    3677  /// algorithm producing a flow of maximum value in a directed
    37   /// graph. The Edmonds-Karp algorithm is slower than the Preflow algorithm
    38   /// but it has an advantage of the step-by-step execution control with
    39   /// feasible flow solutions. The \e source node, the \e target node, the \e
    40   /// capacity of the edges and the \e starting \e flow value of the
    41   /// edges should be passed to the algorithm through the
    42   /// constructor.
     78  /// graphs. The Edmonds-Karp algorithm is slower than the Preflow
     79  /// algorithm but it has an advantage of the step-by-step execution
     80  /// control with feasible flow solutions. The \e source node, the \e
     81  /// target node, the \e capacity of the edges and the \e starting \e
     82  /// flow value of the edges should be passed to the algorithm
     83  /// through the constructor.
    4384  ///
    44   /// The time complexity of the algorithm is \f$ O(n * e^2) \f$ in
     85  /// The time complexity of the algorithm is \f$ O(nm^2) \f$ in
    4586  /// worst case.  Always try the preflow algorithm instead of this if
    4687  /// you just want to compute the optimal flow.
    4788  ///
    4889  /// \param _Graph The directed graph type the algorithm runs on.
    49   /// \param _Number The number type of the capacities and the flow values.
    5090  /// \param _CapacityMap The capacity map type.
    51   /// \param _FlowMap The flow map type.
    52   /// \param _Tolerance The tolerance class to handle computation problems.
     91  /// \param _Traits Traits class to set various data types used by
     92  /// the algorithm.  The default traits class is \ref
     93  /// EdmondsKarpDefaultTraits.  See \ref EdmondsKarpDefaultTraits for the
     94  /// documentation of a Edmonds-Karp traits class.
    5395  ///
    5496  /// \author Balazs Dezso
    5597#ifdef DOXYGEN
    56   template <typename _Graph, typename _Number,
    57             typename _CapacityMap, typename _FlowMap, typename _Tolerance>
    58 #else
    59   template <typename _Graph, typename _Number,
    60             typename _CapacityMap = typename _Graph::template EdgeMap<_Number>,
    61             typename _FlowMap = typename _Graph::template EdgeMap<_Number>,
    62             typename _Tolerance = Tolerance<_Number> >
     98  template <typename _Graph, typename _CapacityMap, typename _Traits>
     99#else
     100  template <typename _Graph,
     101            typename _CapacityMap = typename _Graph::template EdgeMap<int>,
     102            typename _Traits = EdmondsKarpDefaultTraits<_Graph, _CapacityMap> >
    63103#endif
    64104  class EdmondsKarp {
    65105  public:
     106
     107    typedef _Traits Traits;
     108    typedef typename Traits::Graph Graph;
     109    typedef typename Traits::CapacityMap CapacityMap;
     110    typedef typename Traits::Value Value;
     111
     112    typedef typename Traits::FlowMap FlowMap;
     113    typedef typename Traits::Tolerance Tolerance;
    66114
    67115    /// \brief \ref Exception for the case when the source equals the target.
     
    77125
    78126
    79     /// \brief The graph type the algorithm runs on.
    80     typedef _Graph Graph;
    81     /// \brief The value type of the algorithms.
    82     typedef _Number Number;
    83     /// \brief The capacity map on the edges.
    84     typedef _CapacityMap CapacityMap;
    85     /// \brief The flow map on the edges.
    86     typedef _FlowMap FlowMap;
    87     /// \brief The tolerance used by the algorithm.
    88     typedef _Tolerance Tolerance;
    89 
    90     typedef ResGraphAdaptor<Graph, Number, CapacityMap,
    91                             FlowMap, Tolerance> ResGraph;
    92 
    93127  private:
    94128
    95     typedef typename Graph::Node Node;
    96     typedef typename Graph::Edge Edge;
     129    GRAPH_TYPEDEFS(typename Graph);
     130    typedef typename Graph::template NodeMap<Edge> PredMap;
    97131   
    98     typedef typename Graph::NodeIt NodeIt;
    99     typedef typename Graph::EdgeIt EdgeIt;
    100     typedef typename Graph::InEdgeIt InEdgeIt;
    101     typedef typename Graph::OutEdgeIt OutEdgeIt;
     132    const Graph& _graph;
     133    const CapacityMap* _capacity;
     134
     135    Node _source, _target;
     136
     137    FlowMap* _flow;
     138    bool _local_flow;
     139
     140    PredMap* _pred;
     141    std::vector<Node> _queue;
     142   
     143    Tolerance _tolerance;
     144    Value _flow_value;
     145
     146    void createStructures() {
     147      if (!_flow) {
     148        _flow = Traits::createFlowMap(_graph);
     149        _local_flow = true;
     150      }
     151      if (!_pred) {
     152        _pred = new PredMap(_graph);
     153      }
     154      _queue.resize(countNodes(_graph));
     155    }
     156
     157    void destroyStructures() {
     158      if (_local_flow) {
     159        delete _flow;
     160      }
     161      if (_pred) {
     162        delete _pred;
     163      }
     164    }
    102165   
    103166  public:
    104167
     168    ///\name Named template parameters
     169
     170    ///@{
     171
     172    template <typename _FlowMap>
     173    struct DefFlowMapTraits : public Traits {
     174      typedef _FlowMap FlowMap;
     175      static FlowMap *createFlowMap(const Graph&) {
     176        throw UninitializedParameter();
     177      }
     178    };
     179
     180    /// \brief \ref named-templ-param "Named parameter" for setting
     181    /// FlowMap type
     182    ///
     183    /// \ref named-templ-param "Named parameter" for setting FlowMap
     184    /// type
     185    template <typename _FlowMap>
     186    struct DefFlowMap
     187      : public EdmondsKarp<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > {
     188      typedef EdmondsKarp<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> >
     189      Create;
     190    };
     191
     192
     193    /// @}
     194
    105195    /// \brief The constructor of the class.
    106196    ///
    107197    /// The constructor of the class.
    108198    /// \param graph The directed graph the algorithm runs on.
     199    /// \param capacity The capacity of the edges.
    109200    /// \param source The source node.
    110201    /// \param target The target node.
    111     /// \param capacity The capacity of the edges.
    112     /// \param flow The flow of the edges.
    113     /// \param tolerance Tolerance class.
    114     EdmondsKarp(const Graph& graph, Node source, Node target,
    115                 const CapacityMap& capacity, FlowMap& flow,
    116                 const Tolerance& tolerance = Tolerance())
    117       : _graph(graph), _capacity(capacity), _flow(flow),
    118         _tolerance(tolerance), _resgraph(graph, capacity, flow, tolerance),
    119         _source(source), _target(target)
     202    EdmondsKarp(const Graph& graph, const CapacityMap& capacity,
     203                Node source, Node target)
     204      : _graph(graph), _capacity(&capacity), _source(source), _target(target),
     205        _flow(0), _local_flow(false), _pred(0), _tolerance(), _flow_value()
    120206    {
    121207      if (_source == _target) {
     
    123209      }
    124210    }
     211
     212    /// \brief Destrcutor.
     213    ///
     214    /// Destructor.
     215    ~EdmondsKarp() {
     216      destroyStructures();
     217    }
     218
     219    /// \brief Sets the capacity map.
     220    ///
     221    /// Sets the capacity map.
     222    /// \return \c (*this)
     223    EdmondsKarp& capacityMap(const CapacityMap& map) {
     224      _capacity = &map;
     225      return *this;
     226    }
     227
     228    /// \brief Sets the flow map.
     229    ///
     230    /// Sets the flow map.
     231    /// \return \c (*this)
     232    EdmondsKarp& flowMap(FlowMap& map) {
     233      if (_local_flow) {
     234        delete _flow;
     235        _local_flow = false;
     236      }
     237      _flow = &map;
     238      return *this;
     239    }
     240
     241    /// \brief Returns the flow map.
     242    ///
     243    /// \return The flow map.
     244    const FlowMap& flowMap() {
     245      return *_flow;
     246    }
     247
     248    /// \brief Sets the source node.
     249    ///
     250    /// Sets the source node.
     251    /// \return \c (*this)
     252    EdmondsKarp& source(const Node& node) {
     253      _source = node;
     254      return *this;
     255    }
     256
     257    /// \brief Sets the target node.
     258    ///
     259    /// Sets the target node.
     260    /// \return \c (*this)
     261    EdmondsKarp& target(const Node& node) {
     262      _target = node;
     263      return *this;
     264    }
     265
     266    /// \brief Sets the tolerance used by algorithm.
     267    ///
     268    /// Sets the tolerance used by algorithm.
     269    EdmondsKarp& tolerance(const Tolerance& tolerance) const {
     270      _tolerance = tolerance;
     271      return *this;
     272    }
     273
     274    /// \brief Returns the tolerance used by algorithm.
     275    ///
     276    /// Returns the tolerance used by algorithm.
     277    const Tolerance& tolerance() const {
     278      return tolerance;
     279    }
     280
     281    /// \name Execution control The simplest way to execute the
     282    /// algorithm is to use the \c run() member functions.
     283    /// \n
     284    /// If you need more control on initial solution or
     285    /// execution then you have to call one \ref init() function and then
     286    /// the start() or multiple times the \c augment() member function. 
     287   
     288    ///@{
    125289
    126290    /// \brief Initializes the algorithm
     
    128292    /// It sets the flow to empty flow.
    129293    void init() {
     294      createStructures();
    130295      for (EdgeIt it(_graph); it != INVALID; ++it) {
    131         _flow.set(it, 0);
    132       }
    133       _value = 0;
     296        _flow->set(it, 0);
     297      }
     298      _flow_value = 0;
    134299    }
    135300   
    136301    /// \brief Initializes the algorithm
    137302    ///
    138     /// If the flow map initially flow this let the flow map
    139     /// unchanged but the flow value will be set by the flow
    140     /// on the outedges from the source.
    141     void flowInit() {
    142       _value = 0;
     303    /// Initializes the flow to the \c flowMap. The \c flowMap should
     304    /// contain a feasible flow, ie. in each node excluding the source
     305    /// and the target the incoming flow should be equal to the
     306    /// outgoing flow.
     307    template <typename FlowMap>
     308    void flowInit(const FlowMap& flowMap) {
     309      createStructures();
     310      for (EdgeIt e(_graph); e != INVALID; ++e) {
     311        _flow->set(e, flowMap[e]);
     312      }
     313      _flow_value = 0;
    143314      for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
    144         _value += _flow[jt];
     315        _flow_value += (*_flow)[jt];
    145316      }
    146317      for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
    147         _value -= _flow[jt];
     318        _flow_value -= (*_flow)[jt];
    148319      }
    149320    }
     
    151322    /// \brief Initializes the algorithm
    152323    ///
    153     /// If the flow map initially flow this let the flow map
    154     /// unchanged but the flow value will be set by the flow
    155     /// on the outedges from the source. It also checks that
    156     /// the flow map really contains a flow.
    157     /// \return %True when the flow map really a flow.
    158     bool checkedFlowInit() {
    159       _value = 0;
    160       for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
    161         _value += _flow[jt];
    162       }
    163       for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
    164         _value -= _flow[jt];
     324    /// Initializes the flow to the \c flowMap. The \c flowMap should
     325    /// contain a feasible flow, ie. in each node excluding the source
     326    /// and the target the incoming flow should be equal to the
     327    /// outgoing flow. 
     328    /// \return %False when the given flowMap does not contain
     329    /// feasible flow.
     330    template <typename FlowMap>
     331    bool checkedFlowInit(const FlowMap& flowMap) {
     332      createStructures();
     333      for (EdgeIt e(_graph); e != INVALID; ++e) {
     334        _flow->set(e, flowMap[e]);
    165335      }
    166336      for (NodeIt it(_graph); it != INVALID; ++it) {
    167337        if (it == _source || it == _target) continue;
    168         Number outFlow = 0;
     338        Value outFlow = 0;
    169339        for (OutEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
    170           outFlow += _flow[jt];
     340          outFlow += (*_flow)[jt];
    171341        }
    172         Number inFlow = 0;
     342        Value inFlow = 0;
    173343        for (InEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
    174           inFlow += _flow[jt];
     344          inFlow += (*_flow)[jt];
    175345        }
    176346        if (_tolerance.different(outFlow, inFlow)) {
     
    179349      }
    180350      for (EdgeIt it(_graph); it != INVALID; ++it) {
    181         if (_tolerance.less(_flow[it], 0)) return false;
    182         if (_tolerance.less(_capacity[it], _flow[it])) return false;
     351        if (_tolerance.less((*_flow)[it], 0)) return false;
     352        if (_tolerance.less((*_capacity)[it], (*_flow)[it])) return false;
     353      }
     354      _flow_value = 0;
     355      for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
     356        _flow_value += (*_flow)[jt];
     357      }
     358      for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
     359        _flow_value -= (*_flow)[jt];
    183360      }
    184361      return true;
     
    196373    /// current flow is a feasible and optimal solution.
    197374    bool augment() {
    198       typename Bfs<ResGraph>
    199       ::template DefDistMap<NullMap<Node, int> >
    200       ::Create bfs(_resgraph);
    201 
    202       NullMap<Node, int> distMap;
    203       bfs.distMap(distMap);
     375      for (NodeIt n(_graph); n != INVALID; ++n) {
     376        _pred->set(n, INVALID);
     377      }
    204378     
    205       bfs.init();
    206       bfs.addSource(_source);
    207       bfs.start(_target);
    208 
    209       if (!bfs.reached(_target)) {
    210         return false;
    211       }
    212       Number min = _resgraph.rescap(bfs.predEdge(_target));
    213       for (Node it = bfs.predNode(_target); it != _source;
    214            it = bfs.predNode(it)) {
    215         if (min > _resgraph.rescap(bfs.predEdge(it))) {
    216           min = _resgraph.rescap(bfs.predEdge(it));
    217         }
    218       }
    219       for (Node it = _target; it != _source; it = bfs.predNode(it)) {
    220         _resgraph.augment(bfs.predEdge(it), min);
    221       }
    222       _value += min;
    223       return true;
     379      int first = 0, last = 1;
     380     
     381      _queue[0] = _source;
     382      _pred->set(_source, OutEdgeIt(_graph, _source));
     383
     384      while (first != last && (*_pred)[_target] == INVALID) {
     385        Node n = _queue[first++];
     386       
     387        for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
     388          Value rem = (*_capacity)[e] - (*_flow)[e];
     389          Node t = _graph.target(e);
     390          if (_tolerance.positive(rem) && (*_pred)[t] == INVALID) {
     391            _pred->set(t, e);
     392            _queue[last++] = t;
     393          }
     394        }
     395        for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
     396          Value rem = (*_flow)[e];
     397          Node t = _graph.source(e);
     398          if (_tolerance.positive(rem) && (*_pred)[t] == INVALID) {
     399            _pred->set(t, e);
     400            _queue[last++] = t;
     401          }
     402        }
     403      }
     404
     405      if ((*_pred)[_target] != INVALID) {
     406        Node n = _target;
     407        Edge e = (*_pred)[n];
     408
     409        Value prem = (*_capacity)[e] - (*_flow)[e];
     410        n = _graph.source(e);
     411        while (n != _source) {
     412          e = (*_pred)[n];
     413          if (_graph.target(e) == n) {
     414            Value rem = (*_capacity)[e] - (*_flow)[e];
     415            if (rem < prem) prem = rem;
     416            n = _graph.source(e);
     417          } else {
     418            Value rem = (*_flow)[e];
     419            if (rem < prem) prem = rem;
     420            n = _graph.target(e);   
     421          }
     422        }
     423
     424        n = _target;
     425        e = (*_pred)[n];
     426
     427        _flow->set(e, (*_flow)[e] + prem);
     428        n = _graph.source(e);
     429        while (n != _source) {
     430          e = (*_pred)[n];
     431          if (_graph.target(e) == n) {
     432            _flow->set(e, (*_flow)[e] + prem);
     433            n = _graph.source(e);
     434          } else {
     435            _flow->set(e, (*_flow)[e] - prem);
     436            n = _graph.target(e);   
     437          }
     438        }
     439
     440        _flow_value += prem;   
     441        return true;
     442      } else {
     443        return false;
     444      }
    224445    }
    225446
     
    229450    void start() {
    230451      while (augment()) {}
    231     }
    232 
    233     /// \brief Gives back the current flow value.
    234     ///
    235     /// Gives back the current flow _value.
    236     Number flowValue() const {
    237       return _value;
    238452    }
    239453
     
    251465    }
    252466
     467    /// @}
     468
     469    /// \name Query Functions
     470    /// The result of the %Dijkstra algorithm can be obtained using these
     471    /// functions.\n
     472    /// Before the use of these functions,
     473    /// either run() or start() must be called.
     474   
     475    ///@{
     476
     477    /// \brief Returns the value of the maximum flow.
     478    ///
     479    /// Returns the value of the maximum flow by returning the excess
     480    /// of the target node \c t. This value equals to the value of
     481    /// the maximum flow already after the first phase.
     482    Value flowValue() const {
     483      return _flow_value;
     484    }
     485
     486
     487    /// \brief Returns the flow on the edge.
     488    ///
     489    /// Sets the \c flowMap to the flow on the edges. This method can
     490    /// be called after the second phase of algorithm.
     491    Value flow(const Edge& edge) const {
     492      return (*_flow)[edge];
     493    }
     494
     495    /// \brief Returns true when the node is on the source side of minimum cut.
     496    ///
     497
     498    /// Returns true when the node is on the source side of minimum
     499    /// cut. This method can be called both after running \ref
     500    /// startFirstPhase() and \ref startSecondPhase().
     501    bool minCut(const Node& node) const {
     502      return (*_pred)[node] != INVALID;
     503    }
     504
    253505    /// \brief Returns a minimum value cut.
    254506    ///
     
    257509    /// \retval cut Write node bool map.
    258510    template <typename CutMap>
    259     void minCut(CutMap& cut) const {
    260       minMinCut(cut);
    261     }
    262 
    263     /// \brief Returns the inclusionwise minimum of the minimum value cuts.
    264     ///
    265     /// Sets \c cut to the characteristic vector of the minimum value cut
    266     /// which is inclusionwise minimum. It is computed by processing a
    267     /// bfs from the source node \c source in the residual graph. 
    268     /// \retval cut Write node bool map.
    269     template <typename CutMap>
    270     void minMinCut(CutMap& cut) const {
    271 
    272       typename Bfs<ResGraph>
    273       ::template DefDistMap<NullMap<Node, int> >
    274       ::template DefProcessedMap<CutMap>
    275       ::Create bfs(_resgraph);
    276 
    277       NullMap<Node, int> distMap;
    278       bfs.distMap(distMap);
    279 
    280       bfs.processedMap(cut);
    281      
    282       bfs.run(_source);
    283     }
    284 
    285     /// \brief Returns the inclusionwise minimum of the minimum value cuts.
    286     ///
    287     /// Sets \c cut to the characteristic vector of the minimum value cut
    288     /// which is inclusionwise minimum. It is computed by processing a
    289     /// bfs from the source node \c source in the residual graph. 
    290     /// \retval cut Write node bool map.
    291     template <typename CutMap>
    292     void maxMinCut(CutMap& cut) const {
    293 
    294       typedef RevGraphAdaptor<const ResGraph> RevGraph;
    295 
    296       RevGraph revgraph(_resgraph);
    297 
    298       typename Bfs<RevGraph>
    299       ::template DefDistMap<NullMap<Node, int> >
    300       ::template DefPredMap<NullMap<Node, Edge> >
    301       ::template DefProcessedMap<NotWriteMap<CutMap> >
    302       ::Create bfs(revgraph);
    303 
    304       NullMap<Node, int> distMap;
    305       bfs.distMap(distMap);
    306 
    307       NullMap<Node, Edge> predMap;
    308       bfs.predMap(predMap);
    309 
    310       NotWriteMap<CutMap> notcut(cut);
    311       bfs.processedMap(notcut);
    312      
    313       bfs.run(_target);
    314     }
    315 
    316     /// \brief Returns the source node.
    317     ///
    318     /// Returns the source node.
    319     ///
    320     Node source() const {
    321       return _source;
    322     }
    323 
    324     /// \brief Returns the target node.
    325     ///
    326     /// Returns the target node.
    327     ///
    328     Node target() const {
    329       return _target;
    330     }
    331 
    332     /// \brief Returns a reference to capacity map.
    333     ///
    334     /// Returns a reference to capacity map.
    335     ///
    336     const CapacityMap &capacityMap() const {
    337       return *_capacity;
    338     }
    339      
    340     /// \brief Returns a reference to flow map.
    341     ///
    342     /// Returns a reference to flow map.
    343     ///
    344     const FlowMap &flowMap() const {
    345       return *_flow;
    346     }
    347 
    348     /// \brief Returns the tolerance used by algorithm.
    349     ///
    350     /// Returns the tolerance used by algorithm.
    351     const Tolerance& tolerance() const {
    352       return tolerance;
    353     }
    354    
    355   private:
    356    
    357     const Graph& _graph;
    358     const CapacityMap& _capacity;
    359     FlowMap& _flow;
    360     Tolerance _tolerance;
    361    
    362     ResGraph _resgraph;
    363     Node _source, _target;
    364     Number _value;
    365    
     511    void minCutMap(CutMap& cutMap) const {
     512      for (NodeIt n(_graph); n != INVALID; ++n) {
     513        cutMap.set(n, (*_pred)[n] != INVALID);
     514      }
     515      cutMap.set(_source, true);
     516    }   
     517
     518    /// @}
     519
    366520  };
    367521
  • lemon/preflow.h

    r2473 r2514  
    2020#define LEMON_PREFLOW_H
    2121
    22 #include <vector>
    23 #include <queue>
    24 
    2522#include <lemon/error.h>
    2623#include <lemon/bits/invalid.h>
     
    2825#include <lemon/maps.h>
    2926#include <lemon/graph_utils.h>
     27#include <lemon/elevator.h>
    3028
    3129/// \file
     
    3331/// \brief Implementation of the preflow algorithm.
    3432
    35 namespace lemon {
    36 
    37   ///\ingroup max_flow
    38   ///\brief %Preflow algorithms class.
     33namespace lemon {
     34 
     35  /// \brief Default traits class of Preflow class.
    3936  ///
    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   ///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.
     37  /// Default traits class of Preflow class.
     38  /// \param _Graph Graph type.
     39  /// \param _CapacityMap Type of capacity map.
     40  template <typename _Graph, typename _CapacityMap>
     41  struct PreflowDefaultTraits {
     42
     43    /// \brief The graph type the algorithm runs on.
     44    typedef _Graph Graph;
     45
     46    /// \brief The type of the map that stores the edge capacities.
     47    ///
     48    /// The type of the map that stores the edge capacities.
     49    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
     50    typedef _CapacityMap CapacityMap;
     51
     52    /// \brief The type of the length of the edges.
     53    typedef typename CapacityMap::Value Value;
     54
     55    /// \brief The map type that stores the flow values.
     56    ///
     57    /// The map type that stores the flow values.
     58    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
     59    typedef typename Graph::template EdgeMap<Value> FlowMap;
     60
     61    /// \brief Instantiates a FlowMap.
     62    ///
     63    /// This function instantiates a \ref FlowMap.
     64    /// \param graph The graph, to which we would like to define the flow map.
     65    static FlowMap* createFlowMap(const Graph& graph) {
     66      return new FlowMap(graph);
     67    }
     68
     69    /// \brief The eleavator type used by Preflow algorithm.
     70    ///
     71    /// The elevator type used by Preflow algorithm.
     72    ///
     73    /// \sa Elevator
     74    /// \sa LinkedElevator
     75    typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
     76   
     77    /// \brief Instantiates an Elevator.
     78    ///
     79    /// This function instantiates a \ref Elevator.
     80    /// \param graph The graph, to which we would like to define the elevator.
     81    /// \param max_level The maximum level of the elevator.
     82    static Elevator* createElevator(const Graph& graph, int max_level) {
     83      return new Elevator(graph, max_level);
     84    }
     85
     86    /// \brief The tolerance used by the algorithm
     87    ///
     88    /// The tolerance used by the algorithm to handle inexact computation.
     89    typedef Tolerance<Value> Tolerance;
     90
     91  };
     92 
     93
     94  /// \ingroup max_flow
    4995  ///
    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.)
     96  /// \brief %Preflow algorithms class.
    5797  ///
    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   ///\param Tol The tolerance type.
     98  /// This class provides an implementation of the Goldberg's \e
     99  /// preflow \e algorithm producing a flow of maximum value in a
     100  /// directed graph. The preflow algorithms are the fastest known max
     101  /// flow algorithms. The current implementation use a mixture of the
     102  /// \e "highest label" and the \e "bound decrease" heuristics.
     103  /// The worst case time complexity of the algorithm is \f$O(n^3)\f$.
    63104  ///
    64   ///\author Jacint Szabo
    65   ///\todo Second template parameter is superfluous
    66   template <typename Graph, typename Num,
    67             typename CapacityMap=typename Graph::template EdgeMap<Num>,
    68             typename FlowMap=typename Graph::template EdgeMap<Num>,
    69             typename Tol=Tolerance<Num> >
     105  /// The algorithm consists from two phases. After the first phase
     106  /// the maximal flow value and the minimum cut can be obtained. The
     107  /// second phase constructs the feasible maximum flow on each edge.
     108  ///
     109  /// \param _Graph The directed graph type the algorithm runs on.
     110  /// \param _CapacityMap The flow map type.
     111  /// \param _Traits Traits class to set various data types used by
     112  /// the algorithm.  The default traits class is \ref
     113  /// PreflowDefaultTraits.  See \ref PreflowDefaultTraits for the
     114  /// documentation of a %Preflow traits class.
     115  ///
     116  ///\author Jacint Szabo and Balazs Dezso
     117#ifdef DOXYGEN
     118  template <typename _Graph, typename _CapacityMap, typename _Traits>
     119#else
     120  template <typename _Graph,
     121            typename _CapacityMap = typename _Graph::template EdgeMap<int>,
     122            typename _Traits = PreflowDefaultTraits<_Graph, _CapacityMap> >
     123#endif
    70124  class Preflow {
    71   protected:
    72     typedef typename Graph::Node Node;
    73     typedef typename Graph::NodeIt NodeIt;
    74     typedef typename Graph::EdgeIt EdgeIt;
    75     typedef typename Graph::OutEdgeIt OutEdgeIt;
    76     typedef typename Graph::InEdgeIt InEdgeIt;
    77 
    78     typedef typename Graph::template NodeMap<Node> NNMap;
    79     typedef typename std::vector<Node> VecNode;
    80 
    81     const Graph* _g;
    82     Node _source;
    83     Node _target;
    84     const CapacityMap* _capacity;
    85     FlowMap* _flow;
    86 
    87     Tol _surely;
     125  public:
    88126   
    89     int _node_num;      //the number of nodes of G
    90    
    91     typename Graph::template NodeMap<int> level; 
    92     typename Graph::template NodeMap<Num> excess;
    93 
    94     // constants used for heuristics
    95     static const int H0=20;
    96     static const int H1=1;
    97 
    98   public:
    99 
    100     ///\ref Exception for the case when s=t.
    101 
    102     ///\ref Exception for the case when the source equals the target.
    103     class InvalidArgument : public lemon::LogicError {
     127    typedef _Traits Traits;
     128    typedef typename Traits::Graph Graph;
     129    typedef typename Traits::CapacityMap CapacityMap;
     130    typedef typename Traits::Value Value;
     131
     132    typedef typename Traits::FlowMap FlowMap;
     133    typedef typename Traits::Elevator Elevator;
     134    typedef typename Traits::Tolerance Tolerance;
     135
     136    /// \brief \ref Exception for uninitialized parameters.
     137    ///
     138    /// This error represents problems in the initialization
     139    /// of the parameters of the algorithms.
     140    class UninitializedParameter : public lemon::UninitializedParameter {
    104141    public:
    105142      virtual const char* what() const throw() {
    106         return "lemon::Preflow::InvalidArgument";
     143        return "lemon::Preflow::UninitializedParameter";
    107144      }
    108145    };
     146
     147  private:
     148
     149    GRAPH_TYPEDEFS(typename Graph);
     150
     151    const Graph& _graph;
     152    const CapacityMap* _capacity;
     153
     154    int _node_num;
     155
     156    Node _source, _target;
     157
     158    FlowMap* _flow;
     159    bool _local_flow;
     160
     161    Elevator* _level;
     162    bool _local_level;
     163
     164    typedef typename Graph::template NodeMap<Value> ExcessMap;
     165    ExcessMap* _excess;
     166
     167    Tolerance _tolerance;
     168
     169    bool _phase;
     170
     171    void createStructures() {
     172      _node_num = countNodes(_graph);
     173
     174      if (!_flow) {
     175        _flow = Traits::createFlowMap(_graph);
     176        _local_flow = true;
     177      }
     178      if (!_level) {
     179        _level = Traits::createElevator(_graph, _node_num);
     180        _local_level = true;
     181      }
     182      if (!_excess) {
     183        _excess = new ExcessMap(_graph);
     184      }
     185    }
     186
     187    void destroyStructures() {
     188      if (_local_flow) {
     189        delete _flow;
     190      }
     191      if (_local_level) {
     192        delete _level;
     193      }
     194      if (_excess) {
     195        delete _excess;
     196      }
     197    }
     198
     199  public:
     200
     201    typedef Preflow Create;
     202
     203    ///\name Named template parameters
     204
     205    ///@{
     206
     207    template <typename _FlowMap>
     208    struct DefFlowMapTraits : public Traits {
     209      typedef _FlowMap FlowMap;
     210      static FlowMap *createFlowMap(const Graph&) {
     211        throw UninitializedParameter();
     212      }
     213    };
     214
     215    /// \brief \ref named-templ-param "Named parameter" for setting
     216    /// FlowMap type
     217    ///
     218    /// \ref named-templ-param "Named parameter" for setting FlowMap
     219    /// type
     220    template <typename _FlowMap>
     221    struct DefFlowMap
     222      : public Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > {
     223      typedef Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > Create;
     224    };
     225
     226    template <typename _Elevator>
     227    struct DefElevatorTraits : public Traits {
     228      typedef _Elevator Elevator;
     229      static Elevator *createElevator(const Graph&, int) {
     230        throw UninitializedParameter();
     231      }
     232    };
     233
     234    /// \brief \ref named-templ-param "Named parameter" for setting
     235    /// Elevator type
     236    ///
     237    /// \ref named-templ-param "Named parameter" for setting Elevator
     238    /// type
     239    template <typename _Elevator>
     240    struct DefElevator
     241      : public Preflow<Graph, CapacityMap, DefElevatorTraits<_Elevator> > {
     242      typedef Preflow<Graph, CapacityMap, DefElevatorTraits<_Elevator> > Create;
     243    };
     244
     245    template <typename _Elevator>
     246    struct DefStandardElevatorTraits : public Traits {
     247      typedef _Elevator Elevator;
     248      static Elevator *createElevator(const Graph& graph, int max_level) {
     249        return new Elevator(graph, max_level);
     250      }
     251    };
     252
     253    /// \brief \ref named-templ-param "Named parameter" for setting
     254    /// Elevator type
     255    ///
     256    /// \ref named-templ-param "Named parameter" for setting Elevator
     257    /// type. The Elevator should be standard constructor interface, ie.
     258    /// the graph and the maximum level should be passed to it.
     259    template <typename _Elevator>
     260    struct DefStandardElevator
     261      : public Preflow<Graph, CapacityMap,
     262                       DefStandardElevatorTraits<_Elevator> > {
     263      typedef Preflow<Graph, CapacityMap,
     264                      DefStandardElevatorTraits<_Elevator> > Create;
     265    };   
     266
     267    /// @}
     268
     269    /// \brief The constructor of the class.
     270    ///
     271    /// The constructor of the class.
     272    /// \param graph The directed graph the algorithm runs on.
     273    /// \param capacity The capacity of the edges.
     274    /// \param source The source node.
     275    /// \param target The target node.
     276    Preflow(const Graph& graph, const CapacityMap& capacity,
     277               Node source, Node target)
     278      : _graph(graph), _capacity(&capacity),
     279        _node_num(0), _source(source), _target(target),
     280        _flow(0), _local_flow(false),
     281        _level(0), _local_level(false),
     282        _excess(0), _tolerance(), _phase() {}
     283
     284    /// \brief Destrcutor.
     285    ///
     286    /// Destructor.
     287    ~Preflow() {
     288      destroyStructures();
     289    }
     290
     291    /// \brief Sets the capacity map.
     292    ///
     293    /// Sets the capacity map.
     294    /// \return \c (*this)
     295    Preflow& capacityMap(const CapacityMap& map) {
     296      _capacity = &map;
     297      return *this;
     298    }
     299
     300    /// \brief Sets the flow map.
     301    ///
     302    /// Sets the flow map.
     303    /// \return \c (*this)
     304    Preflow& flowMap(FlowMap& map) {
     305      if (_local_flow) {
     306        delete _flow;
     307        _local_flow = false;
     308      }
     309      _flow = &map;
     310      return *this;
     311    }
     312
     313    /// \brief Returns the flow map.
     314    ///
     315    /// \return The flow map.
     316    const FlowMap& flowMap() {
     317      return *_flow;
     318    }
     319
     320    /// \brief Sets the elevator.
     321    ///
     322    /// Sets the elevator.
     323    /// \return \c (*this)
     324    Preflow& elevator(Elevator& elevator) {
     325      if (_local_level) {
     326        delete _level;
     327        _local_level = false;
     328      }
     329      _level = &elevator;
     330      return *this;
     331    }
     332
     333    /// \brief Returns the elevator.
     334    ///
     335    /// \return The elevator.
     336    const Elevator& elevator() {
     337      return *_level;
     338    }
     339
     340    /// \brief Sets the source node.
     341    ///
     342    /// Sets the source node.
     343    /// \return \c (*this)
     344    Preflow& source(const Node& node) {
     345      _source = node;
     346      return *this;
     347    }
     348
     349    /// \brief Sets the target node.
     350    ///
     351    /// Sets the target node.
     352    /// \return \c (*this)
     353    Preflow& target(const Node& node) {
     354      _target = node;
     355      return *this;
     356    }
     357 
     358    /// \brief Sets the tolerance used by algorithm.
     359    ///
     360    /// Sets the tolerance used by algorithm.
     361    Preflow& tolerance(const Tolerance& tolerance) const {
     362      _tolerance = tolerance;
     363      return *this;
     364    }
     365
     366    /// \brief Returns the tolerance used by algorithm.
     367    ///
     368    /// Returns the tolerance used by algorithm.
     369    const Tolerance& tolerance() const {
     370      return tolerance;
     371    }
     372
     373    /// \name Execution control The simplest way to execute the
     374    /// algorithm is to use one of the member functions called \c
     375    /// run(). 
     376    /// \n
     377    /// If you need more control on initial solution or
     378    /// execution then you have to call one \ref init() function and then
     379    /// the startFirstPhase() and if you need the startSecondPhase(). 
    109380   
     381    ///@{
     382
     383    /// \brief Initializes the internal data structures.
     384    ///
     385    /// Initializes the internal data structures.
     386    ///
     387    void init() {
     388      createStructures();
     389
     390      _phase = true;
     391      for (NodeIt n(_graph); n != INVALID; ++n) {
     392        _excess->set(n, 0);
     393      }
     394
     395      for (EdgeIt e(_graph); e != INVALID; ++e) {
     396        _flow->set(e, 0);
     397      }
     398
     399      typename Graph::template NodeMap<bool> reached(_graph, false);
     400
     401      _level->initStart();
     402      _level->initAddItem(_target);
     403
     404      std::vector<Node> queue;
     405      reached.set(_source, true);
     406
     407      queue.push_back(_target);
     408      reached.set(_target, true);
     409      while (!queue.empty()) {
     410        std::vector<Node> nqueue;
     411        for (int i = 0; i < int(queue.size()); ++i) {
     412          Node n = queue[i];
     413          for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
     414            Node u = _graph.source(e);
     415            if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
     416              reached.set(u, true);
     417              _level->initAddItem(u);
     418              nqueue.push_back(u);
     419            }
     420          }
     421        }
     422        queue.swap(nqueue);
     423      }
     424      _level->initFinish();
     425
     426      for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
     427        if (_tolerance.positive((*_capacity)[e])) {
     428          Node u = _graph.target(e);
     429          if ((*_level)[u] == _level->maxLevel()) continue;
     430          _flow->set(e, (*_capacity)[e]);
     431          _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
     432          if (u != _target && !_level->active(u)) {
     433            _level->activate(u);
     434          }
     435        }
     436      }
     437    }
     438
     439    /// \brief Initializes the internal data structures.
     440    ///
     441    /// Initializes the internal data structures and sets the initial
     442    /// flow to the given \c flowMap. The \c flowMap should contain a
     443    /// flow or at least a preflow, ie. in each node excluding the
     444    /// target the incoming flow should greater or equal to the
     445    /// outgoing flow.
     446    /// \return %False when the given \c flowMap is not a preflow.
     447    template <typename FlowMap>
     448    bool flowInit(const FlowMap& flowMap) {
     449      createStructures();
     450
     451      for (EdgeIt e(_graph); e != INVALID; ++e) {
     452        _flow->set(e, flowMap[e]);
     453      }
     454
     455      for (NodeIt n(_graph); n != INVALID; ++n) {
     456        Value excess = 0;
     457        for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
     458          excess += (*_flow)[e];
     459        }
     460        for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
     461          excess -= (*_flow)[e];
     462        }
     463        if (excess < 0 && n != _source) return false;
     464        _excess->set(n, excess);
     465      }
     466
     467      typename Graph::template NodeMap<bool> reached(_graph, false);
     468
     469      _level->initStart();
     470      _level->initAddItem(_target);
     471
     472      std::vector<Node> queue;
     473      reached.set(_source, true);
     474
     475      queue.push_back(_target);
     476      reached.set(_target, true);
     477      while (!queue.empty()) {
     478        _level->initNewLevel();
     479        std::vector<Node> nqueue;
     480        for (int i = 0; i < int(queue.size()); ++i) {
     481          Node n = queue[i];
     482          for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
     483            Node u = _graph.source(e);
     484            if (!reached[u] &&
     485                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
     486              reached.set(u, true);
     487              _level->initAddItem(u);
     488              nqueue.push_back(u);
     489            }
     490          }
     491          for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
     492            Node v = _graph.target(e);
     493            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
     494              reached.set(v, true);
     495              _level->initAddItem(v);
     496              nqueue.push_back(v);
     497            }
     498          }
     499        }
     500        queue.swap(nqueue);
     501      }
     502      _level->initFinish();
     503
     504      for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
     505        Value rem = (*_capacity)[e] - (*_flow)[e];
     506        if (_tolerance.positive(rem)) {
     507          Node u = _graph.target(e);
     508          if ((*_level)[u] == _level->maxLevel()) continue;
     509          _flow->set(e, (*_capacity)[e]);
     510          _excess->set(u, (*_excess)[u] + rem);
     511          if (u != _target && !_level->active(u)) {
     512            _level->activate(u);
     513          }
     514        }
     515      }
     516      for (InEdgeIt e(_graph, _source); e != INVALID; ++e) {
     517        Value rem = (*_flow)[e];
     518        if (_tolerance.positive(rem)) {
     519          Node v = _graph.source(e);
     520          if ((*_level)[v] == _level->maxLevel()) continue;
     521          _flow->set(e, 0);
     522          _excess->set(v, (*_excess)[v] + rem);
     523          if (v != _target && !_level->active(v)) {
     524            _level->activate(v);
     525          }
     526        }
     527      }
     528      return true;
     529    }
     530
     531    /// \brief Starts the first phase of the preflow algorithm.
     532    ///
     533    /// The preflow algorithm consists of two phases, this method runs
     534    /// the first phase. After the first phase the maximum flow value
     535    /// and a minimum value cut can already be computed, although a
     536    /// maximum flow is not yet obtained. So after calling this method
     537    /// \ref flowValue() returns the value of a maximum flow and \ref
     538    /// minCut() returns a minimum cut.     
     539    /// \pre One of the \ref init() functions should be called.
     540    void startFirstPhase() {
     541      _phase = true;
     542
     543      Node n = _level->highestActive();
     544      int level = _level->highestActiveLevel();
     545      while (n != INVALID) {
     546        int num = _node_num;
     547
     548        while (num > 0 && n != INVALID) {
     549          Value excess = (*_excess)[n];
     550          int new_level = _level->maxLevel();
     551
     552          for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
     553            Value rem = (*_capacity)[e] - (*_flow)[e];
     554            if (!_tolerance.positive(rem)) continue;
     555            Node v = _graph.target(e);
     556            if ((*_level)[v] < level) {
     557              if (!_level->active(v) && v != _target) {
     558                _level->activate(v);
     559              }
     560              if (!_tolerance.less(rem, excess)) {
     561                _flow->set(e, (*_flow)[e] + excess);
     562                _excess->set(v, (*_excess)[v] + excess);
     563                excess = 0;
     564                goto no_more_push_1;
     565              } else {
     566                excess -= rem;
     567                _excess->set(v, (*_excess)[v] + rem);
     568                _flow->set(e, (*_capacity)[e]);
     569              }
     570            } else if (new_level > (*_level)[v]) {
     571              new_level = (*_level)[v];
     572            }
     573          }
     574
     575          for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
     576            Value rem = (*_flow)[e];
     577            if (!_tolerance.positive(rem)) continue;
     578            Node v = _graph.source(e);
     579            if ((*_level)[v] < level) {
     580              if (!_level->active(v) && v != _target) {
     581                _level->activate(v);
     582              }
     583              if (!_tolerance.less(rem, excess)) {
     584                _flow->set(e, (*_flow)[e] - excess);
     585                _excess->set(v, (*_excess)[v] + excess);
     586                excess = 0;
     587                goto no_more_push_1;
     588              } else {
     589                excess -= rem;
     590                _excess->set(v, (*_excess)[v] + rem);
     591                _flow->set(e, 0);
     592              }
     593            } else if (new_level > (*_level)[v]) {
     594              new_level = (*_level)[v];
     595            }
     596          }
     597
     598        no_more_push_1:
     599
     600          _excess->set(n, excess);
     601
     602          if (excess != 0) {
     603            if (new_level + 1 < _level->maxLevel()) {
     604              _level->liftHighestActive(new_level + 1);
     605            } else {
     606              _level->liftHighestActiveToTop();
     607            }
     608            if (_level->emptyLevel(level)) {
     609              _level->liftToTop(level);
     610            }
     611          } else {
     612            _level->deactivate(n);
     613          }
     614         
     615          n = _level->highestActive();
     616          level = _level->highestActiveLevel();
     617          --num;
     618        }
     619       
     620        num = _node_num * 20;
     621        while (num > 0 && n != INVALID) {
     622          Value excess = (*_excess)[n];
     623          int new_level = _level->maxLevel();
     624
     625          for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
     626            Value rem = (*_capacity)[e] - (*_flow)[e];
     627            if (!_tolerance.positive(rem)) continue;
     628            Node v = _graph.target(e);
     629            if ((*_level)[v] < level) {
     630              if (!_level->active(v) && v != _target) {
     631                _level->activate(v);
     632              }
     633              if (!_tolerance.less(rem, excess)) {
     634                _flow->set(e, (*_flow)[e] + excess);
     635                _excess->set(v, (*_excess)[v] + excess);
     636                excess = 0;
     637                goto no_more_push_2;
     638              } else {
     639                excess -= rem;
     640                _excess->set(v, (*_excess)[v] + rem);
     641                _flow->set(e, (*_capacity)[e]);
     642              }
     643            } else if (new_level > (*_level)[v]) {
     644              new_level = (*_level)[v];
     645            }
     646          }
     647
     648          for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
     649            Value rem = (*_flow)[e];
     650            if (!_tolerance.positive(rem)) continue;
     651            Node v = _graph.source(e);
     652            if ((*_level)[v] < level) {
     653              if (!_level->active(v) && v != _target) {
     654                _level->activate(v);
     655              }
     656              if (!_tolerance.less(rem, excess)) {
     657                _flow->set(e, (*_flow)[e] - excess);
     658                _excess->set(v, (*_excess)[v] + excess);
     659                excess = 0;
     660                goto no_more_push_2;
     661              } else {
     662                excess -= rem;
     663                _excess->set(v, (*_excess)[v] + rem);
     664                _flow->set(e, 0);
     665              }
     666            } else if (new_level > (*_level)[v]) {
     667              new_level = (*_level)[v];
     668            }
     669          }
     670
     671        no_more_push_2:
     672
     673          _excess->set(n, excess);
     674
     675          if (excess != 0) {
     676            if (new_level + 1 < _level->maxLevel()) {
     677              _level->liftActiveOn(level, new_level + 1);
     678            } else {
     679              _level->liftActiveToTop(level);
     680            }
     681            if (_level->emptyLevel(level)) {
     682              _level->liftToTop(level);
     683            }
     684          } else {
     685            _level->deactivate(n);
     686          }
     687
     688          while (level >= 0 && _level->activeFree(level)) {
     689            --level;
     690          }
     691          if (level == -1) {
     692            n = _level->highestActive();
     693            level = _level->highestActiveLevel();
     694          } else {
     695            n = _level->activeOn(level);
     696          }
     697          --num;
     698        }
     699      }
     700    }
     701
     702    /// \brief Starts the second phase of the preflow algorithm.
     703    ///
     704    /// The preflow algorithm consists of two phases, this method runs
     705    /// the second phase. After calling \ref init() and \ref
     706    /// startFirstPhase() and then \ref startSecondPhase(), \ref
     707    /// flowMap() return a maximum flow, \ref flowValue() returns the
     708    /// value of a maximum flow, \ref minCut() returns a minimum cut
     709    /// \pre The \ref init() and startFirstPhase() functions should be
     710    /// called before.
     711    void startSecondPhase() {
     712      _phase = false;
     713
     714      typename Graph::template NodeMap<bool> reached(_graph, false);
     715      for (NodeIt n(_graph); n != INVALID; ++n) {
     716        reached.set(n, (*_level)[n] < _level->maxLevel());
     717      }
     718
     719      _level->initStart();
     720      _level->initAddItem(_source);
     721 
     722      std::vector<Node> queue;
     723      queue.push_back(_source);
     724      reached.set(_source, true);
     725
     726      while (!queue.empty()) {
     727        _level->initNewLevel();
     728        std::vector<Node> nqueue;
     729        for (int i = 0; i < int(queue.size()); ++i) {
     730          Node n = queue[i];
     731          for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
     732            Node v = _graph.target(e);
     733            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
     734              reached.set(v, true);
     735              _level->initAddItem(v);
     736              nqueue.push_back(v);
     737            }
     738          }
     739          for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
     740            Node u = _graph.source(e);
     741            if (!reached[u] &&
     742                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
     743              reached.set(u, true);
     744              _level->initAddItem(u);
     745              nqueue.push_back(u);
     746            }
     747          }
     748        }
     749        queue.swap(nqueue);
     750      }
     751      _level->initFinish();
     752
     753      for (NodeIt n(_graph); n != INVALID; ++n) {
     754        if ((*_excess)[n] > 0 && _target != n) {
     755          _level->activate(n);
     756        }
     757      }
     758
     759      Node n;
     760      while ((n = _level->highestActive()) != INVALID) {
     761        Value excess = (*_excess)[n];
     762        int level = _level->highestActiveLevel();
     763        int new_level = _level->maxLevel();
     764
     765        for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
     766          Value rem = (*_capacity)[e] - (*_flow)[e];
     767          if (!_tolerance.positive(rem)) continue;
     768          Node v = _graph.target(e);
     769          if ((*_level)[v] < level) {
     770            if (!_level->active(v) && v != _source) {
     771              _level->activate(v);
     772            }
     773            if (!_tolerance.less(rem, excess)) {
     774              _flow->set(e, (*_flow)[e] + excess);
     775              _excess->set(v, (*_excess)[v] + excess);
     776              excess = 0;
     777              goto no_more_push;
     778            } else {
     779              excess -= rem;
     780              _excess->set(v, (*_excess)[v] + rem);
     781              _flow->set(e, (*_capacity)[e]);
     782            }
     783          } else if (new_level > (*_level)[v]) {
     784            new_level = (*_level)[v];
     785          }
     786        }
     787
     788        for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
     789          Value rem = (*_flow)[e];
     790          if (!_tolerance.positive(rem)) continue;
     791          Node v = _graph.source(e);
     792          if ((*_level)[v] < level) {
     793            if (!_level->active(v) && v != _source) {
     794              _level->activate(v);
     795            }
     796            if (!_tolerance.less(rem, excess)) {
     797              _flow->set(e, (*_flow)[e] - excess);
     798              _excess->set(v, (*_excess)[v] + excess);
     799              excess = 0;
     800              goto no_more_push;
     801            } else {
     802              excess -= rem;
     803              _excess->set(v, (*_excess)[v] + rem);
     804              _flow->set(e, 0);
     805            }
     806          } else if (new_level > (*_level)[v]) {
     807            new_level = (*_level)[v];
     808          }
     809        }
     810
     811      no_more_push:
     812
     813        _excess->set(n, excess);
     814     
     815        if (excess != 0) {
     816          if (new_level + 1 < _level->maxLevel()) {
     817            _level->liftHighestActive(new_level + 1);
     818          } else {
     819            // Calculation error
     820            _level->liftHighestActiveToTop();
     821          }
     822          if (_level->emptyLevel(level)) {
     823            // Calculation error
     824            _level->liftToTop(level);
     825          }
     826        } else {
     827          _level->deactivate(n);
     828        }
     829
     830      }
     831    }
     832
     833    /// \brief Runs the preflow algorithm. 
     834    ///
     835    /// Runs the preflow algorithm.
     836    /// \note pf.run() is just a shortcut of the following code.
     837    /// \code
     838    ///   pf.init();
     839    ///   pf.startFirstPhase();
     840    ///   pf.startSecondPhase();
     841    /// \endcode
     842    void run() {
     843      init();
     844      startFirstPhase();
     845      startSecondPhase();
     846    }
     847
     848    /// \brief Runs the preflow algorithm to compute the minimum cut. 
     849    ///
     850    /// Runs the preflow algorithm to compute the minimum cut.
     851    /// \note pf.runMinCut() is just a shortcut of the following code.
     852    /// \code
     853    ///   pf.init();
     854    ///   pf.startFirstPhase();
     855    /// \endcode
     856    void runMinCut() {
     857      init();
     858      startFirstPhase();
     859    }
     860
     861    /// @}
     862
     863    /// \name Query Functions
     864    /// The result of the %Dijkstra algorithm can be obtained using these
     865    /// functions.\n
     866    /// Before the use of these functions,
     867    /// either run() or start() must be called.
    110868   
    111     ///Indicates the property of the starting flow map.
    112    
    113     ///Indicates the property of the starting flow map.
    114     ///
    115     enum FlowEnum{
    116       ///indicates an unspecified edge map. \c flow will be
    117       ///set to the constant zero flow in the beginning of
    118       ///the algorithm in this case.
    119       NO_FLOW,
    120       ///constant zero flow
    121       ZERO_FLOW,
    122       ///any flow, i.e. the sum of the in-flows equals to
    123       ///the sum of the out-flows in every node except the \c source and
    124       ///the \c target.
    125       GEN_FLOW,
    126       ///any preflow, i.e. the sum of the in-flows is at
    127       ///least the sum of the out-flows in every node except the \c source.
    128       PRE_FLOW
    129     };
    130 
    131     ///Indicates the state of the preflow algorithm.
    132 
    133     ///Indicates the state of the preflow algorithm.
    134     ///
    135     enum StatusEnum {
    136       ///before running the algorithm or
    137       ///at an unspecified state.
    138       AFTER_NOTHING,
    139       ///right after running \ref phase1()
    140       AFTER_PREFLOW_PHASE_1,     
    141       ///after running \ref phase2()
    142       AFTER_PREFLOW_PHASE_2
    143     };
    144    
    145   protected:
    146     FlowEnum flow_prop;
    147     StatusEnum status; // Do not needle this flag only if necessary.
    148    
    149   public:
    150     ///The constructor of the class.
    151 
    152     ///The constructor of the class.
    153     ///\param _gr The directed graph the algorithm runs on.
    154     ///\param _s The source node.
    155     ///\param _t The target node.
    156     ///\param _cap The capacity of the edges.
    157     ///\param _f The flow of the edges.
    158     ///\param _sr Tol class.
    159     ///Except the graph, all of these parameters can be reset by
    160     ///calling \ref source, \ref target, \ref capacityMap and \ref
    161     ///flowMap, resp.
    162     Preflow(const Graph& _gr, Node _s, Node _t,
    163             const CapacityMap& _cap, FlowMap& _f,
    164             const Tol &_sr=Tol()) :
    165         _g(&_gr), _source(_s), _target(_t), _capacity(&_cap),
    166         _flow(&_f), _surely(_sr),
    167         _node_num(countNodes(_gr)), level(_gr), excess(_gr,0),
    168         flow_prop(NO_FLOW), status(AFTER_NOTHING) {
    169         if ( _source==_target )
    170           throw InvalidArgument();
    171     }
    172    
    173     ///Give a reference to the tolerance handler class
    174 
    175     ///Give a reference to the tolerance handler class
    176     ///\sa Tolerance
    177     Tol &tolerance() { return _surely; }
    178 
    179     ///Runs the preflow algorithm. 
    180 
    181     ///Runs the preflow algorithm.
    182     ///
    183     void run() {
    184       phase1(flow_prop);
    185       phase2();
    186     }
    187    
    188     ///Runs the preflow algorithm. 
    189    
    190     ///Runs the preflow algorithm.
    191     ///\pre The starting flow map must be
    192     /// - a constant zero flow if \c fp is \c ZERO_FLOW,
    193     /// - an arbitrary flow if \c fp is \c GEN_FLOW,
    194     /// - an arbitrary preflow if \c fp is \c PRE_FLOW,
    195     /// - any map if \c fp is NO_FLOW.
    196     ///If the starting flow map is a flow or a preflow then
    197     ///the algorithm terminates faster.
    198     void run(FlowEnum fp) {
    199       flow_prop=fp;
    200       run();
    201     }
    202      
    203     ///Runs the first phase of the preflow algorithm.
    204 
    205     ///The preflow algorithm consists of two phases, this method runs
    206     ///the first phase. After the first phase the maximum flow value
    207     ///and a minimum value cut can already be computed, although a
    208     ///maximum flow is not yet obtained. So after calling this method
    209     ///\ref flowValue returns the value of a maximum flow and \ref
    210     ///minCut returns a minimum cut.     
    211     ///\warning \ref minMinCut and \ref maxMinCut do not give minimum
    212     ///value cuts unless calling \ref phase2. 
    213     ///\warning A real flow map (i.e. not \ref lemon::NullMap "NullMap")
    214     ///is needed for this phase.
    215     ///\pre The starting flow must be
    216     ///- a constant zero flow if \c fp is \c ZERO_FLOW,
    217     ///- an arbitary flow if \c fp is \c GEN_FLOW,
    218     ///- an arbitary preflow if \c fp is \c PRE_FLOW,
    219     ///- any map if \c fp is NO_FLOW.
    220     void phase1(FlowEnum fp)
    221     {
    222       flow_prop=fp;
    223       phase1();
    224     }
    225 
    226    
    227     ///Runs the first phase of the preflow algorithm.
    228 
    229     ///The preflow algorithm consists of two phases, this method runs
    230     ///the first phase. After the first phase the maximum flow value
    231     ///and a minimum value cut can already be computed, although a
    232     ///maximum flow is not yet obtained. So after calling this method
    233     ///\ref flowValue returns the value of a maximum flow and \ref
    234     ///minCut returns a minimum cut.
    235     ///\warning \ref minMinCut() and \ref maxMinCut() do not
    236     ///give minimum value cuts unless calling \ref phase2().
    237     ///\warning A real flow map (i.e. not \ref lemon::NullMap "NullMap")
    238     ///is needed for this phase.
    239     void phase1()
    240     {
    241       int heur0=int(H0*_node_num);  //time while running 'bound decrease'
    242       int heur1=int(H1*_node_num);  //time while running 'highest label'
    243       int heur=heur1;         //starting time interval (#of relabels)
    244       int numrelabel=0;
    245 
    246       bool what_heur=1;
    247       //It is 0 in case 'bound decrease' and 1 in case 'highest label'
    248 
    249       bool end=false;
    250       //Needed for 'bound decrease', true means no active
    251       //nodes are above bound b.
    252 
    253       int k=_node_num-2;  //bound on the highest level under n containing a node
    254       int b=k;    //bound on the highest level under n containing an active node
    255 
    256       VecNode first(_node_num, INVALID);
    257       NNMap next(*_g, INVALID);
    258 
    259       NNMap left(*_g, INVALID);
    260       NNMap right(*_g, INVALID);
    261       VecNode level_list(_node_num,INVALID);
    262       //List of the nodes in level i<n, set to n.
    263 
    264       preflowPreproc(first, next, level_list, left, right);
    265 
    266       //Push/relabel on the highest level active nodes.
    267       while ( true ) {
    268         if ( b == 0 ) {
    269           if ( !what_heur && !end && k > 0 ) {
    270             b=k;
    271             end=true;
    272           } else break;
    273         }
    274 
    275         if ( first[b]==INVALID ) --b;
    276         else {
    277           end=false;
    278           Node w=first[b];
    279           first[b]=next[w];
    280           int newlevel=push(w, next, first);
    281           if ( excess[w] != 0 ) {
    282             relabel(w, newlevel, first, next, level_list,
    283                     left, right, b, k, what_heur);
    284           }
    285 
    286           ++numrelabel;
    287           if ( numrelabel >= heur ) {
    288             numrelabel=0;
    289             if ( what_heur ) {
    290               what_heur=0;
    291               heur=heur0;
    292               end=false;
    293             } else {
    294               what_heur=1;
    295               heur=heur1;
    296               b=k;
    297             }
    298           }
    299         }
    300       }
    301       flow_prop=PRE_FLOW;
    302       status=AFTER_PREFLOW_PHASE_1;
    303     }
    304     // Heuristics:
    305     //   2 phase
    306     //   gap
    307     //   list 'level_list' on the nodes on level i implemented by hand
    308     //   stack 'active' on the active nodes on level i     
    309     //   runs heuristic 'highest label' for H1*n relabels
    310     //   runs heuristic 'bound decrease' for H0*n relabels,
    311     //        starts with 'highest label'
    312     //   Parameters H0 and H1 are initialized to 20 and 1.
    313 
    314 
    315     ///Runs the second phase of the preflow algorithm.
    316 
    317     ///The preflow algorithm consists of two phases, this method runs
    318     ///the second phase. After calling \ref phase1() and then
    319     ///\ref phase2(),
    320     /// \ref flowMap() return a maximum flow, \ref flowValue
    321     ///returns the value of a maximum flow, \ref minCut returns a
    322     ///minimum cut, while the methods \ref minMinCut and \ref
    323     ///maxMinCut return the inclusionwise minimum and maximum cuts of
    324     ///minimum value, resp.  \pre \ref phase1 must be called before.
    325     ///
    326     /// \todo The inexact computation can cause positive excess on a set of
    327     /// unpushable nodes. We may have to watch the empty level in this case
    328     /// due to avoid the terrible long running time.
    329     void phase2()
    330     {
    331 
    332       int k=_node_num-2;  //bound on the highest level under n containing a node
    333       int b=k;    //bound on the highest level under n of an active node
    334 
    335    
    336       VecNode first(_node_num, INVALID);
    337       NNMap next(*_g, INVALID);
    338       level.set(_source,0);
    339       std::queue<Node> bfs_queue;
    340       bfs_queue.push(_source);
    341 
    342       while ( !bfs_queue.empty() ) {
    343 
    344         Node v=bfs_queue.front();
    345         bfs_queue.pop();
    346         int l=level[v]+1;
    347 
    348         for(InEdgeIt e(*_g,v); e!=INVALID; ++e) {
    349           if ( !_surely.positive((*_capacity)[e] - (*_flow)[e])) continue;
    350           Node u=_g->source(e);
    351           if ( level[u] >= _node_num ) {
    352             bfs_queue.push(u);
    353             level.set(u, l);
    354             if ( excess[u] != 0 ) {
    355               next.set(u,first[l]);
    356               first[l]=u;
    357             }
    358           }
    359         }
    360 
    361         for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) {
    362           if ( !_surely.positive((*_flow)[e]) ) continue;
    363           Node u=_g->target(e);
    364           if ( level[u] >= _node_num ) {
    365             bfs_queue.push(u);
    366             level.set(u, l);
    367             if ( excess[u] != 0 ) {
    368               next.set(u,first[l]);
    369               first[l]=u;
    370             }
    371           }
    372         }
    373       }
    374       b=_node_num-2;
    375 
    376       while ( true ) {
    377 
    378         if ( b == 0 ) break;
    379         if ( first[b]==INVALID ) --b;
    380         else {
    381           Node w=first[b];
    382           first[b]=next[w];
    383           int newlevel=push(w,next, first);
    384          
    385           if ( newlevel == _node_num) {
    386             excess.set(w, 0);
    387             level.set(w,_node_num);
    388           }
    389           //relabel
    390           if ( excess[w] != 0 ) {
    391             level.set(w,++newlevel);
    392             next.set(w,first[newlevel]);
    393             first[newlevel]=w;
    394             b=newlevel;
    395           }
    396         }
    397       } // while(true)
    398       flow_prop=GEN_FLOW;
    399       status=AFTER_PREFLOW_PHASE_2;
    400     }
    401 
    402     /// Returns the value of the maximum flow.
    403 
     869    ///@{
     870
     871    /// \brief Returns the value of the maximum flow.
     872    ///
    404873    /// Returns the value of the maximum flow by returning the excess
    405874    /// of the target node \c t. This value equals to the value of
    406     /// the maximum flow already after running \ref phase1.
    407     Num flowValue() const {
    408       return excess[_target];
    409     }
    410 
    411 
    412     ///Returns a minimum value cut.
    413 
    414     ///Sets \c M to the characteristic vector of a minimum value
    415     ///cut. This method can be called both after running \ref
    416     ///phase1 and \ref phase2. It is much faster after
    417     ///\ref phase1.  \pre M should be a bool-valued node-map. \pre
    418     ///If \ref minCut() is called after \ref phase2() then M should
    419     ///be initialized to false.
    420     template<typename _CutMap>
    421     void minCut(_CutMap& M) const {
    422       switch ( status ) {
    423         case AFTER_PREFLOW_PHASE_1:
    424         for(NodeIt v(*_g); v!=INVALID; ++v) {
    425           if (level[v] < _node_num) {
    426             M.set(v, false);
    427           } else {
    428             M.set(v, true);
    429           }
    430         }
    431         break;
    432         case AFTER_PREFLOW_PHASE_2:
    433         minMinCut(M);
    434         break;
    435         case AFTER_NOTHING:
    436         break;
    437       }
    438     }
    439 
    440     ///Returns the inclusionwise minimum of the minimum value cuts.
    441 
    442     ///Sets \c M to the characteristic vector of the minimum value cut
    443     ///which is inclusionwise minimum. It is computed by processing a
    444     ///bfs from the source node \c s in the residual graph.  \pre M
    445     ///should be a node map of bools initialized to false.  \pre \ref
    446     ///phase2 should already be run.
    447     template<typename _CutMap>
    448     void minMinCut(_CutMap& M) const {
    449 
    450       std::queue<Node> queue;
    451       M.set(_source,true);
    452       queue.push(_source);
    453      
    454       while (!queue.empty()) {
    455         Node w=queue.front();
    456         queue.pop();
    457        
    458         for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
    459           Node v=_g->target(e);
    460           if (!M[v] && _surely.positive((*_capacity)[e] -(*_flow)[e]) ) {
    461             queue.push(v);
    462             M.set(v, true);
    463           }
    464         }
    465        
    466         for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
    467           Node v=_g->source(e);
    468           if (!M[v] && _surely.positive((*_flow)[e]) ) {
    469             queue.push(v);
    470             M.set(v, true);
    471           }
    472         }
    473       }
     875    /// the maximum flow already after the first phase.
     876    Value flowValue() const {
     877      return (*_excess)[_target];
     878    }
     879
     880    /// \brief Returns true when the node is on the source side of minimum cut.
     881    ///
     882    /// Returns true when the node is on the source side of minimum
     883    /// cut. This method can be called both after running \ref
     884    /// startFirstPhase() and \ref startSecondPhase().
     885    bool minCut(const Node& node) const {
     886      return ((*_level)[node] == _level->maxLevel()) == _phase;
     887    }
     888 
     889    /// \brief Returns a minimum value cut.
     890    ///
     891    /// Sets the \c cutMap to the characteristic vector of a minimum value
     892    /// cut. This method can be called both after running \ref
     893    /// startFirstPhase() and \ref startSecondPhase(). The result after second
     894    /// phase could be changed slightly if inexact computation is used.
     895    /// \pre The \c cutMap should be a bool-valued node-map.
     896    template <typename CutMap>
     897    void minCutMap(CutMap& cutMap) const {
     898      for (NodeIt n(_graph); n != INVALID; ++n) {
     899        cutMap.set(n, minCut(n));
     900      }
     901    }
     902
     903    /// \brief Returns the flow on the edge.
     904    ///
     905    /// Sets the \c flowMap to the flow on the edges. This method can
     906    /// be called after the second phase of algorithm.
     907    Value flow(const Edge& edge) const {
     908      return (*_flow)[edge];
    474909    }
    475910   
    476     ///Returns the inclusionwise maximum of the minimum value cuts.
    477 
    478     ///Sets \c M to the characteristic vector of the minimum value cut
    479     ///which is inclusionwise maximum. It is computed by processing a
    480     ///backward bfs from the target node \c t in the residual graph.
    481     ///\pre \ref phase2() or run() should already be run.
    482     template<typename _CutMap>
    483     void maxMinCut(_CutMap& M) const {
    484 
    485       for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true);
    486 
    487       std::queue<Node> queue;
    488 
    489       M.set(_target,false);
    490       queue.push(_target);
    491 
    492       while (!queue.empty()) {
    493         Node w=queue.front();
    494         queue.pop();
    495 
    496         for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
    497           Node v=_g->source(e);
    498           if (M[v] && _surely.positive((*_capacity)[e] - (*_flow)[e]) ) {
    499             queue.push(v);
    500             M.set(v, false);
    501           }
    502         }
    503 
    504         for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
    505           Node v=_g->target(e);
    506           if (M[v] && _surely.positive((*_flow)[e]) ) {
    507             queue.push(v);
    508             M.set(v, false);
    509           }
    510         }
    511       }
    512     }
    513 
    514     ///Sets the source node to \c _s.
    515 
    516     ///Sets the source node to \c _s.
    517     ///
    518     void source(Node _s) {
    519       _source=_s;
    520       if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
    521       status=AFTER_NOTHING;
    522     }
    523 
    524     ///Returns the source node.
    525 
    526     ///Returns the source node.
    527     ///
    528     Node source() const {
    529       return _source;
    530     }
    531 
    532     ///Sets the target node to \c _t.
    533 
    534     ///Sets the target node to \c _t.
    535     ///
    536     void target(Node _t) {
    537       _target=_t;
    538       if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
    539       status=AFTER_NOTHING;
    540     }
    541 
    542     ///Returns the target node.
    543 
    544     ///Returns the target node.
    545     ///
    546     Node target() const {
    547       return _target;
    548     }
    549 
    550     /// Sets the edge map of the capacities to _cap.
    551 
    552     /// Sets the edge map of the capacities to _cap.
    553     ///
    554     void capacityMap(const CapacityMap& _cap) {
    555       _capacity=&_cap;
    556       status=AFTER_NOTHING;
    557     }
    558     /// Returns a reference to capacity map.
    559 
    560     /// Returns a reference to capacity map.
    561     ///
    562     const CapacityMap &capacityMap() const {
    563       return *_capacity;
    564     }
    565 
    566     /// Sets the edge map of the flows to _flow.
    567 
    568     /// Sets the edge map of the flows to _flow.
    569     ///
    570     void flowMap(FlowMap& _f) {
    571       _flow=&_f;
    572       flow_prop=NO_FLOW;
    573       status=AFTER_NOTHING;
    574     }
    575      
    576     /// Returns a reference to flow map.
    577 
    578     /// Returns a reference to flow map.
    579     ///
    580     const FlowMap &flowMap() const {
    581       return *_flow;
    582     }
    583 
    584   private:
    585 
    586     int push(Node w, NNMap& next, VecNode& first) {
    587 
    588       int lev=level[w];
    589       Num exc=excess[w];
    590       int newlevel=_node_num;       //bound on the next level of w
    591 
    592       for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
    593         if ( !_surely.positive((*_capacity)[e] - (*_flow)[e])) continue;
    594         Node v=_g->target(e);
    595        
    596         if( lev > level[v] ) { //Push is allowed now
    597          
    598           if ( excess[v] == 0 && v!=_target && v!=_source ) {
    599             next.set(v,first[level[v]]);
    600             first[level[v]]=v;
    601           }
    602 
    603           Num cap=(*_capacity)[e];
    604           Num flo=(*_flow)[e];
    605           Num remcap=cap-flo;
    606          
    607           if ( ! _surely.less(remcap, exc) ) { //A nonsaturating push.
    608            
    609             _flow->set(e, flo+exc);
    610             excess.set(v, excess[v]+exc);
    611             exc=0;
    612             break;
    613 
    614           } else { //A saturating push.
    615             _flow->set(e, cap);
    616             excess.set(v, excess[v]+remcap);
    617             exc-=remcap;
    618           }
    619         } else if ( newlevel > level[v] ) newlevel = level[v];
    620       } //for out edges wv
    621 
    622       if ( exc != 0 ) {
    623         for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
    624          
    625           if ( !_surely.positive((*_flow)[e]) ) continue;
    626           Node v=_g->source(e);
    627          
    628           if( lev > level[v] ) { //Push is allowed now
    629 
    630             if ( excess[v] == 0 && v!=_target && v!=_source ) {
    631               next.set(v,first[level[v]]);
    632               first[level[v]]=v;
    633             }
    634 
    635             Num flo=(*_flow)[e];
    636 
    637             if ( !_surely.less(flo, exc) ) { //A nonsaturating push.
    638 
    639               _flow->set(e, flo-exc);
    640               excess.set(v, excess[v]+exc);
    641               exc=0;
    642               break;
    643             } else {  //A saturating push.
    644 
    645               excess.set(v, excess[v]+flo);
    646               exc-=flo;
    647               _flow->set(e,0);
    648             }
    649           } else if ( newlevel > level[v] ) newlevel = level[v];
    650         } //for in edges vw
    651 
    652       } // if w still has excess after the out edge for cycle
    653 
    654       excess.set(w, exc);
    655      
    656       return newlevel;
    657     }
    658    
    659    
    660    
    661     void preflowPreproc(VecNode& first, NNMap& next,
    662                         VecNode& level_list, NNMap& left, NNMap& right)
    663     {
    664       for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num);
    665       std::queue<Node> bfs_queue;
    666      
    667       if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
    668         //Reverse_bfs from t in the residual graph,
    669         //to find the starting level.
    670         level.set(_target,0);
    671         bfs_queue.push(_target);
    672        
    673         while ( !bfs_queue.empty() ) {
    674          
    675           Node v=bfs_queue.front();
    676           bfs_queue.pop();
    677           int l=level[v]+1;
    678          
    679           for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
    680             if ( !_surely.positive((*_capacity)[e] - (*_flow)[e] )) continue;
    681             Node w=_g->source(e);
    682             if ( level[w] == _node_num && w != _source ) {
    683               bfs_queue.push(w);
    684               Node z=level_list[l];
    685               if ( z!=INVALID ) left.set(z,w);
    686               right.set(w,z);
    687               level_list[l]=w;
    688               level.set(w, l);
    689             }
    690           }
    691          
    692           for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
    693             if ( !_surely.positive((*_flow)[e]) ) continue;
    694             Node w=_g->target(e);
    695             if ( level[w] == _node_num && w != _source ) {
    696               bfs_queue.push(w);
    697               Node z=level_list[l];
    698               if ( z!=INVALID ) left.set(z,w);
    699               right.set(w,z);
    700               level_list[l]=w;
    701               level.set(w, l);
    702             }
    703           }
    704         } //while
    705       } //if
    706 
    707 
    708       switch (flow_prop) {
    709         case NO_FLOW: 
    710         for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0);
    711         case ZERO_FLOW:
    712         for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
    713        
    714         //Reverse_bfs from t, to find the starting level.
    715         level.set(_target,0);
    716         bfs_queue.push(_target);
    717        
    718         while ( !bfs_queue.empty() ) {
    719          
    720           Node v=bfs_queue.front();
    721           bfs_queue.pop();
    722           int l=level[v]+1;
    723          
    724           for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
    725             Node w=_g->source(e);
    726             if ( level[w] == _node_num && w != _source ) {
    727               bfs_queue.push(w);
    728               Node z=level_list[l];
    729               if ( z!=INVALID ) left.set(z,w);
    730               right.set(w,z);
    731               level_list[l]=w;
    732               level.set(w, l);
    733             }
    734           }
    735         }
    736        
    737         //the starting flow
    738         for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
    739           Num c=(*_capacity)[e];
    740           if ( !_surely.positive(c) ) 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, c);
    748             excess.set(w, excess[w]+c);
    749           }
    750         }
    751         break;
    752 
    753         case GEN_FLOW:
    754         for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
    755         {
    756           Num exc=0;
    757           for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e];
    758           for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e];
    759           if (!_surely.positive(exc)) {
    760             exc = 0;
    761           }
    762           excess.set(_target,exc);
    763         }
    764 
    765         //the starting flow
    766         for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e)  {
    767           Num rem=(*_capacity)[e]-(*_flow)[e];
    768           if ( !_surely.positive(rem) ) continue;
    769           Node w=_g->target(e);
    770           if ( level[w] < _node_num ) {
    771             if ( excess[w] == 0 && w!=_target ) { //putting into the stack
    772               next.set(w,first[level[w]]);
    773               first[level[w]]=w;
    774             }   
    775             _flow->set(e, (*_capacity)[e]);
    776             excess.set(w, excess[w]+rem);
    777           }
    778         }
    779        
    780         for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) {
    781           if ( !_surely.positive((*_flow)[e]) ) continue;
    782           Node w=_g->source(e);
    783           if ( level[w] < _node_num ) {
    784             if ( excess[w] == 0 && w!=_target ) {
    785               next.set(w,first[level[w]]);
    786               first[level[w]]=w;
    787             } 
    788             excess.set(w, excess[w]+(*_flow)[e]);
    789             _flow->set(e, 0);
    790           }
    791         }
    792         break;
    793 
    794         case PRE_FLOW: 
    795         //the starting flow
    796         for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
    797           Num rem=(*_capacity)[e]-(*_flow)[e];
    798           if ( !_surely.positive(rem) ) continue;
    799           Node w=_g->target(e);
    800           if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]);
    801         }
    802        
    803         for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
    804           if ( !_surely.positive((*_flow)[e]) ) continue;
    805           Node w=_g->source(e);
    806           if ( level[w] < _node_num ) _flow->set(e, 0);
    807         }
    808        
    809         //computing the excess
    810         for(NodeIt w(*_g); w!=INVALID; ++w) {
    811           Num exc=0;
    812           for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e];
    813           for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e];
    814           if (!_surely.positive(exc)) {
    815             exc = 0;
    816           }
    817           excess.set(w,exc);
    818          
    819           //putting the active nodes into the stack
    820           int lev=level[w];
    821             if ( exc != 0 && lev < _node_num && Node(w) != _target ) {
    822               next.set(w,first[lev]);
    823               first[lev]=w;
    824             }
    825         }
    826         break;
    827       } //switch
    828     } //preflowPreproc
    829 
    830 
    831     void relabel(Node w, int newlevel, VecNode& first, NNMap& next,
    832                  VecNode& level_list, NNMap& left,
    833                  NNMap& right, int& b, int& k, bool what_heur )
    834     {
    835 
    836       int lev=level[w];
    837 
    838       Node right_n=right[w];
    839       Node left_n=left[w];
    840 
    841       //unlacing starts
    842       if ( right_n!=INVALID ) {
    843         if ( left_n!=INVALID ) {
    844           right.set(left_n, right_n);
    845           left.set(right_n, left_n);
    846         } else {
    847           level_list[lev]=right_n;
    848           left.set(right_n, INVALID);
    849         }
    850       } else {
    851         if ( left_n!=INVALID ) {
    852           right.set(left_n, INVALID);
    853         } else {
    854           level_list[lev]=INVALID;
    855         }
    856       }
    857       //unlacing ends
    858 
    859       if ( level_list[lev]==INVALID ) {
    860 
    861         //gapping starts
    862         for (int i=lev; i!=k ; ) {
    863           Node v=level_list[++i];
    864           while ( v!=INVALID ) {
    865             level.set(v,_node_num);
    866             v=right[v];
    867           }
    868           level_list[i]=INVALID;
    869           if ( !what_heur ) first[i]=INVALID;
    870         }
    871 
    872         level.set(w,_node_num);
    873         b=lev-1;
    874         k=b;
    875         //gapping ends
    876 
    877       } else {
    878 
    879         if ( newlevel == _node_num ) level.set(w,_node_num);
    880         else {
    881           level.set(w,++newlevel);
    882           next.set(w,first[newlevel]);
    883           first[newlevel]=w;
    884           if ( what_heur ) b=newlevel;
    885           if ( k < newlevel ) ++k;      //now k=newlevel
    886           Node z=level_list[newlevel];
    887           if ( z!=INVALID ) left.set(z,w);
    888           right.set(w,z);
    889           left.set(w,INVALID);
    890           level_list[newlevel]=w;
    891         }
    892       }
    893     } //relabel
    894 
    895   };
    896 
    897   ///\ingroup max_flow
    898   ///\brief Function type interface for Preflow algorithm.
    899   ///
    900   ///Function type interface for Preflow algorithm.
    901   ///\sa Preflow
    902   template<class GR, class CM, class FM>
    903   Preflow<GR,typename CM::Value,CM,FM> preflow(const GR &g,
    904                             typename GR::Node source,
    905                             typename GR::Node target,
    906                             const CM &cap,
    907                             FM &flow
    908                             )
    909   {
    910     return Preflow<GR,typename CM::Value,CM,FM>(g,source,target,cap,flow);
    911   }
    912 
    913 } //namespace lemon
    914 
    915 #endif //LEMON_PREFLOW_H
     911    /// @}   
     912  };
     913}
     914
     915#endif
Note: See TracChangeset for help on using the changeset viewer.