COIN-OR::LEMON - Graph Library

Changeset 2514:57143c09dc20 in lemon-0.x


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)

Files:
3 added
2 deleted
8 edited

Legend:

Unmodified
Added
Removed
  • benchmark/hcube.cc

    r2391 r2514  
    113113   Graph::EdgeMap<int> flow(G);
    114114   
    115     Preflow<Graph,int> MF(G,nodes[0],nodes[1<<dim-1],map,flow);
     115    Preflow<Graph> MF(G,map,nodes[0],nodes[1<<dim-1]);
    116116    for(int i=0;i<10;i++)
    117       MF.run(MF.NO_FLOW);
     117      MF.run();
    118118  }
    119119  PrintTime("PREFLOW",T);
  • demo/disjoint_paths_demo.cc

    r2392 r2514  
    6464  Graph::EdgeMap<int> flow(graph);
    6565
    66   Preflow<Graph, int, Capacity> preflow(graph, source, target, capacity, flow);
     66  Preflow<Graph, Capacity> preflow(graph, capacity, source, target);
    6767 
    6868  preflow.run();
     
    7373    title("edge disjoint paths").scaleToA4().
    7474    copyright("(C) 2003-2007 LEMON Project").drawArrows().
    75     edgeColors(composeMap(functorMap(color), flow)).
     75    edgeColors(composeMap(functorMap(color), preflow.flowMap())).
    7676    coords(coords).run();
    7777
     
    8888  SGraph::EdgeMap<int> sflow(sgraph);
    8989
    90   Preflow<SGraph, int, SCapacity> spreflow(sgraph, SGraph::outNode(source),
    91                                            SGraph::inNode(target),
    92                                            scapacity, sflow);
     90  Preflow<SGraph, SCapacity> spreflow(sgraph, scapacity,
     91                                      SGraph::outNode(source),
     92                                      SGraph::inNode(target));
    9393
    9494  spreflow.run();
     
    100100    title("node disjoint paths").scaleToA4().
    101101    copyright("(C) 2003-2007 LEMON Project").drawArrows().
    102     edgeColors(composeMap(functorMap(color), sflow)).
     102    edgeColors(composeMap(functorMap(color), spreflow.flowMap())).
    103103    coords(SGraph::combinedNodeMap(coords,
    104104                                   shiftMap(coords,
  • demo/sub_graph_adaptor_demo.cc

    r2391 r2514  
    110110
    111111  ConstMap<Edge, int> const_1_map(1);
    112   Graph::EdgeMap<int> flow(g, 0);
    113112  // Max flow between s and t in the graph of tight edges.
    114   Preflow<SubGW, int, ConstMap<Edge, int>, Graph::EdgeMap<int> >
    115     preflow(gw, s, t, const_1_map, flow);
     113  Preflow<SubGW, ConstMap<Edge, int> > preflow(gw, const_1_map, s, t);
    116114  preflow.run();
    117115
     
    121119       << endl;
    122120  for(EdgeIt e(g); e!=INVALID; ++e)
    123     if (flow[e])
     121    if (preflow.flow(e) != 0)
    124122      cout << " " << g.id(e) << ", "
    125123           << g.id(g.source(e)) << "--"
  • doc/groups.dox

    r2500 r2514  
    224224feasible circulations.
    225225
    226 \image html flow.png
    227 \image latex flow.eps "Graph flow" width=\textwidth
     226The maximum flow problem is to find a flow between a single-source and
     227single-target that is maximum. Formally, there is \f$G=(V,A)\f$
     228directed graph, an \f$c_a:A\rightarrow\mathbf{R}^+_0\f$ capacity
     229function and given \f$s, t \in V\f$ source and target node. The
     230maximum flow is the solution of the next optimization problem:
     231
     232\f[ 0 \le f_a \le c_a \f]
     233\f[ \sum_{v\in\delta^{-}(u)}f_{vu}=\sum_{v\in\delta^{+}(u)}f_{uv} \quad u \in V \setminus \{s,t\}\f]
     234\f[ \max \sum_{v\in\delta^{+}(s)}f_{uv} - \sum_{v\in\delta^{-}(s)}f_{vu}\f]
     235
     236The lemon contains several algorithms for solve maximum flow problems:
     237- \ref lemon::EdmondsKarp "Edmonds-Karp"
     238- \ref lemon::Preflow "Goldberg's Preflow algorithm"
     239- \ref lemon::DinitzSleatorTarjan "Dinitz's blocking flow algorithm with dynamic tree"
     240- \ref lemon::GoldbergTarjan "Preflow algorithm with dynamic trees"
     241
     242In most cases the \ref lemon::Preflow "preflow" algorithm provides the
     243fastest method to compute the maximum flow. All impelementations
     244provides functions for query the minimum cut, which is the dual linear
     245programming probelm of the maximum flow.
     246
    228247*/
    229248
  • 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
  • test/preflow_test.cc

    r2391 r2514  
    2929using namespace lemon;
    3030
    31 void check_Preflow()
     31void checkPreflow()
    3232{
    3333  typedef int VType;
     
    3838  typedef concepts::ReadMap<Edge,VType> CapMap;
    3939  typedef concepts::ReadWriteMap<Edge,VType> FlowMap;
    40   typedef concepts::ReadWriteMap<Node,bool> CutMap;
     40  typedef concepts::WriteMap<Node,bool> CutMap;
    4141 
    42   typedef Preflow<Graph, int, CapMap, FlowMap> PType;
    43 
    4442  Graph g;
    4543  Node n;
     44  Edge e;
    4645  CapMap cap;
    4746  FlowMap flow;
    4847  CutMap cut;
    4948
    50   PType preflow_test(g,n,n,cap,flow);
    51 
     49  Preflow<Graph, CapMap>::DefFlowMap<FlowMap>::Create preflow_test(g,cap,n,n);
     50
     51  preflow_test.capacityMap(cap);
     52  flow = preflow_test.flowMap();
     53  preflow_test.flowMap(flow);
     54  preflow_test.source(n);
     55  preflow_test.target(n);
     56 
     57  preflow_test.init();
     58  preflow_test.flowInit(cap);
     59  preflow_test.startFirstPhase();
     60  preflow_test.startSecondPhase();
    5261  preflow_test.run();
     62  preflow_test.runMinCut();
     63
    5364  preflow_test.flowValue();
    54   preflow_test.source(n);
    55   preflow_test.flowMap(flow);
    56 
    57   preflow_test.phase1(PType::NO_FLOW);
    58   preflow_test.minCut(cut);
    59 
    60   preflow_test.phase2();
    61   preflow_test.target(n);
    62   preflow_test.capacityMap(cap);
    63   preflow_test.minMinCut(cut);
    64   preflow_test.maxMinCut(cut);
    65 }
    66 
    67 int cut_value ( SmartGraph& g, SmartGraph::NodeMap<bool>& cut,
    68                 SmartGraph::EdgeMap<int>& cap) {
     65  preflow_test.minCut(n);
     66  preflow_test.minCutMap(cut);
     67  preflow_test.flow(e);
     68
     69}
     70
     71int cutValue (const SmartGraph& g,
     72              const SmartGraph::NodeMap<bool>& cut,
     73              const SmartGraph::EdgeMap<int>& cap) {
    6974 
    7075  int c=0;
     
    7378  }
    7479  return c;
     80}
     81
     82bool checkFlow(const SmartGraph& g,
     83               const SmartGraph::EdgeMap<int>& flow,
     84               const SmartGraph::EdgeMap<int>& cap,
     85               SmartGraph::Node s, SmartGraph::Node t) {
     86
     87  for (SmartGraph::EdgeIt e(g); e != INVALID; ++e) {
     88    if (flow[e] < 0 || flow[e] > cap[e]) return false;
     89  }
     90
     91  for (SmartGraph::NodeIt n(g); n != INVALID; ++n) {
     92    if (n == s || n == t) continue;
     93    int sum = 0;
     94    for (SmartGraph::OutEdgeIt e(g, n); e != INVALID; ++e) {
     95      sum += flow[e];
     96    }
     97    for (SmartGraph::InEdgeIt e(g, n); e != INVALID; ++e) {
     98      sum -= flow[e];
     99    }
     100    if (sum != 0) return false;
     101  }
     102  return true;
    75103}
    76104
     
    86114  typedef Graph::NodeMap<bool> CutMap;
    87115
    88   typedef Preflow<Graph, int> PType;
     116  typedef Preflow<Graph, CapMap> PType;
    89117
    90118  std::string f_name;
     
    103131  readDimacs(file, g, cap, s, t);
    104132
    105   FlowMap flow(g,0);
    106 
    107  
    108 
    109   PType preflow_test(g, s, t, cap, flow);
    110   preflow_test.run(PType::ZERO_FLOW);
     133  PType preflow_test(g, cap, s, t);
     134  preflow_test.run();
    111135   
    112   CutMap min_cut(g,false);
    113   preflow_test.minCut(min_cut);
    114   int min_cut_value=cut_value(g,min_cut,cap);
    115    
    116   CutMap min_min_cut(g,false);
    117   preflow_test.minMinCut(min_min_cut);
    118   int min_min_cut_value=cut_value(g,min_min_cut,cap);
    119    
    120   CutMap max_min_cut(g,false);
    121   preflow_test.maxMinCut(max_min_cut);
    122   int max_min_cut_value=cut_value(g,max_min_cut,cap);
    123 
    124   check(preflow_test.flowValue() == min_cut_value &&
    125         min_cut_value == min_min_cut_value &&
    126         min_min_cut_value == max_min_cut_value,
     136  check(checkFlow(g, preflow_test.flowMap(), cap, s, t),
     137        "The flow is not feasible.");
     138
     139  CutMap min_cut(g);
     140  preflow_test.minCutMap(min_cut);
     141  int min_cut_value=cutValue(g,min_cut,cap);
     142 
     143  check(preflow_test.flowValue() == min_cut_value,
    127144        "The max flow value is not equal to the three min cut values.");
    128145
     146  FlowMap flow(g);
     147  flow = preflow_test.flowMap();
     148
    129149  int flow_value=preflow_test.flowValue();
    130150
    131 
    132 
    133151  for(EdgeIt e(g); e!=INVALID; ++e) cap[e]=2*cap[e];
    134   preflow_test.capacityMap(cap); 
    135 
    136   preflow_test.phase1(PType::PRE_FLOW);
    137 
    138   CutMap min_cut1(g,false);
    139   preflow_test.minCut(min_cut1);
    140   min_cut_value=cut_value(g,min_cut1,cap);
     152  preflow_test.flowInit(flow);
     153  preflow_test.startFirstPhase();
     154
     155  CutMap min_cut1(g);
     156  preflow_test.minCutMap(min_cut1);
     157  min_cut_value=cutValue(g,min_cut1,cap);
    141158   
    142159  check(preflow_test.flowValue() == min_cut_value &&
     
    144161        "The max flow value or the min cut value is wrong.");
    145162
    146   preflow_test.phase2();
    147 
    148   CutMap min_cut2(g,false);
    149   preflow_test.minCut(min_cut2);
    150   min_cut_value=cut_value(g,min_cut2,cap);
     163  preflow_test.startSecondPhase();
     164
     165  check(checkFlow(g, preflow_test.flowMap(), cap, s, t),
     166        "The flow is not feasible.");
     167
     168  CutMap min_cut2(g);
     169  preflow_test.minCutMap(min_cut2);
     170  min_cut_value=cutValue(g,min_cut2,cap);
    151171   
    152   CutMap min_min_cut2(g,false);
    153   preflow_test.minMinCut(min_min_cut2);
    154   min_min_cut_value=cut_value(g,min_min_cut2,cap);
    155  
    156   preflow_test.maxMinCut(max_min_cut);
    157   max_min_cut_value=cut_value(g,max_min_cut,cap);
    158 
    159172  check(preflow_test.flowValue() == min_cut_value &&
    160         min_cut_value == min_min_cut_value &&
    161         min_min_cut_value == max_min_cut_value &&
    162173        min_cut_value == 2*flow_value,
    163174        "The max flow value or the three min cut values were not doubled");
    164175
    165 
    166 
    167   EdgeIt e(g);
    168   for( int i=1; i==10; ++i ) {
    169     flow.set(e,0);
    170     ++e;
    171   }
    172176
    173177  preflow_test.flowMap(flow);
     
    186190  preflow_test.run();
    187191
    188   CutMap min_cut3(g,false);
    189   preflow_test.minCut(min_cut3);
    190   min_cut_value=cut_value(g,min_cut3,cap);
     192  CutMap min_cut3(g);
     193  preflow_test.minCutMap(min_cut3);
     194  min_cut_value=cutValue(g,min_cut3,cap);
    191195   
    192   CutMap min_min_cut3(g,false);
    193   preflow_test.minMinCut(min_min_cut3);
    194   min_min_cut_value=cut_value(g,min_min_cut3,cap);
    195    
    196   preflow_test.maxMinCut(max_min_cut);
    197   max_min_cut_value=cut_value(g,max_min_cut,cap);
    198 
    199   check(preflow_test.flowValue() == min_cut_value &&
    200         min_cut_value == min_min_cut_value &&
    201         min_min_cut_value == max_min_cut_value,
     196
     197  check(preflow_test.flowValue() == min_cut_value,
    202198        "The max flow value or the three min cut values are incorrect.");
    203 }
     199 
     200  return 0;
     201}
Note: See TracChangeset for help on using the changeset viewer.