LEMON/LEMON-main Changeset - r581:aa1804409f29

Repositories » LEMON » LEMON-main

Summary
Changelog
Files
Switch To
- loading...
Options
- Lightweight changelog
- Search
Pull Requests

Changeset - r581:aa1804409f29

« 580:2313ed

582:7a28e2 »

r581:aa1804409f29

Tue, 14 Apr 2009 10:35:38

[Not Reviewed]

0 comments (0 inline)

kpeter (Peter Kovacs)
kpeter@inf.elte.hu

Exploit that the standard maps are reference maps (#190)

0 9 0

default

9 files changed with 345 insertions and 345 deletions:

lemon/circulation.h

lemon/core.h

lemon/elevator.h

lemon/gomory_hu.h

lemon/hao_orlin.h

lemon/max_matching.h

lemon/min_cost_arborescence.h

lemon/preflow.h

test/kruskal_test.cc

↑ Collapse diff ↑

lemon/circulation.h

Show inline comments

@@ -264,490 +264,490 @@
 		    };
 		    template <typename T>
 		    struct SetStandardElevatorTraits : public Traits {
 		      typedef T Elevator;
 		      static Elevator *createElevator(const Digraph& digraph, int max_level) {
 		        return new Elevator(digraph, max_level);
+		      }
 		    };
 		    /// \brief \ref named-templ-param "Named parameter" for setting
 		    /// Elevator type with automatic allocation
 		    ///
 		    /// \ref named-templ-param "Named parameter" for setting Elevator
 		    /// type with automatic allocation.
 		    /// The Elevator should have standard constructor interface to be
 		    /// able to automatically created by the algorithm (i.e. the
 		    /// digraph and the maximum level should be passed to it).
 		    /// However an external elevator object could also be passed to the
 		    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
 		    /// before calling \ref run() or \ref init().
 		    /// \sa SetElevator
 		    template <typename T>
 		    struct SetStandardElevator
 		      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
 		                       SetStandardElevatorTraits<T> > {
 		      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
 		                      SetStandardElevatorTraits<T> > Create;
 		    };
 		    /// @}
 		  protected:
 		    Circulation() {}
 		  public:
 		    /// The constructor of the class.
 		    /// The constructor of the class.
 		    /// \param g The digraph the algorithm runs on.
 		    /// \param lo The lower bound capacity of the arcs.
 		    /// \param up The upper bound capacity of the arcs.
 		    /// \param delta The lower bound for the supply of the nodes.
 		    Circulation(const Digraph &g,const LCapMap &lo,
 		                const UCapMap &up,const DeltaMap &delta)
 		      : _g(g), _node_num(),
 		        _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
 		        _level(0), _local_level(false), _excess(0), _el() {}
 		    /// Destructor.
 		    ~Circulation() {
 		      destroyStructures();
+		    }
 		  private:
 		    void createStructures() {
 		      _node_num = _el = countNodes(_g);
 		      if (!_flow) {
 		        _flow = Traits::createFlowMap(_g);
 		        _local_flow = true;
+		      }
 		      if (!_level) {
 		        _level = Traits::createElevator(_g, _node_num);
 		        _local_level = true;
+		      }
 		      if (!_excess) {
 		        _excess = new ExcessMap(_g);
+		      }
+		    }
 		    void destroyStructures() {
 		      if (_local_flow) {
 		        delete _flow;
+		      }
 		      if (_local_level) {
 		        delete _level;
+		      }
 		      if (_excess) {
 		        delete _excess;
+		      }
+		    }
 		  public:
 		    /// Sets the lower bound capacity map.
 		    /// Sets the lower bound capacity map.
 		    /// \return <tt>(*this)</tt>
 		    Circulation& lowerCapMap(const LCapMap& map) {
 		      _lo = &map;
 		      return *this;
+		    }
 		    /// Sets the upper bound capacity map.
 		    /// Sets the upper bound capacity map.
 		    /// \return <tt>(*this)</tt>
 		    Circulation& upperCapMap(const LCapMap& map) {
 		      _up = &map;
 		      return *this;
+		    }
 		    /// Sets the lower bound map for the supply of the nodes.
 		    /// Sets the lower bound map for the supply of the nodes.
 		    /// \return <tt>(*this)</tt>
 		    Circulation& deltaMap(const DeltaMap& map) {
 		      _delta = &map;
 		      return *this;
+		    }
 		    /// \brief Sets the flow map.
 		    ///
 		    /// Sets the flow map.
 		    /// If you don't use this function before calling \ref run() or
 		    /// \ref init(), an instance will be allocated automatically.
 		    /// The destructor deallocates this automatically allocated map,
 		    /// of course.
 		    /// \return <tt>(*this)</tt>
 		    Circulation& flowMap(FlowMap& map) {
 		      if (_local_flow) {
 		        delete _flow;
 		        _local_flow = false;
+		      }
 		      _flow = &map;
 		      return *this;
+		    }
 		    /// \brief Sets the elevator used by algorithm.
 		    ///
 		    /// Sets the elevator used by algorithm.
 		    /// If you don't use this function before calling \ref run() or
 		    /// \ref init(), an instance will be allocated automatically.
 		    /// The destructor deallocates this automatically allocated elevator,
 		    /// of course.
 		    /// \return <tt>(*this)</tt>
 		    Circulation& elevator(Elevator& elevator) {
 		      if (_local_level) {
 		        delete _level;
 		        _local_level = false;
+		      }
 		      _level = &elevator;
 		      return *this;
+		    }
 		    /// \brief Returns a const reference to the elevator.
 		    ///
 		    /// Returns a const reference to the elevator.
 		    ///
 		    /// \pre Either \ref run() or \ref init() must be called before
 		    /// using this function.
 		    const Elevator& elevator() const {
 		      return *_level;
+		    }
 		    /// \brief Sets the tolerance used by algorithm.
 		    ///
 		    /// Sets the tolerance used by algorithm.
 		    Circulation& tolerance(const Tolerance& tolerance) const {
 		      _tol = tolerance;
 		      return *this;
+		    }
 		    /// \brief Returns a const reference to the tolerance.
 		    ///
 		    /// Returns a const reference to the tolerance.
 		    const Tolerance& tolerance() const {
 		      return tolerance;
+		    }
 		    /// \name Execution Control
 		    /// The simplest way to execute the algorithm is to call \ref run().\n
 		    /// If you need more control on the initial solution or the execution,
 		    /// first you have to call one of the \ref init() functions, then
 		    /// the \ref start() function.
 		    ///@{
 		    /// Initializes the internal data structures.
 		    /// Initializes the internal data structures and sets all flow values
 		    /// to the lower bound.
 		    void init()
+		    {
 		      createStructures();
 		      for(NodeIt n(_g);n!=INVALID;++n) {
-		        _excess->set(n, (*_delta)[n]);
+		        (*_excess)[n] = (*_delta)[n];
+		      }
 		      for (ArcIt e(_g);e!=INVALID;++e) {
 		        _flow->set(e, (*_lo)[e]);
 		        _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
 		        _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
 		        (*_excess)[_g.target(e)] += (*_flow)[e];
 		        (*_excess)[_g.source(e)] -= (*_flow)[e];
+		      }
 		      // global relabeling tested, but in general case it provides
 		      // worse performance for random digraphs
 		      _level->initStart();
 		      for(NodeIt n(_g);n!=INVALID;++n)
 		        _level->initAddItem(n);
 		      _level->initFinish();
 		      for(NodeIt n(_g);n!=INVALID;++n)
 		        if(_tol.positive((*_excess)[n]))
 		          _level->activate(n);
+		    }
 		    /// Initializes the internal data structures using a greedy approach.
 		    /// Initializes the internal data structures using a greedy approach
 		    /// to construct the initial solution.
 		    void greedyInit()
+		    {
 		      createStructures();
 		      for(NodeIt n(_g);n!=INVALID;++n) {
-		        _excess->set(n, (*_delta)[n]);
+		        (*_excess)[n] = (*_delta)[n];
+		      }
 		      for (ArcIt e(_g);e!=INVALID;++e) {
 		        if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
 		          _flow->set(e, (*_up)[e]);
 		          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
 		          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
 		          (*_excess)[_g.target(e)] += (*_up)[e];
 		          (*_excess)[_g.source(e)] -= (*_up)[e];
 		        } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
 		          _flow->set(e, (*_lo)[e]);
 		          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_lo)[e]);
 		          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_lo)[e]);
 		          (*_excess)[_g.target(e)] += (*_lo)[e];
 		          (*_excess)[_g.source(e)] -= (*_lo)[e];
 		        } else {
 		          Value fc = -(*_excess)[_g.target(e)];
 		          _flow->set(e, fc);
 		          _excess->set(_g.target(e), 0);
 		          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
 		          (*_excess)[_g.target(e)] = 0;
 		          (*_excess)[_g.source(e)] -= fc;
+		        }
+		      }
 		      _level->initStart();
 		      for(NodeIt n(_g);n!=INVALID;++n)
 		        _level->initAddItem(n);
 		      _level->initFinish();
 		      for(NodeIt n(_g);n!=INVALID;++n)
 		        if(_tol.positive((*_excess)[n]))
 		          _level->activate(n);
+		    }
 		    ///Executes the algorithm
 		    ///This function executes the algorithm.
 		    ///
 		    ///\return \c true if a feasible circulation is found.
 		    ///
 		    ///\sa barrier()
 		    ///\sa barrierMap()
 		    bool start()
+		    {
 		      Node act;
 		      Node bact=INVALID;
 		      Node last_activated=INVALID;
 		      while((act=_level->highestActive())!=INVALID) {
 		        int actlevel=(*_level)[act];
 		        int mlevel=_node_num;
 		        Value exc=(*_excess)[act];
 		        for(OutArcIt e(_g,act);e!=INVALID; ++e) {
 		          Node v = _g.target(e);
 		          Value fc=(*_up)[e]-(*_flow)[e];
 		          if(!_tol.positive(fc)) continue;
 		          if((*_level)[v]<actlevel) {
 		            if(!_tol.less(fc, exc)) {
 		              _flow->set(e, (*_flow)[e] + exc);
-		              _excess->set(v, (*_excess)[v] + exc);
+		              (*_excess)[v] += exc;
 		              if(!_level->active(v) && _tol.positive((*_excess)[v]))
 		                _level->activate(v);
-		              _excess->set(act,0);
+		              (*_excess)[act] = 0;
 		              _level->deactivate(act);
 		              goto next_l;
+		            }
 		            else {
 		              _flow->set(e, (*_up)[e]);
-		              _excess->set(v, (*_excess)[v] + fc);
+		              (*_excess)[v] += fc;
 		              if(!_level->active(v) && _tol.positive((*_excess)[v]))
 		                _level->activate(v);
 		              exc-=fc;
+		            }
+		          }
 		          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
+		        }
 		        for(InArcIt e(_g,act);e!=INVALID; ++e) {
 		          Node v = _g.source(e);
 		          Value fc=(*_flow)[e]-(*_lo)[e];
 		          if(!_tol.positive(fc)) continue;
 		          if((*_level)[v]<actlevel) {
 		            if(!_tol.less(fc, exc)) {
 		              _flow->set(e, (*_flow)[e] - exc);
-		              _excess->set(v, (*_excess)[v] + exc);
+		              (*_excess)[v] += exc;
 		              if(!_level->active(v) && _tol.positive((*_excess)[v]))
 		                _level->activate(v);
-		              _excess->set(act,0);
+		              (*_excess)[act] = 0;
 		              _level->deactivate(act);
 		              goto next_l;
+		            }
 		            else {
 		              _flow->set(e, (*_lo)[e]);
-		              _excess->set(v, (*_excess)[v] + fc);
+		              (*_excess)[v] += fc;
 		              if(!_level->active(v) && _tol.positive((*_excess)[v]))
 		                _level->activate(v);
 		              exc-=fc;
+		            }
+		          }
 		          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
+		        }
-		        _excess->set(act, exc);
+		        (*_excess)[act] = exc;
 		        if(!_tol.positive(exc)) _level->deactivate(act);
 		        else if(mlevel==_node_num) {
 		          _level->liftHighestActiveToTop();
 		          _el = _node_num;
 		          return false;
+		        }
 		        else {
 		          _level->liftHighestActive(mlevel+1);
 		          if(_level->onLevel(actlevel)==0) {
 		            _el = actlevel;
 		            return false;
+		          }
+		        }
 		      next_l:
+		        ;
+		      }
 		      return true;
+		    }
 		    /// Runs the algorithm.
 		    /// This function runs the algorithm.
 		    ///
 		    /// \return \c true if a feasible circulation is found.
 		    ///
 		    /// \note Apart from the return value, c.run() is just a shortcut of
 		    /// the following code.
 		    /// \code
 		    ///   c.greedyInit();
 		    ///   c.start();
 		    /// \endcode
 		    bool run() {
 		      greedyInit();
 		      return start();
+		    }
 		    /// @}
 		    /// \name Query Functions
 		    /// The results of the circulation algorithm can be obtained using
 		    /// these functions.\n
 		    /// Either \ref run() or \ref start() should be called before
 		    /// using them.
 		    ///@{
 		    /// \brief Returns the flow on the given arc.
 		    ///
 		    /// Returns the flow on the given arc.
 		    ///
 		    /// \pre Either \ref run() or \ref init() must be called before
 		    /// using this function.
 		    Value flow(const Arc& arc) const {
 		      return (*_flow)[arc];
+		    }
 		    /// \brief Returns a const reference to the flow map.
 		    ///
 		    /// Returns a const reference to the arc map storing the found flow.
 		    ///
 		    /// \pre Either \ref run() or \ref init() must be called before
 		    /// using this function.
 		    const FlowMap& flowMap() const {
 		      return *_flow;
+		    }
 		    /**
 		       \brief Returns \c true if the given node is in a barrier.
 		       Barrier is a set \e B of nodes for which
 		       \f[ \sum_{a\in\delta_{out}(B)} upper(a) -
 		           \sum_{a\in\delta_{in}(B)} lower(a) < \sum_{v\in B}delta(v) \f]
 		       holds. The existence of a set with this property prooves that a
 		       feasible circualtion cannot exist.
 		       This function returns \c true if the given node is in the found
 		       barrier. If a feasible circulation is found, the function
 		       gives back \c false for every node.
 		       \pre Either \ref run() or \ref init() must be called before
 		       using this function.
 		       \sa barrierMap()
 		       \sa checkBarrier()
 		    */
 		    bool barrier(const Node& node) const
+		    {
 		      return (*_level)[node] >= _el;
+		    }
 		    /// \brief Gives back a barrier.
 		    ///
 		    /// This function sets \c bar to the characteristic vector of the
 		    /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
 		    /// node map with \c bool (or convertible) value type.
 		    ///
 		    /// If a feasible circulation is found, the function gives back an
 		    /// empty set, so \c bar[v] will be \c false for all nodes \c v.
 		    ///
 		    /// \note This function calls \ref barrier() for each node,
 		    /// so it runs in O(n) time.
 		    ///
 		    /// \pre Either \ref run() or \ref init() must be called before
 		    /// using this function.
 		    ///
 		    /// \sa barrier()
 		    /// \sa checkBarrier()
 		    template<class BarrierMap>
 		    void barrierMap(BarrierMap &bar) const
+		    {
 		      for(NodeIt n(_g);n!=INVALID;++n)
 		        bar.set(n, (*_level)[n] >= _el);
+		    }
 		    /// @}
 		    /// \name Checker Functions
 		    /// The feasibility of the results can be checked using
 		    /// these functions.\n
 		    /// Either \ref run() or \ref start() should be called before
 		    /// using them.
 		    ///@{
 		    ///Check if the found flow is a feasible circulation
 		    ///Check if the found flow is a feasible circulation,
 		    ///
 		    bool checkFlow() const {
 		      for(ArcIt e(_g);e!=INVALID;++e)
 		        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
 		      for(NodeIt n(_g);n!=INVALID;++n)
+		        {
 		          Value dif=-(*_delta)[n];
 		          for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
 		          for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
 		          if(_tol.negative(dif)) return false;
+		        }
 		      return true;
+		    }
 		    ///Check whether or not the last execution provides a barrier
 		    ///Check whether or not the last execution provides a barrier.
 		    ///\sa barrier()
 		    ///\sa barrierMap()
 		    bool checkBarrier() const
+		    {
 		      Value delta=0;
 		      for(NodeIt n(_g);n!=INVALID;++n)
 		        if(barrier(n))
 		          delta-=(*_delta)[n];
 		      for(ArcIt e(_g);e!=INVALID;++e)
+		        {
 		          Node s=_g.source(e);
 		          Node t=_g.target(e);
 		          if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
 		          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
+		        }
 		      return _tol.negative(delta);
+		    }
 		    /// @}
 		  };
+		}
 		#endif

lemon/core.h

Show inline comments

@@ -1126,561 +1126,561 @@
 		  /// the next edge from \c u to \c v after \c prev.
 		  /// \return The found edge or \ref INVALID if there is no such an edge.
 		  ///
 		  /// Thus you can iterate through each edge between \c u and \c v
 		  /// as it follows.
 		  ///\code
 		  /// for(Edge e = findEdge(g,u,v); e != INVALID; e = findEdge(g,u,v,e)) {
 		  ///   ...
 		  /// }
 		  ///\endcode
 		  ///
 		  /// \note \ref ConEdgeIt provides iterator interface for the same
 		  /// functionality.
 		  ///
 		  ///\sa ConEdgeIt
 		  template <typename Graph>
 		  inline typename Graph::Edge
 		  findEdge(const Graph &g, typename Graph::Node u, typename Graph::Node v,
 		            typename Graph::Edge p = INVALID) {
 		    return _core_bits::FindEdgeSelector<Graph>::find(g, u, v, p);
+		  }
 		  /// \brief Iterator for iterating on parallel edges connecting the same nodes.
 		  ///
 		  /// Iterator for iterating on parallel edges connecting the same nodes.
 		  /// It is a higher level interface for the findEdge() function. You can
 		  /// use it the following way:
 		  ///\code
 		  /// for (ConEdgeIt<Graph> it(g, u, v); it != INVALID; ++it) {
 		  ///   ...
 		  /// }
 		  ///\endcode
 		  ///
 		  ///\sa findEdge()
 		  template <typename GR>
 		  class ConEdgeIt : public GR::Edge {
 		  public:
 		    typedef GR Graph;
 		    typedef typename Graph::Edge Parent;
 		    typedef typename Graph::Edge Edge;
 		    typedef typename Graph::Node Node;
 		    /// \brief Constructor.
 		    ///
 		    /// Construct a new ConEdgeIt iterating on the edges that
 		    /// connects nodes \c u and \c v.
 		    ConEdgeIt(const Graph& g, Node u, Node v) : _graph(g), _u(u), _v(v) {
 		      Parent::operator=(findEdge(_graph, _u, _v));
+		    }
 		    /// \brief Constructor.
 		    ///
 		    /// Construct a new ConEdgeIt that continues iterating from edge \c e.
 		    ConEdgeIt(const Graph& g, Edge e) : Parent(e), _graph(g) {}
 		    /// \brief Increment operator.
 		    ///
 		    /// It increments the iterator and gives back the next edge.
 		    ConEdgeIt& operator++() {
 		      Parent::operator=(findEdge(_graph, _u, _v, *this));
 		      return *this;
+		    }
 		  private:
 		    const Graph& _graph;
 		    Node _u, _v;
 		  };
 		  ///Dynamic arc look-up between given endpoints.
 		  ///Using this class, you can find an arc in a digraph from a given
 		  ///source to a given target in amortized time <em>O</em>(log<em>d</em>),
 		  ///where <em>d</em> is the out-degree of the source node.
 		  ///
 		  ///It is possible to find \e all parallel arcs between two nodes with
 		  ///the \c operator() member.
 		  ///
 		  ///This is a dynamic data structure. Consider to use \ref ArcLookUp or
 		  ///\ref AllArcLookUp if your digraph is not changed so frequently.
 		  ///
 		  ///This class uses a self-adjusting binary search tree, the Splay tree
 		  ///of Sleator and Tarjan to guarantee the logarithmic amortized
 		  ///time bound for arc look-ups. This class also guarantees the
 		  ///optimal time bound in a constant factor for any distribution of
 		  ///queries.
 		  ///
 		  ///\tparam GR The type of the underlying digraph.
 		  ///
 		  ///\sa ArcLookUp
 		  ///\sa AllArcLookUp
 		  template <typename GR>
 		  class DynArcLookUp
 		    : protected ItemSetTraits<GR, typename GR::Arc>::ItemNotifier::ObserverBase
+		  {
 		  public:
 		    typedef typename ItemSetTraits<GR, typename GR::Arc>
 		    ::ItemNotifier::ObserverBase Parent;
 		    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
 		    typedef GR Digraph;
 		  protected:
 		    class AutoNodeMap : public ItemSetTraits<GR, Node>::template Map<Arc>::Type {
 		    public:
 		      typedef typename ItemSetTraits<GR, Node>::template Map<Arc>::Type Parent;
 		      AutoNodeMap(const GR& digraph) : Parent(digraph, INVALID) {}
 		      virtual void add(const Node& node) {
 		        Parent::add(node);
 		        Parent::set(node, INVALID);
+		      }
 		      virtual void add(const std::vector<Node>& nodes) {
 		        Parent::add(nodes);
 		        for (int i = 0; i < int(nodes.size()); ++i) {
 		          Parent::set(nodes[i], INVALID);
+		        }
+		      }
 		      virtual void build() {
 		        Parent::build();
 		        Node it;
 		        typename Parent::Notifier* nf = Parent::notifier();
 		        for (nf->first(it); it != INVALID; nf->next(it)) {
 		          Parent::set(it, INVALID);
+		        }
+		      }
 		    };
 		    const Digraph &_g;
 		    AutoNodeMap _head;
 		    typename Digraph::template ArcMap<Arc> _parent;
 		    typename Digraph::template ArcMap<Arc> _left;
 		    typename Digraph::template ArcMap<Arc> _right;
 		    class ArcLess {
 		      const Digraph &g;
 		    public:
 		      ArcLess(const Digraph &_g) : g(_g) {}
 		      bool operator()(Arc a,Arc b) const
+		      {
 		        return g.target(a)<g.target(b);
+		      }
 		    };
 		  public:
 		    ///Constructor
 		    ///Constructor.
 		    ///
 		    ///It builds up the search database.
 		    DynArcLookUp(const Digraph &g)
 		      : _g(g),_head(g),_parent(g),_left(g),_right(g)
+		    {
 		      Parent::attach(_g.notifier(typename Digraph::Arc()));
 		      refresh();
+		    }
 		  protected:
 		    virtual void add(const Arc& arc) {
 		      insert(arc);
+		    }
 		    virtual void add(const std::vector<Arc>& arcs) {
 		      for (int i = 0; i < int(arcs.size()); ++i) {
 		        insert(arcs[i]);
+		      }
+		    }
 		    virtual void erase(const Arc& arc) {
 		      remove(arc);
+		    }
 		    virtual void erase(const std::vector<Arc>& arcs) {
 		      for (int i = 0; i < int(arcs.size()); ++i) {
 		        remove(arcs[i]);
+		      }
+		    }
 		    virtual void build() {
 		      refresh();
+		    }
 		    virtual void clear() {
 		      for(NodeIt n(_g);n!=INVALID;++n) {
-		        _head.set(n, INVALID);
+		        _head[n] = INVALID;
+		      }
+		    }
 		    void insert(Arc arc) {
 		      Node s = _g.source(arc);
 		      Node t = _g.target(arc);
 		      _left.set(arc, INVALID);
 		      _right.set(arc, INVALID);
 		      _left[arc] = INVALID;
 		      _right[arc] = INVALID;
 		      Arc e = _head[s];
 		      if (e == INVALID) {
 		        _head.set(s, arc);
 		        _parent.set(arc, INVALID);
 		        _head[s] = arc;
 		        _parent[arc] = INVALID;
 		        return;
+		      }
 		      while (true) {
 		        if (t < _g.target(e)) {
 		          if (_left[e] == INVALID) {
 		            _left.set(e, arc);
 		            _parent.set(arc, e);
 		            _left[e] = arc;
 		            _parent[arc] = e;
 		            splay(arc);
 		            return;
 		          } else {
 		            e = _left[e];
+		          }
 		        } else {
 		          if (_right[e] == INVALID) {
 		            _right.set(e, arc);
 		            _parent.set(arc, e);
 		            _right[e] = arc;
 		            _parent[arc] = e;
 		            splay(arc);
 		            return;
 		          } else {
 		            e = _right[e];
+		          }
+		        }
+		      }
+		    }
 		    void remove(Arc arc) {
 		      if (_left[arc] == INVALID) {
 		        if (_right[arc] != INVALID) {
-		          _parent.set(_right[arc], _parent[arc]);
+		          _parent[_right[arc]] = _parent[arc];
+		        }
 		        if (_parent[arc] != INVALID) {
 		          if (_left[_parent[arc]] == arc) {
-		            _left.set(_parent[arc], _right[arc]);
+		            _left[_parent[arc]] = _right[arc];
 		          } else {
-		            _right.set(_parent[arc], _right[arc]);
+		            _right[_parent[arc]] = _right[arc];
+		          }
 		        } else {
-		          _head.set(_g.source(arc), _right[arc]);
+		          _head[_g.source(arc)] = _right[arc];
+		        }
 		      } else if (_right[arc] == INVALID) {
-		        _parent.set(_left[arc], _parent[arc]);
+		        _parent[_left[arc]] = _parent[arc];
 		        if (_parent[arc] != INVALID) {
 		          if (_left[_parent[arc]] == arc) {
-		            _left.set(_parent[arc], _left[arc]);
+		            _left[_parent[arc]] = _left[arc];
 		          } else {
-		            _right.set(_parent[arc], _left[arc]);
+		            _right[_parent[arc]] = _left[arc];
+		          }
 		        } else {
-		          _head.set(_g.source(arc), _left[arc]);
+		          _head[_g.source(arc)] = _left[arc];
+		        }
 		      } else {
 		        Arc e = _left[arc];
 		        if (_right[e] != INVALID) {
 		          e = _right[e];
 		          while (_right[e] != INVALID) {
 		            e = _right[e];
+		          }
 		          Arc s = _parent[e];
-		          _right.set(_parent[e], _left[e]);
+		          _right[_parent[e]] = _left[e];
 		          if (_left[e] != INVALID) {
-		            _parent.set(_left[e], _parent[e]);
+		            _parent[_left[e]] = _parent[e];
+		          }
 		          _left.set(e, _left[arc]);
 		          _parent.set(_left[arc], e);
 		          _right.set(e, _right[arc]);
 		          _parent.set(_right[arc], e);
 		          _left[e] = _left[arc];
 		          _parent[_left[arc]] = e;
 		          _right[e] = _right[arc];
 		          _parent[_right[arc]] = e;
-		          _parent.set(e, _parent[arc]);
+		          _parent[e] = _parent[arc];
 		          if (_parent[arc] != INVALID) {
 		            if (_left[_parent[arc]] == arc) {
-		              _left.set(_parent[arc], e);
+		              _left[_parent[arc]] = e;
 		            } else {
-		              _right.set(_parent[arc], e);
+		              _right[_parent[arc]] = e;
+		            }
+		          }
 		          splay(s);
 		        } else {
 		          _right.set(e, _right[arc]);
 		          _parent.set(_right[arc], e);
-		          _parent.set(e, _parent[arc]);
+		          _right[e] = _right[arc];
 		          _parent[_right[arc]] = e;
 		          _parent[e] = _parent[arc];
 		          if (_parent[arc] != INVALID) {
 		            if (_left[_parent[arc]] == arc) {
-		              _left.set(_parent[arc], e);
+		              _left[_parent[arc]] = e;
 		            } else {
-		              _right.set(_parent[arc], e);
+		              _right[_parent[arc]] = e;
+		            }
 		          } else {
-		            _head.set(_g.source(arc), e);
+		            _head[_g.source(arc)] = e;
+		          }
+		        }
+		      }
+		    }
 		    Arc refreshRec(std::vector<Arc> &v,int a,int b)
+		    {
 		      int m=(a+b)/2;
 		      Arc me=v[m];
 		      if (a < m) {
 		        Arc left = refreshRec(v,a,m-1);
 		        _left.set(me, left);
 		        _parent.set(left, me);
 		        _left[me] = left;
 		        _parent[left] = me;
 		      } else {
-		        _left.set(me, INVALID);
+		        _left[me] = INVALID;
+		      }
 		      if (m < b) {
 		        Arc right = refreshRec(v,m+1,b);
 		        _right.set(me, right);
 		        _parent.set(right, me);
 		        _right[me] = right;
 		        _parent[right] = me;
 		      } else {
-		        _right.set(me, INVALID);
+		        _right[me] = INVALID;
+		      }
 		      return me;
+		    }
 		    void refresh() {
 		      for(NodeIt n(_g);n!=INVALID;++n) {
 		        std::vector<Arc> v;
 		        for(OutArcIt a(_g,n);a!=INVALID;++a) v.push_back(a);
 		        if (!v.empty()) {
 		          std::sort(v.begin(),v.end(),ArcLess(_g));
 		          Arc head = refreshRec(v,0,v.size()-1);
 		          _head.set(n, head);
 		          _parent.set(head, INVALID);
 		          _head[n] = head;
 		          _parent[head] = INVALID;
+		        }
-		        else _head.set(n, INVALID);
+		        else _head[n] = INVALID;
+		      }
+		    }
 		    void zig(Arc v) {
 		      Arc w = _parent[v];
 		      _parent.set(v, _parent[w]);
 		      _parent.set(w, v);
 		      _left.set(w, _right[v]);
 		      _right.set(v, w);
 		      _parent[v] = _parent[w];
 		      _parent[w] = v;
 		      _left[w] = _right[v];
 		      _right[v] = w;
 		      if (_parent[v] != INVALID) {
 		        if (_right[_parent[v]] == w) {
-		          _right.set(_parent[v], v);
+		          _right[_parent[v]] = v;
 		        } else {
-		          _left.set(_parent[v], v);
+		          _left[_parent[v]] = v;
+		        }
+		      }
 		      if (_left[w] != INVALID){
-		        _parent.set(_left[w], w);
+		        _parent[_left[w]] = w;
+		      }
+		    }
 		    void zag(Arc v) {
 		      Arc w = _parent[v];
 		      _parent.set(v, _parent[w]);
 		      _parent.set(w, v);
 		      _right.set(w, _left[v]);
 		      _left.set(v, w);
 		      _parent[v] = _parent[w];
 		      _parent[w] = v;
 		      _right[w] = _left[v];
 		      _left[v] = w;
 		      if (_parent[v] != INVALID){
 		        if (_left[_parent[v]] == w) {
-		          _left.set(_parent[v], v);
+		          _left[_parent[v]] = v;
 		        } else {
-		          _right.set(_parent[v], v);
+		          _right[_parent[v]] = v;
+		        }
+		      }
 		      if (_right[w] != INVALID){
-		        _parent.set(_right[w], w);
+		        _parent[_right[w]] = w;
+		      }
+		    }
 		    void splay(Arc v) {
 		      while (_parent[v] != INVALID) {
 		        if (v == _left[_parent[v]]) {
 		          if (_parent[_parent[v]] == INVALID) {
 		            zig(v);
 		          } else {
 		            if (_parent[v] == _left[_parent[_parent[v]]]) {
 		              zig(_parent[v]);
 		              zig(v);
 		            } else {
 		              zig(v);
 		              zag(v);
+		            }
+		          }
 		        } else {
 		          if (_parent[_parent[v]] == INVALID) {
 		            zag(v);
 		          } else {
 		            if (_parent[v] == _left[_parent[_parent[v]]]) {
 		              zag(v);
 		              zig(v);
 		            } else {
 		              zag(_parent[v]);
 		              zag(v);
+		            }
+		          }
+		        }
+		      }
 		      _head[_g.source(v)] = v;
+		    }
 		  public:
 		    ///Find an arc between two nodes.
 		    ///Find an arc between two nodes.
 		    ///\param s The source node.
 		    ///\param t The target node.
 		    ///\param p The previous arc between \c s and \c t. It it is INVALID or
 		    ///not given, the operator finds the first appropriate arc.
 		    ///\return An arc from \c s to \c t after \c p or
 		    ///\ref INVALID if there is no more.
 		    ///
 		    ///For example, you can count the number of arcs from \c u to \c v in the
 		    ///following way.
 		    ///\code
 		    ///DynArcLookUp<ListDigraph> ae(g);
 		    ///...
 		    ///int n = 0;
 		    ///for(Arc a = ae(u,v); a != INVALID; a = ae(u,v,a)) n++;
 		    ///\endcode
 		    ///
 		    ///Finding the arcs take at most <em>O</em>(log<em>d</em>)
 		    ///amortized time, specifically, the time complexity of the lookups
 		    ///is equal to the optimal search tree implementation for the
 		    ///current query distribution in a constant factor.
 		    ///
 		    ///\note This is a dynamic data structure, therefore the data
 		    ///structure is updated after each graph alteration. Thus although
 		    ///this data structure is theoretically faster than \ref ArcLookUp
 		    ///and \ref AllArcLookUp, it often provides worse performance than
 		    ///them.
 		    Arc operator()(Node s, Node t, Arc p = INVALID) const  {
 		      if (p == INVALID) {
 		        Arc a = _head[s];
 		        if (a == INVALID) return INVALID;
 		        Arc r = INVALID;
 		        while (true) {
 		          if (_g.target(a) < t) {
 		            if (_right[a] == INVALID) {
 		              const_cast<DynArcLookUp&>(*this).splay(a);
 		              return r;
 		            } else {
 		              a = _right[a];
+		            }
 		          } else {
 		            if (_g.target(a) == t) {
 		              r = a;
+		            }
 		            if (_left[a] == INVALID) {
 		              const_cast<DynArcLookUp&>(*this).splay(a);
 		              return r;
 		            } else {
 		              a = _left[a];
+		            }
+		          }
+		        }
 		      } else {
 		        Arc a = p;
 		        if (_right[a] != INVALID) {
 		          a = _right[a];
 		          while (_left[a] != INVALID) {
 		            a = _left[a];
+		          }
 		          const_cast<DynArcLookUp&>(*this).splay(a);
 		        } else {
 		          while (_parent[a] != INVALID && _right[_parent[a]] ==  a) {
 		            a = _parent[a];
+		          }
 		          if (_parent[a] == INVALID) {
 		            return INVALID;
 		          } else {
 		            a = _parent[a];
 		            const_cast<DynArcLookUp&>(*this).splay(a);
+		          }
+		        }
 		        if (_g.target(a) == t) return a;
 		        else return INVALID;
+		      }
+		    }
 		  };
 		  ///Fast arc look-up between given endpoints.
 		  ///Using this class, you can find an arc in a digraph from a given
 		  ///source to a given target in time <em>O</em>(log<em>d</em>),
 		  ///where <em>d</em> is the out-degree of the source node.
 		  ///
 		  ///It is not possible to find \e all parallel arcs between two nodes.
 		  ///Use \ref AllArcLookUp for this purpose.
 		  ///
 		  ///\warning This class is static, so you should call refresh() (or at
 		  ///least refresh(Node)) to refresh this data structure whenever the
 		  ///digraph changes. This is a time consuming (superlinearly proportional
 		  ///(<em>O</em>(<em>m</em> log<em>m</em>)) to the number of arcs).
 		  ///
 		  ///\tparam GR The type of the underlying digraph.
 		  ///
 		  ///\sa DynArcLookUp
 		  ///\sa AllArcLookUp
 		  template<class GR>
 		  class ArcLookUp
+		  {
 		  public:
 		    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
 		    typedef GR Digraph;
 		  protected:
 		    const Digraph &_g;
 		    typename Digraph::template NodeMap<Arc> _head;
 		    typename Digraph::template ArcMap<Arc> _left;
 		    typename Digraph::template ArcMap<Arc> _right;
 		    class ArcLess {
 		      const Digraph &g;
 		    public:
 		      ArcLess(const Digraph &_g) : g(_g) {}
 		      bool operator()(Arc a,Arc b) const
+		      {
 		        return g.target(a)<g.target(b);
+		      }
 		    };
 		  public:
 		    ///Constructor
 		    ///Constructor.
 		    ///
 		    ///It builds up the search database, which remains valid until the digraph
 		    ///changes.
 		    ArcLookUp(const Digraph &g) :_g(g),_head(g),_left(g),_right(g) {refresh();}
 		  private:
 		    Arc refreshRec(std::vector<Arc> &v,int a,int b)
+		    {
 		      int m=(a+b)/2;
 		      Arc me=v[m];
 		      _left[me] = a<m?refreshRec(v,a,m-1):INVALID;
 		      _right[me] = m<b?refreshRec(v,m+1,b):INVALID;
 		      return me;
+		    }
 		  public:
 		    ///Refresh the search data structure at a node.
 		    ///Build up the search database of node \c n.
 		    ///
 		    ///It runs in time <em>O</em>(<em>d</em> log<em>d</em>), where <em>d</em>
 		    ///is the number of the outgoing arcs of \c n.
 		    void refresh(Node n)
+		    {
 		      std::vector<Arc> v;
 		      for(OutArcIt e(_g,n);e!=INVALID;++e) v.push_back(e);
 		      if(v.size()) {
 		        std::sort(v.begin(),v.end(),ArcLess(_g));
 		        _head[n]=refreshRec(v,0,v.size()-1);
+		      }

lemon/elevator.h

Show inline comments

 		/* -*- mode: C++; indent-tabs-mode: nil; -*-
+		 *
 		 * This file is a part of LEMON, a generic C++ optimization library.
+		 *
 		 * Copyright (C) 2003-2009
 		 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
 		 * (Egervary Research Group on Combinatorial Optimization, EGRES).
+		 *
 		 * Permission to use, modify and distribute this software is granted
 		 * provided that this copyright notice appears in all copies. For
 		 * precise terms see the accompanying LICENSE file.
+		 *
 		 * This software is provided "AS IS" with no warranty of any kind,
 		 * express or implied, and with no claim as to its suitability for any
 		 * purpose.
+		 *
 		 */
 		#ifndef LEMON_ELEVATOR_H
 		#define LEMON_ELEVATOR_H
 		///\ingroup auxdat
 		///\file
 		///\brief Elevator class
 		///
 		///Elevator class implements an efficient data structure
 		///for labeling items in push-relabel type algorithms.
 		///
 		#include <lemon/core.h>
 		#include <lemon/bits/traits.h>
 		namespace lemon {
 		  ///Class for handling "labels" in push-relabel type algorithms.
 		  ///A class for handling "labels" in push-relabel type algorithms.
 		  ///
 		  ///\ingroup auxdat
 		  ///Using this class you can assign "labels" (nonnegative integer numbers)
 		  ///to the edges or nodes of a graph, manipulate and query them through
 		  ///operations typically arising in "push-relabel" type algorithms.
 		  ///
 		  ///Each item is either \em active or not, and you can also choose a
 		  ///highest level active item.
 		  ///
 		  ///\sa LinkedElevator
 		  ///
 		  ///\param GR Type of the underlying graph.
 		  ///\param Item Type of the items the data is assigned to (\c GR::Node,
 		  ///\c GR::Arc or \c GR::Edge).
 		  template<class GR, class Item>
 		  class Elevator
+		  {
 		  public:
 		    typedef Item Key;
 		    typedef int Value;
 		  private:
 		    typedef Item *Vit;
 		    typedef typename ItemSetTraits<GR,Item>::template Map<Vit>::Type VitMap;
 		    typedef typename ItemSetTraits<GR,Item>::template Map<int>::Type IntMap;
 		    const GR &_g;
 		    int _max_level;
 		    int _item_num;
 		    VitMap _where;
 		    IntMap _level;
 		    std::vector<Item> _items;
 		    std::vector<Vit> _first;
 		    std::vector<Vit> _last_active;
 		    int _highest_active;
 		    void copy(Item i, Vit p)
+		    {
-		      _where.set(*p=i,p);
+		      _where[*p=i] = p;
+		    }
 		    void copy(Vit s, Vit p)
+		    {
 		      if(s!=p)
+		        {
 		          Item i=*s;
 		          *p=i;
-		          _where.set(i,p);
+		          _where[i] = p;
+		        }
+		    }
 		    void swap(Vit i, Vit j)
+		    {
 		      Item ti=*i;
 		      Vit ct = _where[ti];
 		      _where.set(ti,_where[*i=*j]);
 		      _where.set(*j,ct);
 		      _where[ti] = _where[*i=*j];
 		      _where[*j] = ct;
 		      *j=ti;
+		    }
 		  public:
 		    ///Constructor with given maximum level.
 		    ///Constructor with given maximum level.
 		    ///
 		    ///\param graph The underlying graph.
 		    ///\param max_level The maximum allowed level.
 		    ///Set the range of the possible labels to <tt>[0..max_level]</tt>.
 		    Elevator(const GR &graph,int max_level) :
 		      _g(graph),
 		      _max_level(max_level),
 		      _item_num(_max_level),
 		      _where(graph),
 		      _level(graph,0),
 		      _items(_max_level),
 		      _first(_max_level+2),
 		      _last_active(_max_level+2),
 		      _highest_active(-1) {}
 		    ///Constructor.
 		    ///Constructor.
 		    ///
 		    ///\param graph The underlying graph.
 		    ///Set the range of the possible labels to <tt>[0..max_level]</tt>,
 		    ///where \c max_level is equal to the number of labeled items in the graph.
 		    Elevator(const GR &graph) :
 		      _g(graph),
 		      _max_level(countItems<GR, Item>(graph)),
 		      _item_num(_max_level),
 		      _where(graph),
 		      _level(graph,0),
 		      _items(_max_level),
 		      _first(_max_level+2),
 		      _last_active(_max_level+2),
 		      _highest_active(-1)
+		    {
+		    }
 		    ///Activate item \c i.
 		    ///Activate item \c i.
 		    ///\pre Item \c i shouldn't be active before.
 		    void activate(Item i)
+		    {
 		      const int l=_level[i];
 		      swap(_where[i],++_last_active[l]);
 		      if(l>_highest_active) _highest_active=l;
+		    }
 		    ///Deactivate item \c i.
 		    ///Deactivate item \c i.
 		    ///\pre Item \c i must be active before.
 		    void deactivate(Item i)
+		    {
 		      swap(_where[i],_last_active[_level[i]]--);
 		      while(_highest_active>=0 &&
 		            _last_active[_highest_active]<_first[_highest_active])
 		        _highest_active--;
+		    }
 		    ///Query whether item \c i is active
 		    bool active(Item i) const { return _where[i]<=_last_active[_level[i]]; }
 		    ///Return the level of item \c i.
 		    int operator[](Item i) const { return _level[i]; }
 		    ///Return the number of items on level \c l.
 		    int onLevel(int l) const
+		    {
 		      return _first[l+1]-_first[l];
+		    }
 		    ///Return true if level \c l is empty.
 		    bool emptyLevel(int l) const
+		    {
 		      return _first[l+1]-_first[l]==0;
+		    }
 		    ///Return the number of items above level \c l.
 		    int aboveLevel(int l) const
+		    {
 		      return _first[_max_level+1]-_first[l+1];
+		    }
 		    ///Return the number of active items on level \c l.
 		    int activesOnLevel(int l) const
+		    {
 		      return _last_active[l]-_first[l]+1;
+		    }
 		    ///Return true if there is no active item on level \c l.
 		    bool activeFree(int l) const
+		    {
 		      return _last_active[l]<_first[l];
+		    }
 		    ///Return the maximum allowed level.
 		    int maxLevel() const
+		    {
 		      return _max_level;
+		    }
 		    ///\name Highest Active Item
 		    ///Functions for working with the highest level
 		    ///active item.
 		    ///@{
 		    ///Return a highest level active item.
 		    ///Return a highest level active item or INVALID if there is no active
 		    ///item.
 		    Item highestActive() const
+		    {
 		      return _highest_active>=0?*_last_active[_highest_active]:INVALID;
+		    }
 		    ///Return the highest active level.
 		    ///Return the level of the highest active item or -1 if there is no active
 		    ///item.
 		    int highestActiveLevel() const
+		    {
 		      return _highest_active;
+		    }
 		    ///Lift the highest active item by one.
 		    ///Lift the item returned by highestActive() by one.
 		    ///
 		    void liftHighestActive()
+		    {
 		      Item it = *_last_active[_highest_active];
-		      _level.set(it,_level[it]+1);
+		      ++_level[it];
 		      swap(_last_active[_highest_active]--,_last_active[_highest_active+1]);
 		      --_first[++_highest_active];
+		    }
 		    ///Lift the highest active item to the given level.
 		    ///Lift the item returned by highestActive() to level \c new_level.
 		    ///
 		    ///\warning \c new_level must be strictly higher
 		    ///than the current level.
 		    ///
 		    void liftHighestActive(int new_level)
+		    {
 		      const Item li = *_last_active[_highest_active];
 		      copy(--_first[_highest_active+1],_last_active[_highest_active]--);
 		      for(int l=_highest_active+1;l<new_level;l++)
+		        {
 		          copy(--_first[l+1],_first[l]);
 		          --_last_active[l];
+		        }
 		      copy(li,_first[new_level]);
-		      _level.set(li,new_level);
+		      _level[li] = new_level;
 		      _highest_active=new_level;
+		    }
 		    ///Lift the highest active item to the top level.
 		    ///Lift the item returned by highestActive() to the top level and
 		    ///deactivate it.
 		    void liftHighestActiveToTop()
+		    {
 		      const Item li = *_last_active[_highest_active];
 		      copy(--_first[_highest_active+1],_last_active[_highest_active]--);
 		      for(int l=_highest_active+1;l<_max_level;l++)
+		        {
 		          copy(--_first[l+1],_first[l]);
 		          --_last_active[l];
+		        }
 		      copy(li,_first[_max_level]);
 		      --_last_active[_max_level];
-		      _level.set(li,_max_level);
+		      _level[li] = _max_level;
 		      while(_highest_active>=0 &&
 		            _last_active[_highest_active]<_first[_highest_active])
 		        _highest_active--;
+		    }
 		    ///@}
 		    ///\name Active Item on Certain Level
 		    ///Functions for working with the active items.
 		    ///@{
 		    ///Return an active item on level \c l.
 		    ///Return an active item on level \c l or \ref INVALID if there is no such
 		    ///an item. (\c l must be from the range [0...\c max_level].
 		    Item activeOn(int l) const
+		    {
 		      return _last_active[l]>=_first[l]?*_last_active[l]:INVALID;
+		    }
 		    ///Lift the active item returned by \c activeOn(level) by one.
 		    ///Lift the active item returned by \ref activeOn() "activeOn(level)"
 		    ///by one.
 		    Item liftActiveOn(int level)
+		    {
 		      Item it =*_last_active[level];
-		      _level.set(it,_level[it]+1);
+		      ++_level[it];
 		      swap(_last_active[level]--, --_first[level+1]);
 		      if (level+1>_highest_active) ++_highest_active;
+		    }
 		    ///Lift the active item returned by \c activeOn(level) to the given level.
 		    ///Lift the active item returned by \ref activeOn() "activeOn(level)"
 		    ///to the given level.
 		    void liftActiveOn(int level, int new_level)
+		    {
 		      const Item ai = *_last_active[level];
 		      copy(--_first[level+1], _last_active[level]--);
 		      for(int l=level+1;l<new_level;l++)
+		        {
 		          copy(_last_active[l],_first[l]);
 		          copy(--_first[l+1], _last_active[l]--);
+		        }
 		      copy(ai,_first[new_level]);
-		      _level.set(ai,new_level);
+		      _level[ai] = new_level;
 		      if (new_level>_highest_active) _highest_active=new_level;
+		    }
 		    ///Lift the active item returned by \c activeOn(level) to the top level.
 		    ///Lift the active item returned by \ref activeOn() "activeOn(level)"
 		    ///to the top level and deactivate it.
 		    void liftActiveToTop(int level)
+		    {
 		      const Item ai = *_last_active[level];
 		      copy(--_first[level+1],_last_active[level]--);
 		      for(int l=level+1;l<_max_level;l++)
+		        {
 		          copy(_last_active[l],_first[l]);
 		          copy(--_first[l+1], _last_active[l]--);
+		        }
 		      copy(ai,_first[_max_level]);
 		      --_last_active[_max_level];
-		      _level.set(ai,_max_level);
+		      _level[ai] = _max_level;
 		      if (_highest_active==level) {
 		        while(_highest_active>=0 &&
 		              _last_active[_highest_active]<_first[_highest_active])
 		          _highest_active--;
+		      }
+		    }
 		    ///@}
 		    ///Lift an active item to a higher level.
 		    ///Lift an active item to a higher level.
 		    ///\param i The item to be lifted. It must be active.
 		    ///\param new_level The new level of \c i. It must be strictly higher
 		    ///than the current level.
 		    ///
 		    void lift(Item i, int new_level)
+		    {
 		      const int lo = _level[i];
 		      const Vit w = _where[i];
 		      copy(_last_active[lo],w);
 		      copy(--_first[lo+1],_last_active[lo]--);
 		      for(int l=lo+1;l<new_level;l++)
+		        {
 		          copy(_last_active[l],_first[l]);
 		          copy(--_first[l+1],_last_active[l]--);
+		        }
 		      copy(i,_first[new_level]);
-		      _level.set(i,new_level);
+		      _level[i] = new_level;
 		      if(new_level>_highest_active) _highest_active=new_level;
+		    }
 		    ///Move an inactive item to the top but one level (in a dirty way).
 		    ///This function moves an inactive item from the top level to the top
 		    ///but one level (in a dirty way).
 		    ///\warning It makes the underlying datastructure corrupt, so use it
 		    ///only if you really know what it is for.
 		    ///\pre The item is on the top level.
 		    void dirtyTopButOne(Item i) {
-		      _level.set(i,_max_level - 1);
+		      _level[i] = _max_level - 1;
+		    }
 		    ///Lift all items on and above the given level to the top level.
 		    ///This function lifts all items on and above level \c l to the top
 		    ///level and deactivates them.
 		    void liftToTop(int l)
+		    {
 		      const Vit f=_first[l];
 		      const Vit tl=_first[_max_level];
 		      for(Vit i=f;i!=tl;++i)
-		        _level.set(*i,_max_level);
+		        _level[*i] = _max_level;
 		      for(int i=l;i<=_max_level;i++)
+		        {
 		          _first[i]=f;
 		          _last_active[i]=f-1;
+		        }
 		      for(_highest_active=l-1;
 		          _highest_active>=0 &&
 		            _last_active[_highest_active]<_first[_highest_active];
 		          _highest_active--) ;
+		    }
 		  private:
 		    int _init_lev;
 		    Vit _init_num;
 		  public:
 		    ///\name Initialization
 		    ///Using these functions you can initialize the levels of the items.
 		    ///\n
 		    ///The initialization must be started with calling \c initStart().
 		    ///Then the items should be listed level by level starting with the
 		    ///lowest one (level 0) using \c initAddItem() and \c initNewLevel().
 		    ///Finally \c initFinish() must be called.
 		    ///The items not listed are put on the highest level.
 		    ///@{
 		    ///Start the initialization process.
 		    void initStart()
+		    {
 		      _init_lev=0;
 		      _init_num=&_items[0];
 		      _first[0]=&_items[0];
 		      _last_active[0]=&_items[0]-1;
 		      Vit n=&_items[0];
 		      for(typename ItemSetTraits<GR,Item>::ItemIt i(_g);i!=INVALID;++i)
+		        {
 		          *n=i;
 		          _where.set(i,n);
 		          _level.set(i,_max_level);
 		          _where[i] = n;
 		          _level[i] = _max_level;
 		          ++n;
+		        }
+		    }
 		    ///Add an item to the current level.
 		    void initAddItem(Item i)
+		    {
 		      swap(_where[i],_init_num);
-		      _level.set(i,_init_lev);
+		      _level[i] = _init_lev;
 		      ++_init_num;
+		    }
 		    ///Start a new level.
 		    ///Start a new level.
 		    ///It shouldn't be used before the items on level 0 are listed.
 		    void initNewLevel()
+		    {
 		      _init_lev++;
 		      _first[_init_lev]=_init_num;
 		      _last_active[_init_lev]=_init_num-1;
+		    }
 		    ///Finalize the initialization process.
 		    void initFinish()
+		    {
 		      for(_init_lev++;_init_lev<=_max_level;_init_lev++)
+		        {
 		          _first[_init_lev]=_init_num;
 		          _last_active[_init_lev]=_init_num-1;
+		        }
 		      _first[_max_level+1]=&_items[0]+_item_num;
 		      _last_active[_max_level+1]=&_items[0]+_item_num-1;
 		      _highest_active = -1;
+		    }
 		    ///@}
 		  };
 		  ///Class for handling "labels" in push-relabel type algorithms.
 		  ///A class for handling "labels" in push-relabel type algorithms.
 		  ///
 		  ///\ingroup auxdat
 		  ///Using this class you can assign "labels" (nonnegative integer numbers)
 		  ///to the edges or nodes of a graph, manipulate and query them through
 		  ///operations typically arising in "push-relabel" type algorithms.
 		  ///
 		  ///Each item is either \em active or not, and you can also choose a
 		  ///highest level active item.
 		  ///
 		  ///\sa Elevator
 		  ///
 		  ///\param GR Type of the underlying graph.
 		  ///\param Item Type of the items the data is assigned to (\c GR::Node,
 		  ///\c GR::Arc or \c GR::Edge).
 		  template <class GR, class Item>
 		  class LinkedElevator {
 		  public:
 		    typedef Item Key;
 		    typedef int Value;
 		  private:
 		    typedef typename ItemSetTraits<GR,Item>::
 		    template Map<Item>::Type ItemMap;
 		    typedef typename ItemSetTraits<GR,Item>::
 		    template Map<int>::Type IntMap;
 		    typedef typename ItemSetTraits<GR,Item>::
 		    template Map<bool>::Type BoolMap;
 		    const GR &_graph;
 		    int _max_level;
 		    int _item_num;
 		    std::vector<Item> _first, _last;
 		    ItemMap _prev, _next;
 		    int _highest_active;
 		    IntMap _level;
 		    BoolMap _active;
 		  public:
 		    ///Constructor with given maximum level.
 		    ///Constructor with given maximum level.
 		    ///
 		    ///\param graph The underlying graph.
 		    ///\param max_level The maximum allowed level.
 		    ///Set the range of the possible labels to <tt>[0..max_level]</tt>.
 		    LinkedElevator(const GR& graph, int max_level)
 		      : _graph(graph), _max_level(max_level), _item_num(_max_level),
 		        _first(_max_level + 1), _last(_max_level + 1),
 		        _prev(graph), _next(graph),
 		        _highest_active(-1), _level(graph), _active(graph) {}
 		    ///Constructor.
 		    ///Constructor.
 		    ///
 		    ///\param graph The underlying graph.
 		    ///Set the range of the possible labels to <tt>[0..max_level]</tt>,
 		    ///where \c max_level is equal to the number of labeled items in the graph.
 		    LinkedElevator(const GR& graph)
 		      : _graph(graph), _max_level(countItems<GR, Item>(graph)),
 		        _item_num(_max_level),
 		        _first(_max_level + 1), _last(_max_level + 1),
 		        _prev(graph, INVALID), _next(graph, INVALID),
 		        _highest_active(-1), _level(graph), _active(graph) {}
 		    ///Activate item \c i.
 		    ///Activate item \c i.
 		    ///\pre Item \c i shouldn't be active before.
 		    void activate(Item i) {
-		      _active.set(i, true);
+		      _active[i] = true;
 		      int level = _level[i];
 		      if (level > _highest_active) {
 		        _highest_active = level;
+		      }
 		      if (_prev[i] == INVALID || _active[_prev[i]]) return;
 		      //unlace
-		      _next.set(_prev[i], _next[i]);
+		      _next[_prev[i]] = _next[i];
 		      if (_next[i] != INVALID) {
-		        _prev.set(_next[i], _prev[i]);
+		        _prev[_next[i]] = _prev[i];
 		      } else {
 		        _last[level] = _prev[i];
+		      }
 		      //lace
 		      _next.set(i, _first[level]);
 		      _prev.set(_first[level], i);
-		      _prev.set(i, INVALID);
+		      _next[i] = _first[level];
 		      _prev[_first[level]] = i;
 		      _prev[i] = INVALID;
 		      _first[level] = i;
+		    }
 		    ///Deactivate item \c i.
 		    ///Deactivate item \c i.
 		    ///\pre Item \c i must be active before.
 		    void deactivate(Item i) {
-		      _active.set(i, false);
+		      _active[i] = false;
 		      int level = _level[i];
 		      if (_next[i] == INVALID || !_active[_next[i]])
 		        goto find_highest_level;
 		      //unlace
-		      _prev.set(_next[i], _prev[i]);
+		      _prev[_next[i]] = _prev[i];
 		      if (_prev[i] != INVALID) {
-		        _next.set(_prev[i], _next[i]);
+		        _next[_prev[i]] = _next[i];
 		      } else {
 		        _first[_level[i]] = _next[i];
+		      }
 		      //lace
 		      _prev.set(i, _last[level]);
 		      _next.set(_last[level], i);
-		      _next.set(i, INVALID);
+		      _prev[i] = _last[level];
 		      _next[_last[level]] = i;
 		      _next[i] = INVALID;
 		      _last[level] = i;
 		    find_highest_level:
 		      if (level == _highest_active) {
 		        while (_highest_active >= 0 && activeFree(_highest_active))
 		          --_highest_active;
+		      }
+		    }
 		    ///Query whether item \c i is active
 		    bool active(Item i) const { return _active[i]; }
 		    ///Return the level of item \c i.
 		    int operator[](Item i) const { return _level[i]; }
 		    ///Return the number of items on level \c l.
 		    int onLevel(int l) const {
 		      int num = 0;
 		      Item n = _first[l];
 		      while (n != INVALID) {
 		        ++num;
 		        n = _next[n];
+		      }
 		      return num;
+		    }
 		    ///Return true if the level is empty.
 		    bool emptyLevel(int l) const {
 		      return _first[l] == INVALID;
+		    }
 		    ///Return the number of items above level \c l.
 		    int aboveLevel(int l) const {
 		      int num = 0;
 		      for (int level = l + 1; level < _max_level; ++level)
 		        num += onLevel(level);
 		      return num;
+		    }
 		    ///Return the number of active items on level \c l.
 		    int activesOnLevel(int l) const {
 		      int num = 0;
 		      Item n = _first[l];
 		      while (n != INVALID && _active[n]) {
 		        ++num;
 		        n = _next[n];
+		      }
 		      return num;
+		    }
 		    ///Return true if there is no active item on level \c l.
 		    bool activeFree(int l) const {
 		      return _first[l] == INVALID || !_active[_first[l]];
+		    }
 		    ///Return the maximum allowed level.
 		    int maxLevel() const {
 		      return _max_level;
+		    }
 		    ///\name Highest Active Item
 		    ///Functions for working with the highest level
 		    ///active item.
 		    ///@{
 		    ///Return a highest level active item.
 		    ///Return a highest level active item or INVALID if there is no active
 		    ///item.
 		    Item highestActive() const {
 		      return _highest_active >= 0 ? _first[_highest_active] : INVALID;
+		    }
 		    ///Return the highest active level.
 		    ///Return the level of the highest active item or -1 if there is no active
 		    ///item.
 		    int highestActiveLevel() const {
 		      return _highest_active;
+		    }
 		    ///Lift the highest active item by one.
 		    ///Lift the item returned by highestActive() by one.
 		    ///
 		    void liftHighestActive() {
 		      Item i = _first[_highest_active];
 		      if (_next[i] != INVALID) {
-		        _prev.set(_next[i], INVALID);
+		        _prev[_next[i]] = INVALID;
 		        _first[_highest_active] = _next[i];
 		      } else {
 		        _first[_highest_active] = INVALID;
 		        _last[_highest_active] = INVALID;
+		      }
-		      _level.set(i, ++_highest_active);
+		      _level[i] = ++_highest_active;
 		      if (_first[_highest_active] == INVALID) {
 		        _first[_highest_active] = i;
 		        _last[_highest_active] = i;
 		        _prev.set(i, INVALID);
 		        _next.set(i, INVALID);
 		        _prev[i] = INVALID;
 		        _next[i] = INVALID;
 		      } else {
 		        _prev.set(_first[_highest_active], i);
 		        _next.set(i, _first[_highest_active]);
 		        _prev[_first[_highest_active]] = i;
 		        _next[i] = _first[_highest_active];
 		        _first[_highest_active] = i;
+		      }
+		    }
 		    ///Lift the highest active item to the given level.
 		    ///Lift the item returned by highestActive() to level \c new_level.
 		    ///
 		    ///\warning \c new_level must be strictly higher
 		    ///than the current level.
 		    ///
 		    void liftHighestActive(int new_level) {
 		      Item i = _first[_highest_active];
 		      if (_next[i] != INVALID) {
-		        _prev.set(_next[i], INVALID);
+		        _prev[_next[i]] = INVALID;
 		        _first[_highest_active] = _next[i];
 		      } else {
 		        _first[_highest_active] = INVALID;
 		        _last[_highest_active] = INVALID;
+		      }
-		      _level.set(i, _highest_active = new_level);
+		      _level[i] = _highest_active = new_level;
 		      if (_first[_highest_active] == INVALID) {
 		        _first[_highest_active] = _last[_highest_active] = i;
 		        _prev.set(i, INVALID);
 		        _next.set(i, INVALID);
 		        _prev[i] = INVALID;
 		        _next[i] = INVALID;
 		      } else {
 		        _prev.set(_first[_highest_active], i);
 		        _next.set(i, _first[_highest_active]);
 		        _prev[_first[_highest_active]] = i;
 		        _next[i] = _first[_highest_active];
 		        _first[_highest_active] = i;
+		      }
+		    }
 		    ///Lift the highest active item to the top level.
 		    ///Lift the item returned by highestActive() to the top level and
 		    ///deactivate it.
 		    void liftHighestActiveToTop() {
 		      Item i = _first[_highest_active];
-		      _level.set(i, _max_level);
+		      _level[i] = _max_level;
 		      if (_next[i] != INVALID) {
-		        _prev.set(_next[i], INVALID);
+		        _prev[_next[i]] = INVALID;
 		        _first[_highest_active] = _next[i];
 		      } else {
 		        _first[_highest_active] = INVALID;
 		        _last[_highest_active] = INVALID;
+		      }
 		      while (_highest_active >= 0 && activeFree(_highest_active))
 		        --_highest_active;
+		    }
 		    ///@}
 		    ///\name Active Item on Certain Level
 		    ///Functions for working with the active items.
 		    ///@{
 		    ///Return an active item on level \c l.
 		    ///Return an active item on level \c l or \ref INVALID if there is no such
 		    ///an item. (\c l must be from the range [0...\c max_level].
 		    Item activeOn(int l) const
+		    {
 		      return _active[_first[l]] ? _first[l] : INVALID;
+		    }
 		    ///Lift the active item returned by \c activeOn(l) by one.
 		    ///Lift the active item returned by \ref activeOn() "activeOn(l)"
 		    ///by one.
 		    Item liftActiveOn(int l)
+		    {
 		      Item i = _first[l];
 		      if (_next[i] != INVALID) {
-		        _prev.set(_next[i], INVALID);
+		        _prev[_next[i]] = INVALID;
 		        _first[l] = _next[i];
 		      } else {
 		        _first[l] = INVALID;
 		        _last[l] = INVALID;
+		      }
-		      _level.set(i, ++l);
+		      _level[i] = ++l;
 		      if (_first[l] == INVALID) {
 		        _first[l] = _last[l] = i;
 		        _prev.set(i, INVALID);
 		        _next.set(i, INVALID);
 		        _prev[i] = INVALID;
 		        _next[i] = INVALID;
 		      } else {
 		        _prev.set(_first[l], i);
 		        _next.set(i, _first[l]);
 		        _prev[_first[l]] = i;
 		        _next[i] = _first[l];
 		        _first[l] = i;
+		      }
 		      if (_highest_active < l) {
 		        _highest_active = l;
+		      }
+		    }
 		    ///Lift the active item returned by \c activeOn(l) to the given level.
 		    ///Lift the active item returned by \ref activeOn() "activeOn(l)"
 		    ///to the given level.
 		    void liftActiveOn(int l, int new_level)
+		    {
 		      Item i = _first[l];
 		      if (_next[i] != INVALID) {
-		        _prev.set(_next[i], INVALID);
+		        _prev[_next[i]] = INVALID;
 		        _first[l] = _next[i];
 		      } else {
 		        _first[l] = INVALID;
 		        _last[l] = INVALID;
+		      }
-		      _level.set(i, l = new_level);
+		      _level[i] = l = new_level;
 		      if (_first[l] == INVALID) {
 		        _first[l] = _last[l] = i;
 		        _prev.set(i, INVALID);
 		        _next.set(i, INVALID);
 		        _prev[i] = INVALID;
 		        _next[i] = INVALID;
 		      } else {
 		        _prev.set(_first[l], i);
 		        _next.set(i, _first[l]);
 		        _prev[_first[l]] = i;
 		        _next[i] = _first[l];
 		        _first[l] = i;
+		      }
 		      if (_highest_active < l) {
 		        _highest_active = l;
+		      }
+		    }
 		    ///Lift the active item returned by \c activeOn(l) to the top level.
 		    ///Lift the active item returned by \ref activeOn() "activeOn(l)"
 		    ///to the top level and deactivate it.
 		    void liftActiveToTop(int l)
+		    {
 		      Item i = _first[l];
 		      if (_next[i] != INVALID) {
-		        _prev.set(_next[i], INVALID);
+		        _prev[_next[i]] = INVALID;
 		        _first[l] = _next[i];
 		      } else {
 		        _first[l] = INVALID;
 		        _last[l] = INVALID;
+		      }
-		      _level.set(i, _max_level);
+		      _level[i] = _max_level;
 		      if (l == _highest_active) {
 		        while (_highest_active >= 0 && activeFree(_highest_active))
 		          --_highest_active;
+		      }
+		    }
 		    ///@}
 		    /// \brief Lift an active item to a higher level.
 		    ///
 		    /// Lift an active item to a higher level.
 		    /// \param i The item to be lifted. It must be active.
 		    /// \param new_level The new level of \c i. It must be strictly higher
 		    /// than the current level.
 		    ///
 		    void lift(Item i, int new_level) {
 		      if (_next[i] != INVALID) {
-		        _prev.set(_next[i], _prev[i]);
+		        _prev[_next[i]] = _prev[i];
 		      } else {
 		        _last[new_level] = _prev[i];
+		      }
 		      if (_prev[i] != INVALID) {
-		        _next.set(_prev[i], _next[i]);
+		        _next[_prev[i]] = _next[i];
 		      } else {
 		        _first[new_level] = _next[i];
+		      }
-		      _level.set(i, new_level);
+		      _level[i] = new_level;
 		      if (_first[new_level] == INVALID) {
 		        _first[new_level] = _last[new_level] = i;
 		        _prev.set(i, INVALID);
 		        _next.set(i, INVALID);
 		        _prev[i] = INVALID;
 		        _next[i] = INVALID;
 		      } else {
 		        _prev.set(_first[new_level], i);
 		        _next.set(i, _first[new_level]);
 		        _prev[_first[new_level]] = i;
 		        _next[i] = _first[new_level];
 		        _first[new_level] = i;
+		      }
 		      if (_highest_active < new_level) {
 		        _highest_active = new_level;
+		      }
+		    }
 		    ///Move an inactive item to the top but one level (in a dirty way).
 		    ///This function moves an inactive item from the top level to the top
 		    ///but one level (in a dirty way).
 		    ///\warning It makes the underlying datastructure corrupt, so use it
 		    ///only if you really know what it is for.
 		    ///\pre The item is on the top level.
 		    void dirtyTopButOne(Item i) {
-		      _level.set(i, _max_level - 1);
+		      _level[i] = _max_level - 1;
+		    }
 		    ///Lift all items on and above the given level to the top level.
 		    ///This function lifts all items on and above level \c l to the top
 		    ///level and deactivates them.
 		    void liftToTop(int l)  {
 		      for (int i = l + 1; _first[i] != INVALID; ++i) {
 		        Item n = _first[i];
 		        while (n != INVALID) {
-		          _level.set(n, _max_level);
+		          _level[n] = _max_level;
 		          n = _next[n];
+		        }
 		        _first[i] = INVALID;
 		        _last[i] = INVALID;
+		      }
 		      if (_highest_active > l - 1) {
 		        _highest_active = l - 1;
 		        while (_highest_active >= 0 && activeFree(_highest_active))
 		          --_highest_active;
+		      }
+		    }
 		  private:
 		    int _init_level;
 		  public:
 		    ///\name Initialization
 		    ///Using these functions you can initialize the levels of the items.
 		    ///\n
 		    ///The initialization must be started with calling \c initStart().
 		    ///Then the items should be listed level by level starting with the
 		    ///lowest one (level 0) using \c initAddItem() and \c initNewLevel().
 		    ///Finally \c initFinish() must be called.
 		    ///The items not listed are put on the highest level.
 		    ///@{
 		    ///Start the initialization process.
 		    void initStart() {
 		      for (int i = 0; i <= _max_level; ++i) {
 		        _first[i] = _last[i] = INVALID;
+		      }
 		      _init_level = 0;
 		      for(typename ItemSetTraits<GR,Item>::ItemIt i(_graph);
 		          i != INVALID; ++i) {
 		        _level.set(i, _max_level);
 		        _active.set(i, false);
 		        _level[i] = _max_level;
 		        _active[i] = false;
+		      }
+		    }
 		    ///Add an item to the current level.
 		    void initAddItem(Item i) {
-		      _level.set(i, _init_level);
+		      _level[i] = _init_level;
 		      if (_last[_init_level] == INVALID) {
 		        _first[_init_level] = i;
 		        _last[_init_level] = i;
 		        _prev.set(i, INVALID);
 		        _next.set(i, INVALID);
 		        _prev[i] = INVALID;
 		        _next[i] = INVALID;
 		      } else {
 		        _prev.set(i, _last[_init_level]);
 		        _next.set(i, INVALID);
-		        _next.set(_last[_init_level], i);
+		        _prev[i] = _last[_init_level];
 		        _next[i] = INVALID;
 		        _next[_last[_init_level]] = i;
 		        _last[_init_level] = i;
+		      }
+		    }
 		    ///Start a new level.
 		    ///Start a new level.
 		    ///It shouldn't be used before the items on level 0 are listed.
 		    void initNewLevel() {
 		      ++_init_level;
+		    }
 		    ///Finalize the initialization process.
 		    void initFinish() {
 		      _highest_active = -1;
+		    }
 		    ///@}
 		  };
 		} //END OF NAMESPACE LEMON
 		#endif

lemon/gomory_hu.h

Show inline comments

 		/* -*- C++ -*-
+		 *
 		 * This file is a part of LEMON, a generic C++ optimization library
+		 *
 		 * Copyright (C) 2003-2008
 		 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
 		 * (Egervary Research Group on Combinatorial Optimization, EGRES).
+		 *
 		 * Permission to use, modify and distribute this software is granted
 		 * provided that this copyright notice appears in all copies. For
 		 * precise terms see the accompanying LICENSE file.
+		 *
 		 * This software is provided "AS IS" with no warranty of any kind,
 		 * express or implied, and with no claim as to its suitability for any
 		 * purpose.
+		 *
 		 */
 		#ifndef LEMON_GOMORY_HU_TREE_H
 		#define LEMON_GOMORY_HU_TREE_H
 		#include <limits>
 		#include <lemon/core.h>
 		#include <lemon/preflow.h>
 		#include <lemon/concept_check.h>
 		#include <lemon/concepts/maps.h>
 		/// \ingroup min_cut
 		/// \file
 		/// \brief Gomory-Hu cut tree in graphs.
 		namespace lemon {
 		  /// \ingroup min_cut
 		  ///
 		  /// \brief Gomory-Hu cut tree algorithm
 		  ///
 		  /// The Gomory-Hu tree is a tree on the node set of a given graph, but it
 		  /// may contain edges which are not in the original graph. It has the
 		  /// property that the minimum capacity edge of the path between two nodes
 		  /// in this tree has the same weight as the minimum cut in the graph
 		  /// between these nodes. Moreover the components obtained by removing
 		  /// this edge from the tree determine the corresponding minimum cut.
 		  ///
 		  /// Therefore once this tree is computed, the minimum cut between any pair
 		  /// of nodes can easily be obtained.
 		  ///
 		  /// The algorithm calculates \e n-1 distinct minimum cuts (currently with
 		  /// the \ref Preflow algorithm), therefore the algorithm has
 		  /// \f$(O(n^3\sqrt{e})\f$ overall time complexity. It calculates a
 		  /// rooted Gomory-Hu tree, its structure and the weights can be obtained
 		  /// by \c predNode(), \c predValue() and \c rootDist().
 		  ///
 		  /// The members \c minCutMap() and \c minCutValue() calculate
 		  /// the minimum cut and the minimum cut value between any two nodes
 		  /// in the graph. You can also list (iterate on) the nodes and the
 		  /// edges of the cuts using \c MinCutNodeIt and \c MinCutEdgeIt.
 		  ///
 		  /// \tparam GR The type of the undirected graph the algorithm runs on.
 		  /// \tparam CAP The type of the edge map describing the edge capacities.
 		  /// It is \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>" by default.
 		#ifdef DOXYGEN
 		  template <typename GR,
 			    typename CAP>
 		#else
 		  template <typename GR,
 			    typename CAP = typename GR::template EdgeMap<int> >
 		#endif
 		  class GomoryHu {
 		  public:
 		    /// The graph type
 		    typedef GR Graph;
 		    /// The type of the edge capacity map
 		    typedef CAP Capacity;
 		    /// The value type of capacities
 		    typedef typename Capacity::Value Value;
 		  private:
 		    TEMPLATE_GRAPH_TYPEDEFS(Graph);
 		    const Graph& _graph;
 		    const Capacity& _capacity;
 		    Node _root;
 		    typename Graph::template NodeMap<Node>* _pred;
 		    typename Graph::template NodeMap<Value>* _weight;
 		    typename Graph::template NodeMap<int>* _order;
 		    void createStructures() {
 		      if (!_pred) {
 			_pred = new typename Graph::template NodeMap<Node>(_graph);
+		      }
 		      if (!_weight) {
 			_weight = new typename Graph::template NodeMap<Value>(_graph);
+		      }
 		      if (!_order) {
 			_order = new typename Graph::template NodeMap<int>(_graph);
+		      }
+		    }
 		    void destroyStructures() {
 		      if (_pred) {
 			delete _pred;
+		      }
 		      if (_weight) {
 			delete _weight;
+		      }
 		      if (_order) {
 			delete _order;
+		      }
+		    }
 		  public:
 		    /// \brief Constructor
 		    ///
 		    /// Constructor
 		    /// \param graph The undirected graph the algorithm runs on.
 		    /// \param capacity The edge capacity map.
 		    GomoryHu(const Graph& graph, const Capacity& capacity)
 		      : _graph(graph), _capacity(capacity),
 			_pred(0), _weight(0), _order(0)
+		    {
 		      checkConcept<concepts::ReadMap<Edge, Value>, Capacity>();
+		    }
 		    /// \brief Destructor
 		    ///
 		    /// Destructor
 		    ~GomoryHu() {
 		      destroyStructures();
+		    }
 		  private:
 		    // Initialize the internal data structures
 		    void init() {
 		      createStructures();
 		      _root = NodeIt(_graph);
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 			_pred->set(n, _root);
 			_order->set(n, -1);
 		        (*_pred)[n] = _root;
 		        (*_order)[n] = -1;
+		      }
 		      _pred->set(_root, INVALID);
 		      _weight->set(_root, std::numeric_limits<Value>::max());
 		      (*_pred)[_root] = INVALID;
 		      (*_weight)[_root] = std::numeric_limits<Value>::max();
+		    }
 		    // Start the algorithm
 		    void start() {
 		      Preflow<Graph, Capacity> fa(_graph, _capacity, _root, INVALID);
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 			if (n == _root) continue;
 			Node pn = (*_pred)[n];
 			fa.source(n);
 			fa.target(pn);
 			fa.runMinCut();
-			_weight->set(n, fa.flowValue());
+			(*_weight)[n] = fa.flowValue();
 			for (NodeIt nn(_graph); nn != INVALID; ++nn) {
 			  if (nn != n && fa.minCut(nn) && (*_pred)[nn] == pn) {
-			    _pred->set(nn, n);
+			    (*_pred)[nn] = n;
+			  }
+			}
 			if ((*_pred)[pn] != INVALID && fa.minCut((*_pred)[pn])) {
 			  _pred->set(n, (*_pred)[pn]);
 			  _pred->set(pn, n);
 			  _weight->set(n, (*_weight)[pn]);
 			  _weight->set(pn, fa.flowValue());
 			  (*_pred)[n] = (*_pred)[pn];
 			  (*_pred)[pn] = n;
 			  (*_weight)[n] = (*_weight)[pn];
 			  (*_weight)[pn] = fa.flowValue();
+			}
+		      }
-		      _order->set(_root, 0);
+		      (*_order)[_root] = 0;
 		      int index = 1;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 			std::vector<Node> st;
 			Node nn = n;
 			while ((*_order)[nn] == -1) {
 			  st.push_back(nn);
 			  nn = (*_pred)[nn];
+			}
 			while (!st.empty()) {
-			  _order->set(st.back(), index++);
+			  (*_order)[st.back()] = index++;
 			  st.pop_back();
+			}
+		      }
+		    }
 		  public:
 		    ///\name Execution Control
 		    ///@{
 		    /// \brief Run the Gomory-Hu algorithm.
 		    ///
 		    /// This function runs the Gomory-Hu algorithm.
 		    void run() {
 		      init();
 		      start();
+		    }
 		    /// @}
 		    ///\name Query Functions
 		    ///The results of the algorithm can be obtained using these
 		    ///functions.\n
 		    ///\ref run() "run()" should be called before using them.\n
 		    ///See also \ref MinCutNodeIt and \ref MinCutEdgeIt.
 		    ///@{
 		    /// \brief Return the predecessor node in the Gomory-Hu tree.
 		    ///
 		    /// This function returns the predecessor node in the Gomory-Hu tree.
 		    /// If the node is
 		    /// the root of the Gomory-Hu tree, then it returns \c INVALID.
 		    Node predNode(const Node& node) {
 		      return (*_pred)[node];
+		    }
 		    /// \brief Return the distance from the root node in the Gomory-Hu tree.
 		    ///
 		    /// This function returns the distance of \c node from the root node
 		    /// in the Gomory-Hu tree.
 		    int rootDist(const Node& node) {
 		      return (*_order)[node];
+		    }
 		    /// \brief Return the weight of the predecessor edge in the
 		    /// Gomory-Hu tree.
 		    ///
 		    /// This function returns the weight of the predecessor edge in the
 		    /// Gomory-Hu tree.  If the node is the root, the result is undefined.
 		    Value predValue(const Node& node) {
 		      return (*_weight)[node];
+		    }
 		    /// \brief Return the minimum cut value between two nodes
 		    ///
 		    /// This function returns the minimum cut value between two nodes. The
 		    /// algorithm finds the nearest common ancestor in the Gomory-Hu
 		    /// tree and calculates the minimum weight edge on the paths to
 		    /// the ancestor.
 		    Value minCutValue(const Node& s, const Node& t) const {
 		      Node sn = s, tn = t;
 		      Value value = std::numeric_limits<Value>::max();
 		      while (sn != tn) {
 			if ((*_order)[sn] < (*_order)[tn]) {
 			  if ((*_weight)[tn] <= value) value = (*_weight)[tn];
 			  tn = (*_pred)[tn];
 			} else {
 			  if ((*_weight)[sn] <= value) value = (*_weight)[sn];
 			  sn = (*_pred)[sn];
+			}
+		      }
 		      return value;
+		    }
 		    /// \brief Return the minimum cut between two nodes
 		    ///
 		    /// This function returns the minimum cut between the nodes \c s and \c t
 		    /// in the \c cutMap parameter by setting the nodes in the component of
 		    /// \c s to \c true and the other nodes to \c false.
 		    ///
 		    /// For higher level interfaces, see MinCutNodeIt and MinCutEdgeIt.
 		    template <typename CutMap>
 		    Value minCutMap(const Node& s, ///< The base node.
 		                    const Node& t,
 		                    ///< The node you want to separate from node \c s.
 		                    CutMap& cutMap
 		                    ///< The cut will be returned in this map.
 		                    /// It must be a \c bool (or convertible)
 		                    /// \ref concepts::ReadWriteMap "ReadWriteMap"
 		                    /// on the graph nodes.
 		                    ) const {
 		      Node sn = s, tn = t;
 		      bool s_root=false;
 		      Node rn = INVALID;
 		      Value value = std::numeric_limits<Value>::max();
 		      while (sn != tn) {
 			if ((*_order)[sn] < (*_order)[tn]) {
 			  if ((*_weight)[tn] <= value) {
 			    rn = tn;
 		            s_root = false;
 			    value = (*_weight)[tn];
+			  }
 			  tn = (*_pred)[tn];
 			} else {
 			  if ((*_weight)[sn] <= value) {
 			    rn = sn;
 		            s_root = true;
 			    value = (*_weight)[sn];
+			  }
 			  sn = (*_pred)[sn];
+			}
+		      }
 		      typename Graph::template NodeMap<bool> reached(_graph, false);
-		      reached.set(_root, true);
+		      reached[_root] = true;
 		      cutMap.set(_root, !s_root);
-		      reached.set(rn, true);
+		      reached[rn] = true;
 		      cutMap.set(rn, s_root);
 		      std::vector<Node> st;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 			st.clear();
 		        Node nn = n;
 			while (!reached[nn]) {
 			  st.push_back(nn);
 			  nn = (*_pred)[nn];
+			}
 			while (!st.empty()) {
 			  cutMap.set(st.back(), cutMap[nn]);
 			  st.pop_back();
+			}
+		      }
 		      return value;
+		    }
 		    ///@}
 		    friend class MinCutNodeIt;
 		    /// Iterate on the nodes of a minimum cut
 		    /// This iterator class lists the nodes of a minimum cut found by
 		    /// GomoryHu. Before using it, you must allocate a GomoryHu class,
 		    /// and call its \ref GomoryHu::run() "run()" method.
 		    ///
 		    /// This example counts the nodes in the minimum cut separating \c s from
 		    /// \c t.
 		    /// \code
 		    /// GomoruHu<Graph> gom(g, capacities);
 		    /// gom.run();
 		    /// int cnt=0;
 		    /// for(GomoruHu<Graph>::MinCutNodeIt n(gom,s,t); n!=INVALID; ++n) ++cnt;
 		    /// \endcode
 		    class MinCutNodeIt
+		    {
 		      bool _side;
 		      typename Graph::NodeIt _node_it;
 		      typename Graph::template NodeMap<bool> _cut;
 		    public:
 		      /// Constructor
 		      /// Constructor.
 		      ///
 		      MinCutNodeIt(GomoryHu const &gomory,
 		                   ///< The GomoryHu class. You must call its
 		                   ///  run() method
 		                   ///  before initializing this iterator.
 		                   const Node& s, ///< The base node.
 		                   const Node& t,
 		                   ///< The node you want to separate from node \c s.
 		                   bool side=true
 		                   ///< If it is \c true (default) then the iterator lists
 		                   ///  the nodes of the component containing \c s,
 		                   ///  otherwise it lists the other component.
 		                   /// \note As the minimum cut is not always unique,
 		                   /// \code
 		                   /// MinCutNodeIt(gomory, s, t, true);
 		                   /// \endcode
 		                   /// and
 		                   /// \code
 		                   /// MinCutNodeIt(gomory, t, s, false);
 		                   /// \endcode
 		                   /// does not necessarily give the same set of nodes.
 		                   /// However it is ensured that
 		                   /// \code
 		                   /// MinCutNodeIt(gomory, s, t, true);
 		                   /// \endcode
 		                   /// and
 		                   /// \code
 		                   /// MinCutNodeIt(gomory, s, t, false);
 		                   /// \endcode
 		                   /// together list each node exactly once.
+		                   )
 		        : _side(side), _cut(gomory._graph)
+		      {
 		        gomory.minCutMap(s,t,_cut);
 		        for(_node_it=typename Graph::NodeIt(gomory._graph);
 		            _node_it!=INVALID && _cut[_node_it]!=_side;
 		            ++_node_it) {}
+		      }
 		      /// Conversion to \c Node
 		      /// Conversion to \c Node.
 		      ///
 		      operator typename Graph::Node() const
+		      {
 		        return _node_it;
+		      }
 		      bool operator==(Invalid) { return _node_it==INVALID; }
 		      bool operator!=(Invalid) { return _node_it!=INVALID; }
 		      /// Next node
 		      /// Next node.
 		      ///
 		      MinCutNodeIt &operator++()
+		      {
 		        for(++_node_it;_node_it!=INVALID&&_cut[_node_it]!=_side;++_node_it) {}
 		        return *this;
+		      }
 		      /// Postfix incrementation
 		      /// Postfix incrementation.
 		      ///
 		      /// \warning This incrementation
 		      /// returns a \c Node, not a \c MinCutNodeIt, as one may
 		      /// expect.
 		      typename Graph::Node operator++(int)
+		      {
 		        typename Graph::Node n=*this;
 		        ++(*this);
 		        return n;
+		      }
 		    };
 		    friend class MinCutEdgeIt;
 		    /// Iterate on the edges of a minimum cut
 		    /// This iterator class lists the edges of a minimum cut found by
 		    /// GomoryHu. Before using it, you must allocate a GomoryHu class,
 		    /// and call its \ref GomoryHu::run() "run()" method.
 		    ///
 		    /// This example computes the value of the minimum cut separating \c s from
 		    /// \c t.
 		    /// \code
 		    /// GomoruHu<Graph> gom(g, capacities);
 		    /// gom.run();
 		    /// int value=0;
 		    /// for(GomoruHu<Graph>::MinCutEdgeIt e(gom,s,t); e!=INVALID; ++e)
 		    ///   value+=capacities[e];
 		    /// \endcode
 		    /// the result will be the same as it is returned by
 		    /// \ref GomoryHu::minCutValue() "gom.minCutValue(s,t)"
 		    class MinCutEdgeIt
+		    {
 		      bool _side;
 		      const Graph &_graph;
 		      typename Graph::NodeIt _node_it;
 		      typename Graph::OutArcIt _arc_it;
 		      typename Graph::template NodeMap<bool> _cut;
 		      void step()
+		      {
 		        ++_arc_it;
 		        while(_node_it!=INVALID && _arc_it==INVALID)
+		          {
 		            for(++_node_it;_node_it!=INVALID&&!_cut[_node_it];++_node_it) {}
 		            if(_node_it!=INVALID)
 		              _arc_it=typename Graph::OutArcIt(_graph,_node_it);
+		          }
+		      }
 		    public:
 		      MinCutEdgeIt(GomoryHu const &gomory,
 		                   ///< The GomoryHu class. You must call its
 		                   ///  run() method
 		                   ///  before initializing this iterator.
 		                   const Node& s,  ///< The base node.
 		                   const Node& t,
 		                   ///< The node you want to separate from node \c s.
 		                   bool side=true
 		                   ///< If it is \c true (default) then the listed arcs
 		                   ///  will be oriented from the
 		                   ///  the nodes of the component containing \c s,
 		                   ///  otherwise they will be oriented in the opposite
 		                   ///  direction.
+		                   )
 		        : _graph(gomory._graph), _cut(_graph)
+		      {
 		        gomory.minCutMap(s,t,_cut);
 		        if(!side)
 		          for(typename Graph::NodeIt n(_graph);n!=INVALID;++n)
 		            _cut[n]=!_cut[n];
 		        for(_node_it=typename Graph::NodeIt(_graph);
 		            _node_it!=INVALID && !_cut[_node_it];
 		            ++_node_it) {}
 		        _arc_it = _node_it!=INVALID ?
 		          typename Graph::OutArcIt(_graph,_node_it) : INVALID;
 		        while(_node_it!=INVALID && _arc_it == INVALID)
+		          {
 		            for(++_node_it; _node_it!=INVALID&&!_cut[_node_it]; ++_node_it) {}
 		            if(_node_it!=INVALID)
 		              _arc_it= typename Graph::OutArcIt(_graph,_node_it);
+		          }
 		        while(_arc_it!=INVALID && _cut[_graph.target(_arc_it)]) step();
+		      }
 		      /// Conversion to \c Arc

lemon/hao_orlin.h

Show inline comments

 		/* -*- mode: C++; indent-tabs-mode: nil; -*-
+		 *
 		 * This file is a part of LEMON, a generic C++ optimization library.
+		 *
 		 * Copyright (C) 2003-2009
 		 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
 		 * (Egervary Research Group on Combinatorial Optimization, EGRES).
+		 *
 		 * Permission to use, modify and distribute this software is granted
 		 * provided that this copyright notice appears in all copies. For
 		 * precise terms see the accompanying LICENSE file.
+		 *
 		 * This software is provided "AS IS" with no warranty of any kind,
 		 * express or implied, and with no claim as to its suitability for any
 		 * purpose.
+		 *
 		 */
 		#ifndef LEMON_HAO_ORLIN_H
 		#define LEMON_HAO_ORLIN_H
 		#include <vector>
 		#include <list>
 		#include <limits>
 		#include <lemon/maps.h>
 		#include <lemon/core.h>
 		#include <lemon/tolerance.h>
 		/// \file
 		/// \ingroup min_cut
 		/// \brief Implementation of the Hao-Orlin algorithm.
 		///
 		/// Implementation of the Hao-Orlin algorithm class for testing network
 		/// reliability.
 		namespace lemon {
 		  /// \ingroup min_cut
 		  ///
 		  /// \brief %Hao-Orlin algorithm to find a minimum cut in directed graphs.
 		  ///
 		  /// Hao-Orlin calculates a minimum cut in a directed graph
 		  /// \f$D=(V,A)\f$. It takes a fixed node \f$ source \in V \f$ and
 		  /// consists of two phases: in the first phase it determines a
 		  /// minimum cut with \f$ source \f$ on the source-side (i.e. a set
 		  /// \f$ X\subsetneq V \f$ with \f$ source \in X \f$ and minimal
 		  /// out-degree) and in the second phase it determines a minimum cut
 		  /// with \f$ source \f$ on the sink-side (i.e. a set
 		  /// \f$ X\subsetneq V \f$ with \f$ source \notin X \f$ and minimal
 		  /// out-degree). Obviously, the smaller of these two cuts will be a
 		  /// minimum cut of \f$ D \f$. The algorithm is a modified
 		  /// push-relabel preflow algorithm and our implementation calculates
 		  /// the minimum cut in \f$ O(n^2\sqrt{m}) \f$ time (we use the
 		  /// highest-label rule), or in \f$O(nm)\f$ for unit capacities. The
 		  /// purpose of such algorithm is testing network reliability. For an
 		  /// undirected graph you can run just the first phase of the
 		  /// algorithm or you can use the algorithm of Nagamochi and Ibaraki
 		  /// which solves the undirected problem in
 		  /// \f$ O(nm + n^2 \log n) \f$ time: it is implemented in the
 		  /// NagamochiIbaraki algorithm class.
 		  ///
 		  /// \param GR The digraph class the algorithm runs on.
 		  /// \param CAP An arc map of capacities which can be any numreric type.
 		  /// The default type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
 		  /// \param TOL Tolerance class for handling inexact computations. The
 		  /// default tolerance type is \ref Tolerance "Tolerance<CAP::Value>".
 		#ifdef DOXYGEN
 		  template <typename GR, typename CAP, typename TOL>
 		#else
 		  template <typename GR,
 		            typename CAP = typename GR::template ArcMap<int>,
 		            typename TOL = Tolerance<typename CAP::Value> >
 		#endif
 		  class HaoOrlin {
 		  private:
 		    typedef GR Digraph;
 		    typedef CAP CapacityMap;
 		    typedef TOL Tolerance;
 		    typedef typename CapacityMap::Value Value;
 		    TEMPLATE_GRAPH_TYPEDEFS(Digraph);
 		    const Digraph& _graph;
 		    const CapacityMap* _capacity;
 		    typedef typename Digraph::template ArcMap<Value> FlowMap;
 		    FlowMap* _flow;
 		    Node _source;
 		    int _node_num;
 		    // Bucketing structure
 		    std::vector<Node> _first, _last;
 		    typename Digraph::template NodeMap<Node>* _next;
 		    typename Digraph::template NodeMap<Node>* _prev;
 		    typename Digraph::template NodeMap<bool>* _active;
 		    typename Digraph::template NodeMap<int>* _bucket;
 		    std::vector<bool> _dormant;
 		    std::list<std::list<int> > _sets;
 		    std::list<int>::iterator _highest;
 		    typedef typename Digraph::template NodeMap<Value> ExcessMap;
 		    ExcessMap* _excess;
 		    typedef typename Digraph::template NodeMap<bool> SourceSetMap;
 		    SourceSetMap* _source_set;
 		    Value _min_cut;
 		    typedef typename Digraph::template NodeMap<bool> MinCutMap;
 		    MinCutMap* _min_cut_map;
 		    Tolerance _tolerance;
 		  public:
 		    /// \brief Constructor
 		    ///
 		    /// Constructor of the algorithm class.
 		    HaoOrlin(const Digraph& graph, const CapacityMap& capacity,
 		             const Tolerance& tolerance = Tolerance()) :
 		      _graph(graph), _capacity(&capacity), _flow(0), _source(),
 		      _node_num(), _first(), _last(), _next(0), _prev(0),
 		      _active(0), _bucket(0), _dormant(), _sets(), _highest(),
 		      _excess(0), _source_set(0), _min_cut(), _min_cut_map(0),
 		      _tolerance(tolerance) {}
 		    ~HaoOrlin() {
 		      if (_min_cut_map) {
 		        delete _min_cut_map;
+		      }
 		      if (_source_set) {
 		        delete _source_set;
+		      }
 		      if (_excess) {
 		        delete _excess;
+		      }
 		      if (_next) {
 		        delete _next;
+		      }
 		      if (_prev) {
 		        delete _prev;
+		      }
 		      if (_active) {
 		        delete _active;
+		      }
 		      if (_bucket) {
 		        delete _bucket;
+		      }
 		      if (_flow) {
 		        delete _flow;
+		      }
+		    }
 		  private:
 		    void activate(const Node& i) {
-		      _active->set(i, true);
+		      (*_active)[i] = true;
 		      int bucket = (*_bucket)[i];
 		      if ((*_prev)[i] == INVALID || (*_active)[(*_prev)[i]]) return;
 		      //unlace
-		      _next->set((*_prev)[i], (*_next)[i]);
+		      (*_next)[(*_prev)[i]] = (*_next)[i];
 		      if ((*_next)[i] != INVALID) {
-		        _prev->set((*_next)[i], (*_prev)[i]);
+		        (*_prev)[(*_next)[i]] = (*_prev)[i];
 		      } else {
 		        _last[bucket] = (*_prev)[i];
+		      }
 		      //lace
 		      _next->set(i, _first[bucket]);
 		      _prev->set(_first[bucket], i);
-		      _prev->set(i, INVALID);
+		      (*_next)[i] = _first[bucket];
 		      (*_prev)[_first[bucket]] = i;
 		      (*_prev)[i] = INVALID;
 		      _first[bucket] = i;
+		    }
 		    void deactivate(const Node& i) {
-		      _active->set(i, false);
+		      (*_active)[i] = false;
 		      int bucket = (*_bucket)[i];
 		      if ((*_next)[i] == INVALID || !(*_active)[(*_next)[i]]) return;
 		      //unlace
-		      _prev->set((*_next)[i], (*_prev)[i]);
+		      (*_prev)[(*_next)[i]] = (*_prev)[i];
 		      if ((*_prev)[i] != INVALID) {
-		        _next->set((*_prev)[i], (*_next)[i]);
+		        (*_next)[(*_prev)[i]] = (*_next)[i];
 		      } else {
 		        _first[bucket] = (*_next)[i];
+		      }
 		      //lace
 		      _prev->set(i, _last[bucket]);
 		      _next->set(_last[bucket], i);
-		      _next->set(i, INVALID);
+		      (*_prev)[i] = _last[bucket];
 		      (*_next)[_last[bucket]] = i;
 		      (*_next)[i] = INVALID;
 		      _last[bucket] = i;
+		    }
 		    void addItem(const Node& i, int bucket) {
 		      (*_bucket)[i] = bucket;
 		      if (_last[bucket] != INVALID) {
 		        _prev->set(i, _last[bucket]);
 		        _next->set(_last[bucket], i);
-		        _next->set(i, INVALID);
+		        (*_prev)[i] = _last[bucket];
 		        (*_next)[_last[bucket]] = i;
 		        (*_next)[i] = INVALID;
 		        _last[bucket] = i;
 		      } else {
-		        _prev->set(i, INVALID);
+		        (*_prev)[i] = INVALID;
 		        _first[bucket] = i;
-		        _next->set(i, INVALID);
+		        (*_next)[i] = INVALID;
 		        _last[bucket] = i;
+		      }
+		    }
 		    void findMinCutOut() {
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
-		        _excess->set(n, 0);
+		        (*_excess)[n] = 0;
+		      }
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
-		        _flow->set(a, 0);
+		        (*_flow)[a] = 0;
+		      }
 		      int bucket_num = 0;
 		      std::vector<Node> queue(_node_num);
 		      int qfirst = 0, qlast = 0, qsep = 0;
+		      {
 		        typename Digraph::template NodeMap<bool> reached(_graph, false);
-		        reached.set(_source, true);
+		        reached[_source] = true;
 		        bool first_set = true;
 		        for (NodeIt t(_graph); t != INVALID; ++t) {
 		          if (reached[t]) continue;
 		          _sets.push_front(std::list<int>());
 		          queue[qlast++] = t;
-		          reached.set(t, true);
+		          reached[t] = true;
 		          while (qfirst != qlast) {
 		            if (qsep == qfirst) {
 		              ++bucket_num;
 		              _sets.front().push_front(bucket_num);
 		              _dormant[bucket_num] = !first_set;
 		              _first[bucket_num] = _last[bucket_num] = INVALID;
 		              qsep = qlast;
+		            }
 		            Node n = queue[qfirst++];
 		            addItem(n, bucket_num);
 		            for (InArcIt a(_graph, n); a != INVALID; ++a) {
 		              Node u = _graph.source(a);
 		              if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
-		                reached.set(u, true);
+		                reached[u] = true;
 		                queue[qlast++] = u;
+		              }
+		            }
+		          }
 		          first_set = false;
+		        }
 		        ++bucket_num;
-		        _bucket->set(_source, 0);
+		        (*_bucket)[_source] = 0;
 		        _dormant[0] = true;
+		      }
-		      _source_set->set(_source, true);
+		      (*_source_set)[_source] = true;
 		      Node target = _last[_sets.back().back()];
+		      {
 		        for (OutArcIt a(_graph, _source); a != INVALID; ++a) {
 		          if (_tolerance.positive((*_capacity)[a])) {
 		            Node u = _graph.target(a);
 		            _flow->set(a, (*_capacity)[a]);
 		            _excess->set(u, (*_excess)[u] + (*_capacity)[a]);
 		            (*_flow)[a] = (*_capacity)[a];
 		            (*_excess)[u] += (*_capacity)[a];
 		            if (!(*_active)[u] && u != _source) {
 		              activate(u);
+		            }
+		          }
+		        }
 		        if ((*_active)[target]) {
 		          deactivate(target);
+		        }
 		        _highest = _sets.back().begin();
 		        while (_highest != _sets.back().end() &&
 		               !(*_active)[_first[*_highest]]) {
 		          ++_highest;
+		        }
+		      }
 		      while (true) {
 		        while (_highest != _sets.back().end()) {
 		          Node n = _first[*_highest];
 		          Value excess = (*_excess)[n];
 		          int next_bucket = _node_num;
 		          int under_bucket;
 		          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
 		            under_bucket = -1;
 		          } else {
 		            under_bucket = *(++std::list<int>::iterator(_highest));
+		          }
 		          for (OutArcIt a(_graph, n); a != INVALID; ++a) {
 		            Node v = _graph.target(a);
 		            if (_dormant[(*_bucket)[v]]) continue;
 		            Value rem = (*_capacity)[a] - (*_flow)[a];
 		            if (!_tolerance.positive(rem)) continue;
 		            if ((*_bucket)[v] == under_bucket) {
 		              if (!(*_active)[v] && v != target) {
 		                activate(v);
+		              }
 		              if (!_tolerance.less(rem, excess)) {
 		                _flow->set(a, (*_flow)[a] + excess);
 		                _excess->set(v, (*_excess)[v] + excess);
 		                (*_flow)[a] += excess;
 		                (*_excess)[v] += excess;
 		                excess = 0;
 		                goto no_more_push;
 		              } else {
 		                excess -= rem;
 		                _excess->set(v, (*_excess)[v] + rem);
 		                _flow->set(a, (*_capacity)[a]);
 		                (*_excess)[v] += rem;
 		                (*_flow)[a] = (*_capacity)[a];
+		              }
 		            } else if (next_bucket > (*_bucket)[v]) {
 		              next_bucket = (*_bucket)[v];
+		            }
+		          }
 		          for (InArcIt a(_graph, n); a != INVALID; ++a) {
 		            Node v = _graph.source(a);
 		            if (_dormant[(*_bucket)[v]]) continue;
 		            Value rem = (*_flow)[a];
 		            if (!_tolerance.positive(rem)) continue;
 		            if ((*_bucket)[v] == under_bucket) {
 		              if (!(*_active)[v] && v != target) {
 		                activate(v);
+		              }
 		              if (!_tolerance.less(rem, excess)) {
 		                _flow->set(a, (*_flow)[a] - excess);
 		                _excess->set(v, (*_excess)[v] + excess);
 		                (*_flow)[a] -= excess;
 		                (*_excess)[v] += excess;
 		                excess = 0;
 		                goto no_more_push;
 		              } else {
 		                excess -= rem;
 		                _excess->set(v, (*_excess)[v] + rem);
 		                _flow->set(a, 0);
 		                (*_excess)[v] += rem;
 		                (*_flow)[a] = 0;
+		              }
 		            } else if (next_bucket > (*_bucket)[v]) {
 		              next_bucket = (*_bucket)[v];
+		            }
+		          }
 		        no_more_push:
-		          _excess->set(n, excess);
+		          (*_excess)[n] = excess;
 		          if (excess != 0) {
 		            if ((*_next)[n] == INVALID) {
 		              typename std::list<std::list<int> >::iterator new_set =
 		                _sets.insert(--_sets.end(), std::list<int>());
 		              new_set->splice(new_set->end(), _sets.back(),
 		                              _sets.back().begin(), ++_highest);
 		              for (std::list<int>::iterator it = new_set->begin();
 		                   it != new_set->end(); ++it) {
 		                _dormant[*it] = true;
+		              }
 		              while (_highest != _sets.back().end() &&
 		                     !(*_active)[_first[*_highest]]) {
 		                ++_highest;
+		              }
 		            } else if (next_bucket == _node_num) {
 		              _first[(*_bucket)[n]] = (*_next)[n];
-		              _prev->set((*_next)[n], INVALID);
+		              (*_prev)[(*_next)[n]] = INVALID;
 		              std::list<std::list<int> >::iterator new_set =
 		                _sets.insert(--_sets.end(), std::list<int>());
 		              new_set->push_front(bucket_num);
-		              _bucket->set(n, bucket_num);
+		              (*_bucket)[n] = bucket_num;
 		              _first[bucket_num] = _last[bucket_num] = n;
 		              _next->set(n, INVALID);
 		              _prev->set(n, INVALID);
 		              (*_next)[n] = INVALID;
 		              (*_prev)[n] = INVALID;
 		              _dormant[bucket_num] = true;
 		              ++bucket_num;
 		              while (_highest != _sets.back().end() &&
 		                     !(*_active)[_first[*_highest]]) {
 		                ++_highest;
+		              }
 		            } else {
 		              _first[*_highest] = (*_next)[n];
-		              _prev->set((*_next)[n], INVALID);
+		              (*_prev)[(*_next)[n]] = INVALID;
 		              while (next_bucket != *_highest) {
 		                --_highest;
+		              }
 		              if (_highest == _sets.back().begin()) {
 		                _sets.back().push_front(bucket_num);
 		                _dormant[bucket_num] = false;
 		                _first[bucket_num] = _last[bucket_num] = INVALID;
 		                ++bucket_num;
+		              }
 		              --_highest;
 		              _bucket->set(n, *_highest);
 		              _next->set(n, _first[*_highest]);
 		              (*_bucket)[n] = *_highest;
 		              (*_next)[n] = _first[*_highest];
 		              if (_first[*_highest] != INVALID) {
-		                _prev->set(_first[*_highest], n);
+		                (*_prev)[_first[*_highest]] = n;
 		              } else {
 		                _last[*_highest] = n;
+		              }
 		              _first[*_highest] = n;
+		            }
 		          } else {
 		            deactivate(n);
 		            if (!(*_active)[_first[*_highest]]) {
 		              ++_highest;
 		              if (_highest != _sets.back().end() &&
 		                  !(*_active)[_first[*_highest]]) {
 		                _highest = _sets.back().end();
+		              }
+		            }
+		          }
+		        }
 		        if ((*_excess)[target] < _min_cut) {
 		          _min_cut = (*_excess)[target];
 		          for (NodeIt i(_graph); i != INVALID; ++i) {
-		            _min_cut_map->set(i, true);
+		            (*_min_cut_map)[i] = true;
+		          }
 		          for (std::list<int>::iterator it = _sets.back().begin();
 		               it != _sets.back().end(); ++it) {
 		            Node n = _first[*it];
 		            while (n != INVALID) {
-		              _min_cut_map->set(n, false);
+		              (*_min_cut_map)[n] = false;
 		              n = (*_next)[n];
+		            }
+		          }
+		        }
+		        {
 		          Node new_target;
 		          if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
 		            if ((*_next)[target] == INVALID) {
 		              _last[(*_bucket)[target]] = (*_prev)[target];
 		              new_target = (*_prev)[target];
 		            } else {
-		              _prev->set((*_next)[target], (*_prev)[target]);
+		              (*_prev)[(*_next)[target]] = (*_prev)[target];
 		              new_target = (*_next)[target];
+		            }
 		            if ((*_prev)[target] == INVALID) {
 		              _first[(*_bucket)[target]] = (*_next)[target];
 		            } else {
-		              _next->set((*_prev)[target], (*_next)[target]);
+		              (*_next)[(*_prev)[target]] = (*_next)[target];
+		            }
 		          } else {
 		            _sets.back().pop_back();
 		            if (_sets.back().empty()) {
 		              _sets.pop_back();
 		              if (_sets.empty())
 		                break;
 		              for (std::list<int>::iterator it = _sets.back().begin();
 		                   it != _sets.back().end(); ++it) {
 		                _dormant[*it] = false;
+		              }
+		            }
 		            new_target = _last[_sets.back().back()];
+		          }
-		          _bucket->set(target, 0);
+		          (*_bucket)[target] = 0;
-		          _source_set->set(target, true);
+		          (*_source_set)[target] = true;
 		          for (OutArcIt a(_graph, target); a != INVALID; ++a) {
 		            Value rem = (*_capacity)[a] - (*_flow)[a];
 		            if (!_tolerance.positive(rem)) continue;
 		            Node v = _graph.target(a);
 		            if (!(*_active)[v] && !(*_source_set)[v]) {
 		              activate(v);
+		            }
 		            _excess->set(v, (*_excess)[v] + rem);
 		            _flow->set(a, (*_capacity)[a]);
 		            (*_excess)[v] += rem;
 		            (*_flow)[a] = (*_capacity)[a];
+		          }
 		          for (InArcIt a(_graph, target); a != INVALID; ++a) {
 		            Value rem = (*_flow)[a];
 		            if (!_tolerance.positive(rem)) continue;
 		            Node v = _graph.source(a);
 		            if (!(*_active)[v] && !(*_source_set)[v]) {
 		              activate(v);
+		            }
 		            _excess->set(v, (*_excess)[v] + rem);
 		            _flow->set(a, 0);
 		            (*_excess)[v] += rem;
 		            (*_flow)[a] = 0;
+		          }
 		          target = new_target;
 		          if ((*_active)[target]) {
 		            deactivate(target);
+		          }
 		          _highest = _sets.back().begin();
 		          while (_highest != _sets.back().end() &&
 		                 !(*_active)[_first[*_highest]]) {
 		            ++_highest;
+		          }
+		        }
+		      }
+		    }
 		    void findMinCutIn() {
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
-		        _excess->set(n, 0);
+		        (*_excess)[n] = 0;
+		      }
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
-		        _flow->set(a, 0);
+		        (*_flow)[a] = 0;
+		      }
 		      int bucket_num = 0;
 		      std::vector<Node> queue(_node_num);
 		      int qfirst = 0, qlast = 0, qsep = 0;
+		      {
 		        typename Digraph::template NodeMap<bool> reached(_graph, false);
-		        reached.set(_source, true);
+		        reached[_source] = true;
 		        bool first_set = true;
 		        for (NodeIt t(_graph); t != INVALID; ++t) {
 		          if (reached[t]) continue;
 		          _sets.push_front(std::list<int>());
 		          queue[qlast++] = t;
-		          reached.set(t, true);
+		          reached[t] = true;
 		          while (qfirst != qlast) {
 		            if (qsep == qfirst) {
 		              ++bucket_num;
 		              _sets.front().push_front(bucket_num);
 		              _dormant[bucket_num] = !first_set;
 		              _first[bucket_num] = _last[bucket_num] = INVALID;
 		              qsep = qlast;
+		            }
 		            Node n = queue[qfirst++];
 		            addItem(n, bucket_num);
 		            for (OutArcIt a(_graph, n); a != INVALID; ++a) {
 		              Node u = _graph.target(a);
 		              if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
-		                reached.set(u, true);
+		                reached[u] = true;
 		                queue[qlast++] = u;
+		              }
+		            }
+		          }
 		          first_set = false;
+		        }
 		        ++bucket_num;
-		        _bucket->set(_source, 0);
+		        (*_bucket)[_source] = 0;
 		        _dormant[0] = true;
+		      }
-		      _source_set->set(_source, true);
+		      (*_source_set)[_source] = true;
 		      Node target = _last[_sets.back().back()];
+		      {
 		        for (InArcIt a(_graph, _source); a != INVALID; ++a) {
 		          if (_tolerance.positive((*_capacity)[a])) {
 		            Node u = _graph.source(a);
 		            _flow->set(a, (*_capacity)[a]);
 		            _excess->set(u, (*_excess)[u] + (*_capacity)[a]);
 		            (*_flow)[a] = (*_capacity)[a];
 		            (*_excess)[u] += (*_capacity)[a];
 		            if (!(*_active)[u] && u != _source) {
 		              activate(u);
+		            }
+		          }
+		        }
 		        if ((*_active)[target]) {
 		          deactivate(target);
+		        }
 		        _highest = _sets.back().begin();
 		        while (_highest != _sets.back().end() &&
 		               !(*_active)[_first[*_highest]]) {
 		          ++_highest;
+		        }
+		      }
 		      while (true) {
 		        while (_highest != _sets.back().end()) {
 		          Node n = _first[*_highest];
 		          Value excess = (*_excess)[n];
 		          int next_bucket = _node_num;
 		          int under_bucket;
 		          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
 		            under_bucket = -1;
 		          } else {
 		            under_bucket = *(++std::list<int>::iterator(_highest));
+		          }
 		          for (InArcIt a(_graph, n); a != INVALID; ++a) {
 		            Node v = _graph.source(a);
 		            if (_dormant[(*_bucket)[v]]) continue;
 		            Value rem = (*_capacity)[a] - (*_flow)[a];
 		            if (!_tolerance.positive(rem)) continue;
 		            if ((*_bucket)[v] == under_bucket) {
 		              if (!(*_active)[v] && v != target) {
 		                activate(v);
+		              }
 		              if (!_tolerance.less(rem, excess)) {
 		                _flow->set(a, (*_flow)[a] + excess);
 		                _excess->set(v, (*_excess)[v] + excess);
 		                (*_flow)[a] += excess;
 		                (*_excess)[v] += excess;
 		                excess = 0;
 		                goto no_more_push;
 		              } else {
 		                excess -= rem;
 		                _excess->set(v, (*_excess)[v] + rem);
 		                _flow->set(a, (*_capacity)[a]);
 		                (*_excess)[v] += rem;
 		                (*_flow)[a] = (*_capacity)[a];
+		              }
 		            } else if (next_bucket > (*_bucket)[v]) {
 		              next_bucket = (*_bucket)[v];
+		            }
+		          }
 		          for (OutArcIt a(_graph, n); a != INVALID; ++a) {
 		            Node v = _graph.target(a);
 		            if (_dormant[(*_bucket)[v]]) continue;
 		            Value rem = (*_flow)[a];
 		            if (!_tolerance.positive(rem)) continue;
 		            if ((*_bucket)[v] == under_bucket) {
 		              if (!(*_active)[v] && v != target) {
 		                activate(v);
+		              }
 		              if (!_tolerance.less(rem, excess)) {
 		                _flow->set(a, (*_flow)[a] - excess);
 		                _excess->set(v, (*_excess)[v] + excess);
 		                (*_flow)[a] -= excess;
 		                (*_excess)[v] += excess;
 		                excess = 0;
 		                goto no_more_push;
 		              } else {
 		                excess -= rem;
 		                _excess->set(v, (*_excess)[v] + rem);
 		                _flow->set(a, 0);
 		                (*_excess)[v] += rem;
 		                (*_flow)[a] = 0;
+		              }
 		            } else if (next_bucket > (*_bucket)[v]) {
 		              next_bucket = (*_bucket)[v];
+		            }
+		          }
 		        no_more_push:
-		          _excess->set(n, excess);
+		          (*_excess)[n] = excess;
 		          if (excess != 0) {
 		            if ((*_next)[n] == INVALID) {
 		              typename std::list<std::list<int> >::iterator new_set =
 		                _sets.insert(--_sets.end(), std::list<int>());
 		              new_set->splice(new_set->end(), _sets.back(),
 		                              _sets.back().begin(), ++_highest);
 		              for (std::list<int>::iterator it = new_set->begin();
 		                   it != new_set->end(); ++it) {
 		                _dormant[*it] = true;
+		              }
 		              while (_highest != _sets.back().end() &&
 		                     !(*_active)[_first[*_highest]]) {
 		                ++_highest;
+		              }
 		            } else if (next_bucket == _node_num) {
 		              _first[(*_bucket)[n]] = (*_next)[n];
-		              _prev->set((*_next)[n], INVALID);
+		              (*_prev)[(*_next)[n]] = INVALID;
 		              std::list<std::list<int> >::iterator new_set =
 		                _sets.insert(--_sets.end(), std::list<int>());
 		              new_set->push_front(bucket_num);
-		              _bucket->set(n, bucket_num);
+		              (*_bucket)[n] = bucket_num;
 		              _first[bucket_num] = _last[bucket_num] = n;
 		              _next->set(n, INVALID);
 		              _prev->set(n, INVALID);
 		              (*_next)[n] = INVALID;
 		              (*_prev)[n] = INVALID;
 		              _dormant[bucket_num] = true;
 		              ++bucket_num;
 		              while (_highest != _sets.back().end() &&
 		                     !(*_active)[_first[*_highest]]) {
 		                ++_highest;
+		              }
 		            } else {
 		              _first[*_highest] = (*_next)[n];
-		              _prev->set((*_next)[n], INVALID);
+		              (*_prev)[(*_next)[n]] = INVALID;
 		              while (next_bucket != *_highest) {
 		                --_highest;
+		              }
 		              if (_highest == _sets.back().begin()) {
 		                _sets.back().push_front(bucket_num);
 		                _dormant[bucket_num] = false;
 		                _first[bucket_num] = _last[bucket_num] = INVALID;
 		                ++bucket_num;
+		              }
 		              --_highest;
 		              _bucket->set(n, *_highest);
 		              _next->set(n, _first[*_highest]);
 		              (*_bucket)[n] = *_highest;
 		              (*_next)[n] = _first[*_highest];
 		              if (_first[*_highest] != INVALID) {
-		                _prev->set(_first[*_highest], n);
+		                (*_prev)[_first[*_highest]] = n;
 		              } else {
 		                _last[*_highest] = n;
+		              }
 		              _first[*_highest] = n;
+		            }
 		          } else {
 		            deactivate(n);
 		            if (!(*_active)[_first[*_highest]]) {
 		              ++_highest;
 		              if (_highest != _sets.back().end() &&
 		                  !(*_active)[_first[*_highest]]) {
 		                _highest = _sets.back().end();
+		              }
+		            }
+		          }
+		        }
 		        if ((*_excess)[target] < _min_cut) {
 		          _min_cut = (*_excess)[target];
 		          for (NodeIt i(_graph); i != INVALID; ++i) {
-		            _min_cut_map->set(i, false);
+		            (*_min_cut_map)[i] = false;
+		          }
 		          for (std::list<int>::iterator it = _sets.back().begin();
 		               it != _sets.back().end(); ++it) {
 		            Node n = _first[*it];
 		            while (n != INVALID) {
-		              _min_cut_map->set(n, true);
+		              (*_min_cut_map)[n] = true;
 		              n = (*_next)[n];
+		            }
+		          }
+		        }
+		        {
 		          Node new_target;
 		          if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
 		            if ((*_next)[target] == INVALID) {
 		              _last[(*_bucket)[target]] = (*_prev)[target];
 		              new_target = (*_prev)[target];
 		            } else {
-		              _prev->set((*_next)[target], (*_prev)[target]);
+		              (*_prev)[(*_next)[target]] = (*_prev)[target];
 		              new_target = (*_next)[target];
+		            }
 		            if ((*_prev)[target] == INVALID) {
 		              _first[(*_bucket)[target]] = (*_next)[target];
 		            } else {
-		              _next->set((*_prev)[target], (*_next)[target]);
+		              (*_next)[(*_prev)[target]] = (*_next)[target];
+		            }
 		          } else {
 		            _sets.back().pop_back();
 		            if (_sets.back().empty()) {
 		              _sets.pop_back();
 		              if (_sets.empty())
 		                break;
 		              for (std::list<int>::iterator it = _sets.back().begin();
 		                   it != _sets.back().end(); ++it) {
 		                _dormant[*it] = false;
+		              }
+		            }
 		            new_target = _last[_sets.back().back()];
+		          }
-		          _bucket->set(target, 0);
+		          (*_bucket)[target] = 0;
-		          _source_set->set(target, true);
+		          (*_source_set)[target] = true;
 		          for (InArcIt a(_graph, target); a != INVALID; ++a) {
 		            Value rem = (*_capacity)[a] - (*_flow)[a];
 		            if (!_tolerance.positive(rem)) continue;
 		            Node v = _graph.source(a);
 		            if (!(*_active)[v] && !(*_source_set)[v]) {
 		              activate(v);
+		            }
 		            _excess->set(v, (*_excess)[v] + rem);
 		            _flow->set(a, (*_capacity)[a]);
 		            (*_excess)[v] += rem;
 		            (*_flow)[a] = (*_capacity)[a];
+		          }
 		          for (OutArcIt a(_graph, target); a != INVALID; ++a) {
 		            Value rem = (*_flow)[a];
 		            if (!_tolerance.positive(rem)) continue;
 		            Node v = _graph.target(a);
 		            if (!(*_active)[v] && !(*_source_set)[v]) {
 		              activate(v);
+		            }
 		            _excess->set(v, (*_excess)[v] + rem);
 		            _flow->set(a, 0);
 		            (*_excess)[v] += rem;
 		            (*_flow)[a] = 0;
+		          }
 		          target = new_target;
 		          if ((*_active)[target]) {
 		            deactivate(target);
+		          }
 		          _highest = _sets.back().begin();
 		          while (_highest != _sets.back().end() &&
 		                 !(*_active)[_first[*_highest]]) {
 		            ++_highest;
+		          }
+		        }
+		      }
+		    }
 		  public:
 		    /// \name Execution control
 		    /// The simplest way to execute the algorithm is to use
 		    /// one of the member functions called \ref run().
 		    /// \n
 		    /// If you need more control on the execution,
 		    /// first you must call \ref init(), then the \ref calculateIn() or
 		    /// \ref calculateOut() functions.
 		    /// @{
 		    /// \brief Initializes the internal data structures.
 		    ///
 		    /// Initializes the internal data structures. It creates
 		    /// the maps, residual graph adaptors and some bucket structures
 		    /// for the algorithm.
 		    void init() {
 		      init(NodeIt(_graph));
+		    }
 		    /// \brief Initializes the internal data structures.
 		    ///
 		    /// Initializes the internal data structures. It creates
 		    /// the maps, residual graph adaptor and some bucket structures
 		    /// for the algorithm. Node \c source  is used as the push-relabel
 		    /// algorithm's source.
 		    void init(const Node& source) {
 		      _source = source;
 		      _node_num = countNodes(_graph);
 		      _first.resize(_node_num);
 		      _last.resize(_node_num);
 		      _dormant.resize(_node_num);
 		      if (!_flow) {
 		        _flow = new FlowMap(_graph);
+		      }
 		      if (!_next) {
 		        _next = new typename Digraph::template NodeMap<Node>(_graph);
+		      }
 		      if (!_prev) {
 		        _prev = new typename Digraph::template NodeMap<Node>(_graph);
+		      }
 		      if (!_active) {
 		        _active = new typename Digraph::template NodeMap<bool>(_graph);
+		      }
 		      if (!_bucket) {
 		        _bucket = new typename Digraph::template NodeMap<int>(_graph);
+		      }
 		      if (!_excess) {
 		        _excess = new ExcessMap(_graph);
+		      }
 		      if (!_source_set) {
 		        _source_set = new SourceSetMap(_graph);
+		      }
 		      if (!_min_cut_map) {
 		        _min_cut_map = new MinCutMap(_graph);
+		      }
 		      _min_cut = std::numeric_limits<Value>::max();
+		    }
 		    /// \brief Calculates a minimum cut with \f$ source \f$ on the
 		    /// source-side.
 		    ///
 		    /// Calculates a minimum cut with \f$ source \f$ on the
 		    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with
 		    /// \f$ source \in X \f$ and minimal out-degree).
 		    void calculateOut() {
 		      findMinCutOut();
+		    }
 		    /// \brief Calculates a minimum cut with \f$ source \f$ on the
 		    /// target-side.
 		    ///
 		    /// Calculates a minimum cut with \f$ source \f$ on the
 		    /// target-side (i.e. a set \f$ X\subsetneq V \f$ with
 		    /// \f$ source \in X \f$ and minimal out-degree).
 		    void calculateIn() {
 		      findMinCutIn();
+		    }
 		    /// \brief Runs the algorithm.
 		    ///
 		    /// Runs the algorithm. It finds nodes \c source and \c target
 		    /// arbitrarily and then calls \ref init(), \ref calculateOut()
 		    /// and \ref calculateIn().
 		    void run() {
 		      init();
 		      calculateOut();
 		      calculateIn();
+		    }
 		    /// \brief Runs the algorithm.
 		    ///
 		    /// Runs the algorithm. It uses the given \c source node, finds a
 		    /// proper \c target and then calls the \ref init(), \ref
 		    /// calculateOut() and \ref calculateIn().
 		    void run(const Node& s) {
 		      init(s);
 		      calculateOut();
 		      calculateIn();
+		    }
 		    /// @}
 		    /// \name Query Functions
 		    /// The result of the %HaoOrlin algorithm
 		    /// can be obtained using these functions.
 		    /// \n
 		    /// Before using these functions, either \ref run(), \ref
 		    /// calculateOut() or \ref calculateIn() must be called.
 		    /// @{
 		    /// \brief Returns the value of the minimum value cut.
 		    ///
 		    /// Returns the value of the minimum value cut.
 		    Value minCutValue() const {
 		      return _min_cut;
+		    }
 		    /// \brief Returns a minimum cut.
 		    ///
 		    /// Sets \c nodeMap to the characteristic vector of a minimum
 		    /// value cut: it will give a nonempty set \f$ X\subsetneq V \f$
 		    /// with minimal out-degree (i.e. \c nodeMap will be true exactly
 		    /// for the nodes of \f$ X \f$).  \pre nodeMap should be a
 		    /// bool-valued node-map.
 		    template <typename NodeMap>
 		    Value minCutMap(NodeMap& nodeMap) const {
 		      for (NodeIt it(_graph); it != INVALID; ++it) {
 		        nodeMap.set(it, (*_min_cut_map)[it]);
+		      }
 		      return _min_cut;
+		    }
 		    /// @}
 		  }; //class HaoOrlin
 		} //namespace lemon
 		#endif //LEMON_HAO_ORLIN_H

lemon/max_matching.h

Show inline comments

@@ -93,615 +93,615 @@
 		    EarMap* _ear;
 		    IntNodeMap* _blossom_set_index;
 		    BlossomSet* _blossom_set;
 		    NodeIntMap* _blossom_rep;
 		    IntNodeMap* _tree_set_index;
 		    TreeSet* _tree_set;
 		    NodeQueue _node_queue;
 		    int _process, _postpone, _last;
 		    int _node_num;
 		  private:
 		    void createStructures() {
 		      _node_num = countNodes(_graph);
 		      if (!_matching) {
 		        _matching = new MatchingMap(_graph);
+		      }
 		      if (!_status) {
 		        _status = new StatusMap(_graph);
+		      }
 		      if (!_ear) {
 		        _ear = new EarMap(_graph);
+		      }
 		      if (!_blossom_set) {
 		        _blossom_set_index = new IntNodeMap(_graph);
 		        _blossom_set = new BlossomSet(*_blossom_set_index);
+		      }
 		      if (!_blossom_rep) {
 		        _blossom_rep = new NodeIntMap(_node_num);
+		      }
 		      if (!_tree_set) {
 		        _tree_set_index = new IntNodeMap(_graph);
 		        _tree_set = new TreeSet(*_tree_set_index);
+		      }
 		      _node_queue.resize(_node_num);
+		    }
 		    void destroyStructures() {
 		      if (_matching) {
 		        delete _matching;
+		      }
 		      if (_status) {
 		        delete _status;
+		      }
 		      if (_ear) {
 		        delete _ear;
+		      }
 		      if (_blossom_set) {
 		        delete _blossom_set;
 		        delete _blossom_set_index;
+		      }
 		      if (_blossom_rep) {
 		        delete _blossom_rep;
+		      }
 		      if (_tree_set) {
 		        delete _tree_set_index;
 		        delete _tree_set;
+		      }
+		    }
 		    void processDense(const Node& n) {
 		      _process = _postpone = _last = 0;
 		      _node_queue[_last++] = n;
 		      while (_process != _last) {
 		        Node u = _node_queue[_process++];
 		        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
 		          Node v = _graph.target(a);
 		          if ((*_status)[v] == MATCHED) {
 		            extendOnArc(a);
 		          } else if ((*_status)[v] == UNMATCHED) {
 		            augmentOnArc(a);
 		            return;
+		          }
+		        }
+		      }
 		      while (_postpone != _last) {
 		        Node u = _node_queue[_postpone++];
 		        for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
 		          Node v = _graph.target(a);
 		          if ((*_status)[v] == EVEN) {
 		            if (_blossom_set->find(u) != _blossom_set->find(v)) {
 		              shrinkOnEdge(a);
+		            }
+		          }
 		          while (_process != _last) {
 		            Node w = _node_queue[_process++];
 		            for (OutArcIt b(_graph, w); b != INVALID; ++b) {
 		              Node x = _graph.target(b);
 		              if ((*_status)[x] == MATCHED) {
 		                extendOnArc(b);
 		              } else if ((*_status)[x] == UNMATCHED) {
 		                augmentOnArc(b);
 		                return;
+		              }
+		            }
+		          }
+		        }
+		      }
+		    }
 		    void processSparse(const Node& n) {
 		      _process = _last = 0;
 		      _node_queue[_last++] = n;
 		      while (_process != _last) {
 		        Node u = _node_queue[_process++];
 		        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
 		          Node v = _graph.target(a);
 		          if ((*_status)[v] == EVEN) {
 		            if (_blossom_set->find(u) != _blossom_set->find(v)) {
 		              shrinkOnEdge(a);
+		            }
 		          } else if ((*_status)[v] == MATCHED) {
 		            extendOnArc(a);
 		          } else if ((*_status)[v] == UNMATCHED) {
 		            augmentOnArc(a);
 		            return;
+		          }
+		        }
+		      }
+		    }
 		    void shrinkOnEdge(const Edge& e) {
 		      Node nca = INVALID;
+		      {
 		        std::set<Node> left_set, right_set;
 		        Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
 		        left_set.insert(left);
 		        Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
 		        right_set.insert(right);
 		        while (true) {
 		          if ((*_matching)[left] == INVALID) break;
 		          left = _graph.target((*_matching)[left]);
 		          left = (*_blossom_rep)[_blossom_set->
 		                                 find(_graph.target((*_ear)[left]))];
 		          if (right_set.find(left) != right_set.end()) {
 		            nca = left;
 		            break;
+		          }
 		          left_set.insert(left);
 		          if ((*_matching)[right] == INVALID) break;
 		          right = _graph.target((*_matching)[right]);
 		          right = (*_blossom_rep)[_blossom_set->
 		                                  find(_graph.target((*_ear)[right]))];
 		          if (left_set.find(right) != left_set.end()) {
 		            nca = right;
 		            break;
+		          }
 		          right_set.insert(right);
+		        }
 		        if (nca == INVALID) {
 		          if ((*_matching)[left] == INVALID) {
 		            nca = right;
 		            while (left_set.find(nca) == left_set.end()) {
 		              nca = _graph.target((*_matching)[nca]);
 		              nca =(*_blossom_rep)[_blossom_set->
 		                                   find(_graph.target((*_ear)[nca]))];
+		            }
 		          } else {
 		            nca = left;
 		            while (right_set.find(nca) == right_set.end()) {
 		              nca = _graph.target((*_matching)[nca]);
 		              nca = (*_blossom_rep)[_blossom_set->
 		                                   find(_graph.target((*_ear)[nca]))];
+		            }
+		          }
+		        }
+		      }
+		      {
 		        Node node = _graph.u(e);
 		        Arc arc = _graph.direct(e, true);
 		        Node base = (*_blossom_rep)[_blossom_set->find(node)];
 		        while (base != nca) {
-		          _ear->set(node, arc);
+		          (*_ear)[node] = arc;
 		          Node n = node;
 		          while (n != base) {
 		            n = _graph.target((*_matching)[n]);
 		            Arc a = (*_ear)[n];
 		            n = _graph.target(a);
-		            _ear->set(n, _graph.oppositeArc(a));
+		            (*_ear)[n] = _graph.oppositeArc(a);
+		          }
 		          node = _graph.target((*_matching)[base]);
 		          _tree_set->erase(base);
 		          _tree_set->erase(node);
 		          _blossom_set->insert(node, _blossom_set->find(base));
-		          _status->set(node, EVEN);
+		          (*_status)[node] = EVEN;
 		          _node_queue[_last++] = node;
 		          arc = _graph.oppositeArc((*_ear)[node]);
 		          node = _graph.target((*_ear)[node]);
 		          base = (*_blossom_rep)[_blossom_set->find(node)];
 		          _blossom_set->join(_graph.target(arc), base);
+		        }
+		      }
-		      _blossom_rep->set(_blossom_set->find(nca), nca);
+		      (*_blossom_rep)[_blossom_set->find(nca)] = nca;
+		      {
 		        Node node = _graph.v(e);
 		        Arc arc = _graph.direct(e, false);
 		        Node base = (*_blossom_rep)[_blossom_set->find(node)];
 		        while (base != nca) {
-		          _ear->set(node, arc);
+		          (*_ear)[node] = arc;
 		          Node n = node;
 		          while (n != base) {
 		            n = _graph.target((*_matching)[n]);
 		            Arc a = (*_ear)[n];
 		            n = _graph.target(a);
-		            _ear->set(n, _graph.oppositeArc(a));
+		            (*_ear)[n] = _graph.oppositeArc(a);
+		          }
 		          node = _graph.target((*_matching)[base]);
 		          _tree_set->erase(base);
 		          _tree_set->erase(node);
 		          _blossom_set->insert(node, _blossom_set->find(base));
-		          _status->set(node, EVEN);
+		          (*_status)[node] = EVEN;
 		          _node_queue[_last++] = node;
 		          arc = _graph.oppositeArc((*_ear)[node]);
 		          node = _graph.target((*_ear)[node]);
 		          base = (*_blossom_rep)[_blossom_set->find(node)];
 		          _blossom_set->join(_graph.target(arc), base);
+		        }
+		      }
-		      _blossom_rep->set(_blossom_set->find(nca), nca);
+		      (*_blossom_rep)[_blossom_set->find(nca)] = nca;
+		    }
 		    void extendOnArc(const Arc& a) {
 		      Node base = _graph.source(a);
 		      Node odd = _graph.target(a);
-		      _ear->set(odd, _graph.oppositeArc(a));
+		      (*_ear)[odd] = _graph.oppositeArc(a);
 		      Node even = _graph.target((*_matching)[odd]);
 		      _blossom_rep->set(_blossom_set->insert(even), even);
 		      _status->set(odd, ODD);
-		      _status->set(even, EVEN);
+		      (*_blossom_rep)[_blossom_set->insert(even)] = even;
 		      (*_status)[odd] = ODD;
 		      (*_status)[even] = EVEN;
 		      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
 		      _tree_set->insert(odd, tree);
 		      _tree_set->insert(even, tree);
 		      _node_queue[_last++] = even;
+		    }
 		    void augmentOnArc(const Arc& a) {
 		      Node even = _graph.source(a);
 		      Node odd = _graph.target(a);
 		      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
 		      _matching->set(odd, _graph.oppositeArc(a));
 		      _status->set(odd, MATCHED);
 		      (*_matching)[odd] = _graph.oppositeArc(a);
 		      (*_status)[odd] = MATCHED;
 		      Arc arc = (*_matching)[even];
-		      _matching->set(even, a);
+		      (*_matching)[even] = a;
 		      while (arc != INVALID) {
 		        odd = _graph.target(arc);
 		        arc = (*_ear)[odd];
 		        even = _graph.target(arc);
-		        _matching->set(odd, arc);
+		        (*_matching)[odd] = arc;
 		        arc = (*_matching)[even];
-		        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
+		        (*_matching)[even] = _graph.oppositeArc((*_matching)[odd]);
+		      }
 		      for (typename TreeSet::ItemIt it(*_tree_set, tree);
 		           it != INVALID; ++it) {
 		        if ((*_status)[it] == ODD) {
-		          _status->set(it, MATCHED);
+		          (*_status)[it] = MATCHED;
 		        } else {
 		          int blossom = _blossom_set->find(it);
 		          for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
 		               jt != INVALID; ++jt) {
-		            _status->set(jt, MATCHED);
+		            (*_status)[jt] = MATCHED;
+		          }
 		          _blossom_set->eraseClass(blossom);
+		        }
+		      }
 		      _tree_set->eraseClass(tree);
+		    }
 		  public:
 		    /// \brief Constructor
 		    ///
 		    /// Constructor.
 		    MaxMatching(const Graph& graph)
 		      : _graph(graph), _matching(0), _status(0), _ear(0),
 		        _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
 		        _tree_set_index(0), _tree_set(0) {}
 		    ~MaxMatching() {
 		      destroyStructures();
+		    }
 		    /// \name Execution control
 		    /// The simplest way to execute the algorithm is to use the
 		    /// \c run() member function.
 		    /// \n
 		    /// If you need better control on the execution, you must call
 		    /// \ref init(), \ref greedyInit() or \ref matchingInit()
 		    /// functions first, then you can start the algorithm with the \ref
 		    /// startSparse() or startDense() functions.
 		    ///@{
 		    /// \brief Sets the actual matching to the empty matching.
 		    ///
 		    /// Sets the actual matching to the empty matching.
 		    ///
 		    void init() {
 		      createStructures();
 		      for(NodeIt n(_graph); n != INVALID; ++n) {
 		        _matching->set(n, INVALID);
 		        _status->set(n, UNMATCHED);
 		        (*_matching)[n] = INVALID;
 		        (*_status)[n] = UNMATCHED;
+		      }
+		    }
 		    ///\brief Finds an initial matching in a greedy way
 		    ///
 		    ///It finds an initial matching in a greedy way.
 		    void greedyInit() {
 		      createStructures();
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        _matching->set(n, INVALID);
 		        _status->set(n, UNMATCHED);
 		        (*_matching)[n] = INVALID;
 		        (*_status)[n] = UNMATCHED;
+		      }
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_matching)[n] == INVALID) {
 		          for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
 		            Node v = _graph.target(a);
 		            if ((*_matching)[v] == INVALID && v != n) {
 		              _matching->set(n, a);
 		              _status->set(n, MATCHED);
 		              _matching->set(v, _graph.oppositeArc(a));
 		              _status->set(v, MATCHED);
 		              (*_matching)[n] = a;
 		              (*_status)[n] = MATCHED;
 		              (*_matching)[v] = _graph.oppositeArc(a);
 		              (*_status)[v] = MATCHED;
 		              break;
+		            }
+		          }
+		        }
+		      }
+		    }
 		    /// \brief Initialize the matching from a map containing.
 		    ///
 		    /// Initialize the matching from a \c bool valued \c Edge map. This
 		    /// map must have the property that there are no two incident edges
 		    /// with true value, ie. it contains a matching.
 		    /// \return \c true if the map contains a matching.
 		    template <typename MatchingMap>
 		    bool matchingInit(const MatchingMap& matching) {
 		      createStructures();
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        _matching->set(n, INVALID);
 		        _status->set(n, UNMATCHED);
 		        (*_matching)[n] = INVALID;
 		        (*_status)[n] = UNMATCHED;
+		      }
 		      for(EdgeIt e(_graph); e!=INVALID; ++e) {
 		        if (matching[e]) {
 		          Node u = _graph.u(e);
 		          if ((*_matching)[u] != INVALID) return false;
 		          _matching->set(u, _graph.direct(e, true));
 		          _status->set(u, MATCHED);
 		          (*_matching)[u] = _graph.direct(e, true);
 		          (*_status)[u] = MATCHED;
 		          Node v = _graph.v(e);
 		          if ((*_matching)[v] != INVALID) return false;
 		          _matching->set(v, _graph.direct(e, false));
 		          _status->set(v, MATCHED);
 		          (*_matching)[v] = _graph.direct(e, false);
 		          (*_status)[v] = MATCHED;
+		        }
+		      }
 		      return true;
+		    }
 		    /// \brief Starts Edmonds' algorithm
 		    ///
 		    /// If runs the original Edmonds' algorithm.
 		    void startSparse() {
 		      for(NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_status)[n] == UNMATCHED) {
 		          (*_blossom_rep)[_blossom_set->insert(n)] = n;
 		          _tree_set->insert(n);
-		          _status->set(n, EVEN);
+		          (*_status)[n] = EVEN;
 		          processSparse(n);
+		        }
+		      }
+		    }
 		    /// \brief Starts Edmonds' algorithm.
 		    ///
 		    /// It runs Edmonds' algorithm with a heuristic of postponing
 		    /// shrinks, therefore resulting in a faster algorithm for dense graphs.
 		    void startDense() {
 		      for(NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_status)[n] == UNMATCHED) {
 		          (*_blossom_rep)[_blossom_set->insert(n)] = n;
 		          _tree_set->insert(n);
-		          _status->set(n, EVEN);
+		          (*_status)[n] = EVEN;
 		          processDense(n);
+		        }
+		      }
+		    }
 		    /// \brief Runs Edmonds' algorithm
 		    ///
 		    /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>)
 		    /// or Edmonds' algorithm with a heuristic of
 		    /// postponing shrinks for dense graphs.
 		    void run() {
 		      if (countEdges(_graph) < 2 * countNodes(_graph)) {
 		        greedyInit();
 		        startSparse();
 		      } else {
 		        init();
 		        startDense();
+		      }
+		    }
 		    /// @}
 		    /// \name Primal solution
 		    /// Functions to get the primal solution, ie. the matching.
 		    /// @{
 		    ///\brief Returns the size of the current matching.
 		    ///
 		    ///Returns the size of the current matching. After \ref
 		    ///run() it returns the size of the maximum matching in the graph.
 		    int matchingSize() const {
 		      int size = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_matching)[n] != INVALID) {
 		          ++size;
+		        }
+		      }
 		      return size / 2;
+		    }
 		    /// \brief Returns true when the edge is in the matching.
 		    ///
 		    /// Returns true when the edge is in the matching.
 		    bool matching(const Edge& edge) const {
 		      return edge == (*_matching)[_graph.u(edge)];
+		    }
 		    /// \brief Returns the matching edge incident to the given node.
 		    ///
 		    /// Returns the matching edge of a \c node in the actual matching or
 		    /// INVALID if the \c node is not covered by the actual matching.
 		    Arc matching(const Node& n) const {
 		      return (*_matching)[n];
+		    }
 		    ///\brief Returns the mate of a node in the actual matching.
 		    ///
 		    ///Returns the mate of a \c node in the actual matching or
 		    ///INVALID if the \c node is not covered by the actual matching.
 		    Node mate(const Node& n) const {
 		      return (*_matching)[n] != INVALID ?
 		        _graph.target((*_matching)[n]) : INVALID;
+		    }
 		    /// @}
 		    /// \name Dual solution
 		    /// Functions to get the dual solution, ie. the decomposition.
 		    /// @{
 		    /// \brief Returns the class of the node in the Edmonds-Gallai
 		    /// decomposition.
 		    ///
 		    /// Returns the class of the node in the Edmonds-Gallai
 		    /// decomposition.
 		    Status decomposition(const Node& n) const {
 		      return (*_status)[n];
+		    }
 		    /// \brief Returns true when the node is in the barrier.
 		    ///
 		    /// Returns true when the node is in the barrier.
 		    bool barrier(const Node& n) const {
 		      return (*_status)[n] == ODD;
+		    }
 		    /// @}
 		  };
 		  /// \ingroup matching
 		  ///
 		  /// \brief Weighted matching in general graphs
 		  ///
 		  /// This class provides an efficient implementation of Edmond's
 		  /// maximum weighted matching algorithm. The implementation is based
 		  /// on extensive use of priority queues and provides
 		  /// \f$O(nm\log n)\f$ time complexity.
 		  ///
 		  /// The maximum weighted matching problem is to find undirected
 		  /// edges in the graph with maximum overall weight and no two of
 		  /// them shares their ends. The problem can be formulated with the
 		  /// following linear program.
 		  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
 		  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
 		      \quad \forall B\in\mathcal{O}\f] */
 		  /// \f[x_e \ge 0\quad \forall e\in E\f]
 		  /// \f[\max \sum_{e\in E}x_ew_e\f]
 		  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
 		  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
 		  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
 		  /// subsets of the nodes.
 		  ///
 		  /// The algorithm calculates an optimal matching and a proof of the
 		  /// optimality. The solution of the dual problem can be used to check
 		  /// the result of the algorithm. The dual linear problem is the
 		  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
 		      z_B \ge w_{uv} \quad \forall uv\in E\f] */
 		  /// \f[y_u \ge 0 \quad \forall u \in V\f]
 		  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
 		  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
 		      \frac{\vert B \vert - 1}{2}z_B\f] */
 		  ///
 		  /// The algorithm can be executed with \c run() or the \c init() and
 		  /// then the \c start() member functions. After it the matching can
 		  /// be asked with \c matching() or mate() functions. The dual
 		  /// solution can be get with \c nodeValue(), \c blossomNum() and \c
 		  /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
 		  /// "BlossomIt" nested class, which is able to iterate on the nodes
 		  /// of a blossom. If the value type is integral then the dual
 		  /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
 		  template <typename GR,
 		            typename WM = typename GR::template EdgeMap<int> >
 		  class MaxWeightedMatching {
 		  public:
 		    ///\e
 		    typedef GR Graph;
 		    ///\e
 		    typedef WM WeightMap;
 		    ///\e
 		    typedef typename WeightMap::Value Value;
 		    /// \brief Scaling factor for dual solution
 		    ///
 		    /// Scaling factor for dual solution, it is equal to 4 or 1
 		    /// according to the value type.
 		    static const int dualScale =
 		      std::numeric_limits<Value>::is_integer ? 4 : 1;
 		    typedef typename Graph::template NodeMap<typename Graph::Arc>
 		    MatchingMap;
 		  private:
 		    TEMPLATE_GRAPH_TYPEDEFS(Graph);
 		    typedef typename Graph::template NodeMap<Value> NodePotential;
 		    typedef std::vector<Node> BlossomNodeList;
 		    struct BlossomVariable {
 		      int begin, end;
 		      Value value;
 		      BlossomVariable(int _begin, int _end, Value _value)
 		        : begin(_begin), end(_end), value(_value) {}
 		    };
 		    typedef std::vector<BlossomVariable> BlossomPotential;
 		    const Graph& _graph;
 		    const WeightMap& _weight;
 		    MatchingMap* _matching;
 		    NodePotential* _node_potential;
 		    BlossomPotential _blossom_potential;
 		    BlossomNodeList _blossom_node_list;
 		    int _node_num;
 		    int _blossom_num;
 		    typedef RangeMap<int> IntIntMap;
 		    enum Status {
 		      EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
 		    };
@@ -1359,503 +1359,503 @@
 		              nca =
 		                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
 		              right_path.push_back(nca);
 		              nca =
 		                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
 		              right_path.push_back(nca);
+		            }
 		          } else {
 		            nca = left;
 		            while (right_set.find(nca) == right_set.end()) {
 		              nca =
 		                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
 		              left_path.push_back(nca);
 		              nca =
 		                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
 		              left_path.push_back(nca);
+		            }
+		          }
+		        }
+		      }
 		      std::vector<int> subblossoms;
 		      Arc prev;
 		      prev = _graph.direct(edge, true);
 		      for (int i = 0; left_path[i] != nca; i += 2) {
 		        subblossoms.push_back(left_path[i]);
 		        (*_blossom_data)[left_path[i]].next = prev;
 		        _tree_set->erase(left_path[i]);
 		        subblossoms.push_back(left_path[i + 1]);
 		        (*_blossom_data)[left_path[i + 1]].status = EVEN;
 		        oddToEven(left_path[i + 1], tree);
 		        _tree_set->erase(left_path[i + 1]);
 		        prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
+		      }
 		      int k = 0;
 		      while (right_path[k] != nca) ++k;
 		      subblossoms.push_back(nca);
 		      (*_blossom_data)[nca].next = prev;
 		      for (int i = k - 2; i >= 0; i -= 2) {
 		        subblossoms.push_back(right_path[i + 1]);
 		        (*_blossom_data)[right_path[i + 1]].status = EVEN;
 		        oddToEven(right_path[i + 1], tree);
 		        _tree_set->erase(right_path[i + 1]);
 		        (*_blossom_data)[right_path[i + 1]].next =
 		          (*_blossom_data)[right_path[i + 1]].pred;
 		        subblossoms.push_back(right_path[i]);
 		        _tree_set->erase(right_path[i]);
+		      }
 		      int surface =
 		        _blossom_set->join(subblossoms.begin(), subblossoms.end());
 		      for (int i = 0; i < int(subblossoms.size()); ++i) {
 		        if (!_blossom_set->trivial(subblossoms[i])) {
 		          (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
+		        }
 		        (*_blossom_data)[subblossoms[i]].status = MATCHED;
+		      }
 		      (*_blossom_data)[surface].pot = -2 * _delta_sum;
 		      (*_blossom_data)[surface].offset = 0;
 		      (*_blossom_data)[surface].status = EVEN;
 		      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
 		      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
 		      _tree_set->insert(surface, tree);
 		      _tree_set->erase(nca);
+		    }
 		    void splitBlossom(int blossom) {
 		      Arc next = (*_blossom_data)[blossom].next;
 		      Arc pred = (*_blossom_data)[blossom].pred;
 		      int tree = _tree_set->find(blossom);
 		      (*_blossom_data)[blossom].status = MATCHED;
 		      oddToMatched(blossom);
 		      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
 		        _delta2->erase(blossom);
+		      }
 		      std::vector<int> subblossoms;
 		      _blossom_set->split(blossom, std::back_inserter(subblossoms));
 		      Value offset = (*_blossom_data)[blossom].offset;
 		      int b = _blossom_set->find(_graph.source(pred));
 		      int d = _blossom_set->find(_graph.source(next));
 		      int ib = -1, id = -1;
 		      for (int i = 0; i < int(subblossoms.size()); ++i) {
 		        if (subblossoms[i] == b) ib = i;
 		        if (subblossoms[i] == d) id = i;
 		        (*_blossom_data)[subblossoms[i]].offset = offset;
 		        if (!_blossom_set->trivial(subblossoms[i])) {
 		          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
+		        }
 		        if (_blossom_set->classPrio(subblossoms[i]) !=
 		            std::numeric_limits<Value>::max()) {
 		          _delta2->push(subblossoms[i],
 		                        _blossom_set->classPrio(subblossoms[i]) -
 		                        (*_blossom_data)[subblossoms[i]].offset);
+		        }
+		      }
 		      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
 		        for (int i = (id + 1) % subblossoms.size();
 		             i != ib; i = (i + 2) % subblossoms.size()) {
 		          int sb = subblossoms[i];
 		          int tb = subblossoms[(i + 1) % subblossoms.size()];
 		          (*_blossom_data)[sb].next =
 		            _graph.oppositeArc((*_blossom_data)[tb].next);
+		        }
 		        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
 		          int sb = subblossoms[i];
 		          int tb = subblossoms[(i + 1) % subblossoms.size()];
 		          int ub = subblossoms[(i + 2) % subblossoms.size()];
 		          (*_blossom_data)[sb].status = ODD;
 		          matchedToOdd(sb);
 		          _tree_set->insert(sb, tree);
 		          (*_blossom_data)[sb].pred = pred;
 		          (*_blossom_data)[sb].next =
 		                           _graph.oppositeArc((*_blossom_data)[tb].next);
 		          pred = (*_blossom_data)[ub].next;
 		          (*_blossom_data)[tb].status = EVEN;
 		          matchedToEven(tb, tree);
 		          _tree_set->insert(tb, tree);
 		          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
+		        }
 		        (*_blossom_data)[subblossoms[id]].status = ODD;
 		        matchedToOdd(subblossoms[id]);
 		        _tree_set->insert(subblossoms[id], tree);
 		        (*_blossom_data)[subblossoms[id]].next = next;
 		        (*_blossom_data)[subblossoms[id]].pred = pred;
 		      } else {
 		        for (int i = (ib + 1) % subblossoms.size();
 		             i != id; i = (i + 2) % subblossoms.size()) {
 		          int sb = subblossoms[i];
 		          int tb = subblossoms[(i + 1) % subblossoms.size()];
 		          (*_blossom_data)[sb].next =
 		            _graph.oppositeArc((*_blossom_data)[tb].next);
+		        }
 		        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
 		          int sb = subblossoms[i];
 		          int tb = subblossoms[(i + 1) % subblossoms.size()];
 		          int ub = subblossoms[(i + 2) % subblossoms.size()];
 		          (*_blossom_data)[sb].status = ODD;
 		          matchedToOdd(sb);
 		          _tree_set->insert(sb, tree);
 		          (*_blossom_data)[sb].next = next;
 		          (*_blossom_data)[sb].pred =
 		            _graph.oppositeArc((*_blossom_data)[tb].next);
 		          (*_blossom_data)[tb].status = EVEN;
 		          matchedToEven(tb, tree);
 		          _tree_set->insert(tb, tree);
 		          (*_blossom_data)[tb].pred =
 		            (*_blossom_data)[tb].next =
 		            _graph.oppositeArc((*_blossom_data)[ub].next);
 		          next = (*_blossom_data)[ub].next;
+		        }
 		        (*_blossom_data)[subblossoms[ib]].status = ODD;
 		        matchedToOdd(subblossoms[ib]);
 		        _tree_set->insert(subblossoms[ib], tree);
 		        (*_blossom_data)[subblossoms[ib]].next = next;
 		        (*_blossom_data)[subblossoms[ib]].pred = pred;
+		      }
 		      _tree_set->erase(blossom);
+		    }
 		    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
 		      if (_blossom_set->trivial(blossom)) {
 		        int bi = (*_node_index)[base];
 		        Value pot = (*_node_data)[bi].pot;
-		        _matching->set(base, matching);
+		        (*_matching)[base] = matching;
 		        _blossom_node_list.push_back(base);
-		        _node_potential->set(base, pot);
+		        (*_node_potential)[base] = pot;
 		      } else {
 		        Value pot = (*_blossom_data)[blossom].pot;
 		        int bn = _blossom_node_list.size();
 		        std::vector<int> subblossoms;
 		        _blossom_set->split(blossom, std::back_inserter(subblossoms));
 		        int b = _blossom_set->find(base);
 		        int ib = -1;
 		        for (int i = 0; i < int(subblossoms.size()); ++i) {
 		          if (subblossoms[i] == b) { ib = i; break; }
+		        }
 		        for (int i = 1; i < int(subblossoms.size()); i += 2) {
 		          int sb = subblossoms[(ib + i) % subblossoms.size()];
 		          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
 		          Arc m = (*_blossom_data)[tb].next;
 		          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
 		          extractBlossom(tb, _graph.source(m), m);
+		        }
 		        extractBlossom(subblossoms[ib], base, matching);
 		        int en = _blossom_node_list.size();
 		        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
+		      }
+		    }
 		    void extractMatching() {
 		      std::vector<int> blossoms;
 		      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
 		        blossoms.push_back(c);
+		      }
 		      for (int i = 0; i < int(blossoms.size()); ++i) {
 		        if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
 		          Value offset = (*_blossom_data)[blossoms[i]].offset;
 		          (*_blossom_data)[blossoms[i]].pot += 2 * offset;
 		          for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
 		               n != INVALID; ++n) {
 		            (*_node_data)[(*_node_index)[n]].pot -= offset;
+		          }
 		          Arc matching = (*_blossom_data)[blossoms[i]].next;
 		          Node base = _graph.source(matching);
 		          extractBlossom(blossoms[i], base, matching);
 		        } else {
 		          Node base = (*_blossom_data)[blossoms[i]].base;
 		          extractBlossom(blossoms[i], base, INVALID);
+		        }
+		      }
+		    }
 		  public:
 		    /// \brief Constructor
 		    ///
 		    /// Constructor.
 		    MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
 		      : _graph(graph), _weight(weight), _matching(0),
 		        _node_potential(0), _blossom_potential(), _blossom_node_list(),
 		        _node_num(0), _blossom_num(0),
 		        _blossom_index(0), _blossom_set(0), _blossom_data(0),
 		        _node_index(0), _node_heap_index(0), _node_data(0),
 		        _tree_set_index(0), _tree_set(0),
 		        _delta1_index(0), _delta1(0),
 		        _delta2_index(0), _delta2(0),
 		        _delta3_index(0), _delta3(0),
 		        _delta4_index(0), _delta4(0),
 		        _delta_sum() {}
 		    ~MaxWeightedMatching() {
 		      destroyStructures();
+		    }
 		    /// \name Execution control
 		    /// The simplest way to execute the algorithm is to use the
 		    /// \c run() member function.
 		    ///@{
 		    /// \brief Initialize the algorithm
 		    ///
 		    /// Initialize the algorithm
 		    void init() {
 		      createStructures();
 		      for (ArcIt e(_graph); e != INVALID; ++e) {
-		        _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
+		        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
+		      }
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
-		        _delta1_index->set(n, _delta1->PRE_HEAP);
+		        (*_delta1_index)[n] = _delta1->PRE_HEAP;
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
-		        _delta3_index->set(e, _delta3->PRE_HEAP);
+		        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+		      }
 		      for (int i = 0; i < _blossom_num; ++i) {
 		        _delta2_index->set(i, _delta2->PRE_HEAP);
 		        _delta4_index->set(i, _delta4->PRE_HEAP);
 		        (*_delta2_index)[i] = _delta2->PRE_HEAP;
 		        (*_delta4_index)[i] = _delta4->PRE_HEAP;
+		      }
 		      int index = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        Value max = 0;
 		        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		          if (_graph.target(e) == n) continue;
 		          if ((dualScale * _weight[e]) / 2 > max) {
 		            max = (dualScale * _weight[e]) / 2;
+		          }
+		        }
-		        _node_index->set(n, index);
+		        (*_node_index)[n] = index;
 		        (*_node_data)[index].pot = max;
 		        _delta1->push(n, max);
 		        int blossom =
 		          _blossom_set->insert(n, std::numeric_limits<Value>::max());
 		        _tree_set->insert(blossom);
 		        (*_blossom_data)[blossom].status = EVEN;
 		        (*_blossom_data)[blossom].pred = INVALID;
 		        (*_blossom_data)[blossom].next = INVALID;
 		        (*_blossom_data)[blossom].pot = 0;
 		        (*_blossom_data)[blossom].offset = 0;
 		        ++index;
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        int si = (*_node_index)[_graph.u(e)];
 		        int ti = (*_node_index)[_graph.v(e)];
 		        if (_graph.u(e) != _graph.v(e)) {
 		          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
 		                            dualScale * _weight[e]) / 2);
+		        }
+		      }
+		    }
 		    /// \brief Starts the algorithm
 		    ///
 		    /// Starts the algorithm
 		    void start() {
 		      enum OpType {
 		        D1, D2, D3, D4
 		      };
 		      int unmatched = _node_num;
 		      while (unmatched > 0) {
 		        Value d1 = !_delta1->empty() ?
 		          _delta1->prio() : std::numeric_limits<Value>::max();
 		        Value d2 = !_delta2->empty() ?
 		          _delta2->prio() : std::numeric_limits<Value>::max();
 		        Value d3 = !_delta3->empty() ?
 		          _delta3->prio() : std::numeric_limits<Value>::max();
 		        Value d4 = !_delta4->empty() ?
 		          _delta4->prio() : std::numeric_limits<Value>::max();
 		        _delta_sum = d1; OpType ot = D1;
 		        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
 		        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
 		        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
 		        switch (ot) {
 		        case D1:
+		          {
 		            Node n = _delta1->top();
 		            unmatchNode(n);
 		            --unmatched;
+		          }
 		          break;
 		        case D2:
+		          {
 		            int blossom = _delta2->top();
 		            Node n = _blossom_set->classTop(blossom);
 		            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
 		            extendOnArc(e);
+		          }
 		          break;
 		        case D3:
+		          {
 		            Edge e = _delta3->top();
 		            int left_blossom = _blossom_set->find(_graph.u(e));
 		            int right_blossom = _blossom_set->find(_graph.v(e));
 		            if (left_blossom == right_blossom) {
 		              _delta3->pop();
 		            } else {
 		              int left_tree;
 		              if ((*_blossom_data)[left_blossom].status == EVEN) {
 		                left_tree = _tree_set->find(left_blossom);
 		              } else {
 		                left_tree = -1;
 		                ++unmatched;
+		              }
 		              int right_tree;
 		              if ((*_blossom_data)[right_blossom].status == EVEN) {
 		                right_tree = _tree_set->find(right_blossom);
 		              } else {
 		                right_tree = -1;
 		                ++unmatched;
+		              }
 		              if (left_tree == right_tree) {
 		                shrinkOnEdge(e, left_tree);
 		              } else {
 		                augmentOnEdge(e);
 		                unmatched -= 2;
+		              }
+		            }
 		          } break;
 		        case D4:
 		          splitBlossom(_delta4->top());
 		          break;
+		        }
+		      }
 		      extractMatching();
+		    }
 		    /// \brief Runs %MaxWeightedMatching algorithm.
 		    ///
 		    /// This method runs the %MaxWeightedMatching algorithm.
 		    ///
 		    /// \note mwm.run() is just a shortcut of the following code.
 		    /// \code
 		    ///   mwm.init();
 		    ///   mwm.start();
 		    /// \endcode
 		    void run() {
 		      init();
 		      start();
+		    }
 		    /// @}
 		    /// \name Primal solution
 		    /// Functions to get the primal solution, ie. the matching.
 		    /// @{
 		    /// \brief Returns the weight of the matching.
 		    ///
 		    /// Returns the weight of the matching.
 		    Value matchingValue() const {
 		      Value sum = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_matching)[n] != INVALID) {
 		          sum += _weight[(*_matching)[n]];
+		        }
+		      }
 		      return sum /= 2;
+		    }
 		    /// \brief Returns the cardinality of the matching.
 		    ///
 		    /// Returns the cardinality of the matching.
 		    int matchingSize() const {
 		      int num = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_matching)[n] != INVALID) {
 		          ++num;
+		        }
+		      }
 		      return num /= 2;
+		    }
 		    /// \brief Returns true when the edge is in the matching.
 		    ///
 		    /// Returns true when the edge is in the matching.
 		    bool matching(const Edge& edge) const {
 		      return edge == (*_matching)[_graph.u(edge)];
+		    }
 		    /// \brief Returns the incident matching arc.
 		    ///
 		    /// Returns the incident matching arc from given node. If the
 		    /// node is not matched then it gives back \c INVALID.
 		    Arc matching(const Node& node) const {
 		      return (*_matching)[node];
+		    }
 		    /// \brief Returns the mate of the node.
 		    ///
 		    /// Returns the adjancent node in a mathcing arc. If the node is
 		    /// not matched then it gives back \c INVALID.
 		    Node mate(const Node& node) const {
 		      return (*_matching)[node] != INVALID ?
 		        _graph.target((*_matching)[node]) : INVALID;
+		    }
 		    /// @}
 		    /// \name Dual solution
 		    /// Functions to get the dual solution.
 		    /// @{
 		    /// \brief Returns the value of the dual solution.
 		    ///
 		    /// Returns the value of the dual solution. It should be equal to
 		    /// the primal value scaled by \ref dualScale "dual scale".
 		    Value dualValue() const {
@@ -2552,494 +2552,494 @@
 		              nca =
 		                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
 		              right_path.push_back(nca);
 		              nca =
 		                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
 		              right_path.push_back(nca);
+		            }
 		          } else {
 		            nca = left;
 		            while (right_set.find(nca) == right_set.end()) {
 		              nca =
 		                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
 		              left_path.push_back(nca);
 		              nca =
 		                _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
 		              left_path.push_back(nca);
+		            }
+		          }
+		        }
+		      }
 		      std::vector<int> subblossoms;
 		      Arc prev;
 		      prev = _graph.direct(edge, true);
 		      for (int i = 0; left_path[i] != nca; i += 2) {
 		        subblossoms.push_back(left_path[i]);
 		        (*_blossom_data)[left_path[i]].next = prev;
 		        _tree_set->erase(left_path[i]);
 		        subblossoms.push_back(left_path[i + 1]);
 		        (*_blossom_data)[left_path[i + 1]].status = EVEN;
 		        oddToEven(left_path[i + 1], tree);
 		        _tree_set->erase(left_path[i + 1]);
 		        prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
+		      }
 		      int k = 0;
 		      while (right_path[k] != nca) ++k;
 		      subblossoms.push_back(nca);
 		      (*_blossom_data)[nca].next = prev;
 		      for (int i = k - 2; i >= 0; i -= 2) {
 		        subblossoms.push_back(right_path[i + 1]);
 		        (*_blossom_data)[right_path[i + 1]].status = EVEN;
 		        oddToEven(right_path[i + 1], tree);
 		        _tree_set->erase(right_path[i + 1]);
 		        (*_blossom_data)[right_path[i + 1]].next =
 		          (*_blossom_data)[right_path[i + 1]].pred;
 		        subblossoms.push_back(right_path[i]);
 		        _tree_set->erase(right_path[i]);
+		      }
 		      int surface =
 		        _blossom_set->join(subblossoms.begin(), subblossoms.end());
 		      for (int i = 0; i < int(subblossoms.size()); ++i) {
 		        if (!_blossom_set->trivial(subblossoms[i])) {
 		          (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
+		        }
 		        (*_blossom_data)[subblossoms[i]].status = MATCHED;
+		      }
 		      (*_blossom_data)[surface].pot = -2 * _delta_sum;
 		      (*_blossom_data)[surface].offset = 0;
 		      (*_blossom_data)[surface].status = EVEN;
 		      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
 		      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
 		      _tree_set->insert(surface, tree);
 		      _tree_set->erase(nca);
+		    }
 		    void splitBlossom(int blossom) {
 		      Arc next = (*_blossom_data)[blossom].next;
 		      Arc pred = (*_blossom_data)[blossom].pred;
 		      int tree = _tree_set->find(blossom);
 		      (*_blossom_data)[blossom].status = MATCHED;
 		      oddToMatched(blossom);
 		      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
 		        _delta2->erase(blossom);
+		      }
 		      std::vector<int> subblossoms;
 		      _blossom_set->split(blossom, std::back_inserter(subblossoms));
 		      Value offset = (*_blossom_data)[blossom].offset;
 		      int b = _blossom_set->find(_graph.source(pred));
 		      int d = _blossom_set->find(_graph.source(next));
 		      int ib = -1, id = -1;
 		      for (int i = 0; i < int(subblossoms.size()); ++i) {
 		        if (subblossoms[i] == b) ib = i;
 		        if (subblossoms[i] == d) id = i;
 		        (*_blossom_data)[subblossoms[i]].offset = offset;
 		        if (!_blossom_set->trivial(subblossoms[i])) {
 		          (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
+		        }
 		        if (_blossom_set->classPrio(subblossoms[i]) !=
 		            std::numeric_limits<Value>::max()) {
 		          _delta2->push(subblossoms[i],
 		                        _blossom_set->classPrio(subblossoms[i]) -
 		                        (*_blossom_data)[subblossoms[i]].offset);
+		        }
+		      }
 		      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
 		        for (int i = (id + 1) % subblossoms.size();
 		             i != ib; i = (i + 2) % subblossoms.size()) {
 		          int sb = subblossoms[i];
 		          int tb = subblossoms[(i + 1) % subblossoms.size()];
 		          (*_blossom_data)[sb].next =
 		            _graph.oppositeArc((*_blossom_data)[tb].next);
+		        }
 		        for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
 		          int sb = subblossoms[i];
 		          int tb = subblossoms[(i + 1) % subblossoms.size()];
 		          int ub = subblossoms[(i + 2) % subblossoms.size()];
 		          (*_blossom_data)[sb].status = ODD;
 		          matchedToOdd(sb);
 		          _tree_set->insert(sb, tree);
 		          (*_blossom_data)[sb].pred = pred;
 		          (*_blossom_data)[sb].next =
 		                           _graph.oppositeArc((*_blossom_data)[tb].next);
 		          pred = (*_blossom_data)[ub].next;
 		          (*_blossom_data)[tb].status = EVEN;
 		          matchedToEven(tb, tree);
 		          _tree_set->insert(tb, tree);
 		          (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
+		        }
 		        (*_blossom_data)[subblossoms[id]].status = ODD;
 		        matchedToOdd(subblossoms[id]);
 		        _tree_set->insert(subblossoms[id], tree);
 		        (*_blossom_data)[subblossoms[id]].next = next;
 		        (*_blossom_data)[subblossoms[id]].pred = pred;
 		      } else {
 		        for (int i = (ib + 1) % subblossoms.size();
 		             i != id; i = (i + 2) % subblossoms.size()) {
 		          int sb = subblossoms[i];
 		          int tb = subblossoms[(i + 1) % subblossoms.size()];
 		          (*_blossom_data)[sb].next =
 		            _graph.oppositeArc((*_blossom_data)[tb].next);
+		        }
 		        for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
 		          int sb = subblossoms[i];
 		          int tb = subblossoms[(i + 1) % subblossoms.size()];
 		          int ub = subblossoms[(i + 2) % subblossoms.size()];
 		          (*_blossom_data)[sb].status = ODD;
 		          matchedToOdd(sb);
 		          _tree_set->insert(sb, tree);
 		          (*_blossom_data)[sb].next = next;
 		          (*_blossom_data)[sb].pred =
 		            _graph.oppositeArc((*_blossom_data)[tb].next);
 		          (*_blossom_data)[tb].status = EVEN;
 		          matchedToEven(tb, tree);
 		          _tree_set->insert(tb, tree);
 		          (*_blossom_data)[tb].pred =
 		            (*_blossom_data)[tb].next =
 		            _graph.oppositeArc((*_blossom_data)[ub].next);
 		          next = (*_blossom_data)[ub].next;
+		        }
 		        (*_blossom_data)[subblossoms[ib]].status = ODD;
 		        matchedToOdd(subblossoms[ib]);
 		        _tree_set->insert(subblossoms[ib], tree);
 		        (*_blossom_data)[subblossoms[ib]].next = next;
 		        (*_blossom_data)[subblossoms[ib]].pred = pred;
+		      }
 		      _tree_set->erase(blossom);
+		    }
 		    void extractBlossom(int blossom, const Node& base, const Arc& matching) {
 		      if (_blossom_set->trivial(blossom)) {
 		        int bi = (*_node_index)[base];
 		        Value pot = (*_node_data)[bi].pot;
-		        _matching->set(base, matching);
+		        (*_matching)[base] = matching;
 		        _blossom_node_list.push_back(base);
-		        _node_potential->set(base, pot);
+		        (*_node_potential)[base] = pot;
 		      } else {
 		        Value pot = (*_blossom_data)[blossom].pot;
 		        int bn = _blossom_node_list.size();
 		        std::vector<int> subblossoms;
 		        _blossom_set->split(blossom, std::back_inserter(subblossoms));
 		        int b = _blossom_set->find(base);
 		        int ib = -1;
 		        for (int i = 0; i < int(subblossoms.size()); ++i) {
 		          if (subblossoms[i] == b) { ib = i; break; }
+		        }
 		        for (int i = 1; i < int(subblossoms.size()); i += 2) {
 		          int sb = subblossoms[(ib + i) % subblossoms.size()];
 		          int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
 		          Arc m = (*_blossom_data)[tb].next;
 		          extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
 		          extractBlossom(tb, _graph.source(m), m);
+		        }
 		        extractBlossom(subblossoms[ib], base, matching);
 		        int en = _blossom_node_list.size();
 		        _blossom_potential.push_back(BlossomVariable(bn, en, pot));
+		      }
+		    }
 		    void extractMatching() {
 		      std::vector<int> blossoms;
 		      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
 		        blossoms.push_back(c);
+		      }
 		      for (int i = 0; i < int(blossoms.size()); ++i) {
 		        Value offset = (*_blossom_data)[blossoms[i]].offset;
 		        (*_blossom_data)[blossoms[i]].pot += 2 * offset;
 		        for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
 		             n != INVALID; ++n) {
 		          (*_node_data)[(*_node_index)[n]].pot -= offset;
+		        }
 		        Arc matching = (*_blossom_data)[blossoms[i]].next;
 		        Node base = _graph.source(matching);
 		        extractBlossom(blossoms[i], base, matching);
+		      }
+		    }
 		  public:
 		    /// \brief Constructor
 		    ///
 		    /// Constructor.
 		    MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
 		      : _graph(graph), _weight(weight), _matching(0),
 		        _node_potential(0), _blossom_potential(), _blossom_node_list(),
 		        _node_num(0), _blossom_num(0),
 		        _blossom_index(0), _blossom_set(0), _blossom_data(0),
 		        _node_index(0), _node_heap_index(0), _node_data(0),
 		        _tree_set_index(0), _tree_set(0),
 		        _delta2_index(0), _delta2(0),
 		        _delta3_index(0), _delta3(0),
 		        _delta4_index(0), _delta4(0),
 		        _delta_sum() {}
 		    ~MaxWeightedPerfectMatching() {
 		      destroyStructures();
+		    }
 		    /// \name Execution control
 		    /// The simplest way to execute the algorithm is to use the
 		    /// \c run() member function.
 		    ///@{
 		    /// \brief Initialize the algorithm
 		    ///
 		    /// Initialize the algorithm
 		    void init() {
 		      createStructures();
 		      for (ArcIt e(_graph); e != INVALID; ++e) {
-		        _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
+		        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
-		        _delta3_index->set(e, _delta3->PRE_HEAP);
+		        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+		      }
 		      for (int i = 0; i < _blossom_num; ++i) {
 		        _delta2_index->set(i, _delta2->PRE_HEAP);
 		        _delta4_index->set(i, _delta4->PRE_HEAP);
 		        (*_delta2_index)[i] = _delta2->PRE_HEAP;
 		        (*_delta4_index)[i] = _delta4->PRE_HEAP;
+		      }
 		      int index = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        Value max = - std::numeric_limits<Value>::max();
 		        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		          if (_graph.target(e) == n) continue;
 		          if ((dualScale * _weight[e]) / 2 > max) {
 		            max = (dualScale * _weight[e]) / 2;
+		          }
+		        }
-		        _node_index->set(n, index);
+		        (*_node_index)[n] = index;
 		        (*_node_data)[index].pot = max;
 		        int blossom =
 		          _blossom_set->insert(n, std::numeric_limits<Value>::max());
 		        _tree_set->insert(blossom);
 		        (*_blossom_data)[blossom].status = EVEN;
 		        (*_blossom_data)[blossom].pred = INVALID;
 		        (*_blossom_data)[blossom].next = INVALID;
 		        (*_blossom_data)[blossom].pot = 0;
 		        (*_blossom_data)[blossom].offset = 0;
 		        ++index;
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        int si = (*_node_index)[_graph.u(e)];
 		        int ti = (*_node_index)[_graph.v(e)];
 		        if (_graph.u(e) != _graph.v(e)) {
 		          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
 		                            dualScale * _weight[e]) / 2);
+		        }
+		      }
+		    }
 		    /// \brief Starts the algorithm
 		    ///
 		    /// Starts the algorithm
 		    bool start() {
 		      enum OpType {
 		        D2, D3, D4
 		      };
 		      int unmatched = _node_num;
 		      while (unmatched > 0) {
 		        Value d2 = !_delta2->empty() ?
 		          _delta2->prio() : std::numeric_limits<Value>::max();
 		        Value d3 = !_delta3->empty() ?
 		          _delta3->prio() : std::numeric_limits<Value>::max();
 		        Value d4 = !_delta4->empty() ?
 		          _delta4->prio() : std::numeric_limits<Value>::max();
 		        _delta_sum = d2; OpType ot = D2;
 		        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
 		        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
 		        if (_delta_sum == std::numeric_limits<Value>::max()) {
 		          return false;
+		        }
 		        switch (ot) {
 		        case D2:
+		          {
 		            int blossom = _delta2->top();
 		            Node n = _blossom_set->classTop(blossom);
 		            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
 		            extendOnArc(e);
+		          }
 		          break;
 		        case D3:
+		          {
 		            Edge e = _delta3->top();
 		            int left_blossom = _blossom_set->find(_graph.u(e));
 		            int right_blossom = _blossom_set->find(_graph.v(e));
 		            if (left_blossom == right_blossom) {
 		              _delta3->pop();
 		            } else {
 		              int left_tree = _tree_set->find(left_blossom);
 		              int right_tree = _tree_set->find(right_blossom);
 		              if (left_tree == right_tree) {
 		                shrinkOnEdge(e, left_tree);
 		              } else {
 		                augmentOnEdge(e);
 		                unmatched -= 2;
+		              }
+		            }
 		          } break;
 		        case D4:
 		          splitBlossom(_delta4->top());
 		          break;
+		        }
+		      }
 		      extractMatching();
 		      return true;
+		    }
 		    /// \brief Runs %MaxWeightedPerfectMatching algorithm.
 		    ///
 		    /// This method runs the %MaxWeightedPerfectMatching algorithm.
 		    ///
 		    /// \note mwm.run() is just a shortcut of the following code.
 		    /// \code
 		    ///   mwm.init();
 		    ///   mwm.start();
 		    /// \endcode
 		    bool run() {
 		      init();
 		      return start();
+		    }
 		    /// @}
 		    /// \name Primal solution
 		    /// Functions to get the primal solution, ie. the matching.
 		    /// @{
 		    /// \brief Returns the matching value.
 		    ///
 		    /// Returns the matching value.
 		    Value matchingValue() const {
 		      Value sum = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if ((*_matching)[n] != INVALID) {
 		          sum += _weight[(*_matching)[n]];
+		        }
+		      }
 		      return sum /= 2;
+		    }
 		    /// \brief Returns true when the edge is in the matching.
 		    ///
 		    /// Returns true when the edge is in the matching.
 		    bool matching(const Edge& edge) const {
 		      return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
+		    }
 		    /// \brief Returns the incident matching edge.
 		    ///
 		    /// Returns the incident matching arc from given edge.
 		    Arc matching(const Node& node) const {
 		      return (*_matching)[node];
+		    }
 		    /// \brief Returns the mate of the node.
 		    ///
 		    /// Returns the adjancent node in a mathcing arc.
 		    Node mate(const Node& node) const {
 		      return _graph.target((*_matching)[node]);
+		    }
 		    /// @}
 		    /// \name Dual solution
 		    /// Functions to get the dual solution.
 		    /// @{
 		    /// \brief Returns the value of the dual solution.
 		    ///
 		    /// Returns the value of the dual solution. It should be equal to
 		    /// the primal value scaled by \ref dualScale "dual scale".
 		    Value dualValue() const {
 		      Value sum = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        sum += nodeValue(n);
+		      }
 		      for (int i = 0; i < blossomNum(); ++i) {
 		        sum += blossomValue(i) * (blossomSize(i) / 2);
+		      }
 		      return sum;
+		    }
 		    /// \brief Returns the value of the node.
 		    ///
 		    /// Returns the the value of the node.
 		    Value nodeValue(const Node& n) const {
 		      return (*_node_potential)[n];
+		    }
 		    /// \brief Returns the number of the blossoms in the basis.
 		    ///
 		    /// Returns the number of the blossoms in the basis.
 		    /// \see BlossomIt
 		    int blossomNum() const {
 		      return _blossom_potential.size();
+		    }
 		    /// \brief Returns the number of the nodes in the blossom.
 		    ///
 		    /// Returns the number of the nodes in the blossom.
 		    int blossomSize(int k) const {
 		      return _blossom_potential[k].end - _blossom_potential[k].begin;
+		    }
 		    /// \brief Returns the value of the blossom.
 		    ///
 		    /// Returns the the value of the blossom.

lemon/min_cost_arborescence.h

Show inline comments

@@ -104,691 +104,691 @@
 		  /// the optimality of the solution can be checked.
 		  ///
 		  /// \param GR The digraph type the algorithm runs on. The default value
 		  /// is \ref ListDigraph.
 		  /// \param CM This read-only ArcMap determines the costs of the
 		  /// arcs. It is read once for each arc, so the map may involve in
 		  /// relatively time consuming process to compute the arc cost if
 		  /// it is necessary. The default map type is \ref
 		  /// concepts::Digraph::ArcMap "Digraph::ArcMap<int>".
 		  /// \param TR Traits class to set various data types used
 		  /// by the algorithm. The default traits class is
 		  /// \ref MinCostArborescenceDefaultTraits
 		  /// "MinCostArborescenceDefaultTraits<GR, CM>".  See \ref
 		  /// MinCostArborescenceDefaultTraits for the documentation of a
 		  /// MinCostArborescence traits class.
 		#ifndef DOXYGEN
 		  template <typename GR = ListDigraph,
 		            typename CM = typename GR::template ArcMap<int>,
 		            typename TR =
 		              MinCostArborescenceDefaultTraits<GR, CM> >
 		#else
 		  template <typename GR, typename CM, typedef TR>
 		#endif
 		  class MinCostArborescence {
 		  public:
 		    /// The traits.
 		    typedef TR Traits;
 		    /// The type of the underlying digraph.
 		    typedef typename Traits::Digraph Digraph;
 		    /// The type of the map that stores the arc costs.
 		    typedef typename Traits::CostMap CostMap;
 		    ///The type of the costs of the arcs.
 		    typedef typename Traits::Value Value;
 		    ///The type of the predecessor map.
 		    typedef typename Traits::PredMap PredMap;
 		    ///The type of the map that stores which arcs are in the arborescence.
 		    typedef typename Traits::ArborescenceMap ArborescenceMap;
 		    typedef MinCostArborescence Create;
 		  private:
 		    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
 		    struct CostArc {
 		      Arc arc;
 		      Value value;
 		      CostArc() {}
 		      CostArc(Arc _arc, Value _value) : arc(_arc), value(_value) {}
 		    };
 		    const Digraph *_digraph;
 		    const CostMap *_cost;
 		    PredMap *_pred;
 		    bool local_pred;
 		    ArborescenceMap *_arborescence;
 		    bool local_arborescence;
 		    typedef typename Digraph::template ArcMap<int> ArcOrder;
 		    ArcOrder *_arc_order;
 		    typedef typename Digraph::template NodeMap<int> NodeOrder;
 		    NodeOrder *_node_order;
 		    typedef typename Digraph::template NodeMap<CostArc> CostArcMap;
 		    CostArcMap *_cost_arcs;
 		    struct StackLevel {
 		      std::vector<CostArc> arcs;
 		      int node_level;
 		    };
 		    std::vector<StackLevel> level_stack;
 		    std::vector<Node> queue;
 		    typedef std::vector<typename Digraph::Node> DualNodeList;
 		    DualNodeList _dual_node_list;
 		    struct DualVariable {
 		      int begin, end;
 		      Value value;
 		      DualVariable(int _begin, int _end, Value _value)
 		        : begin(_begin), end(_end), value(_value) {}
 		    };
 		    typedef std::vector<DualVariable> DualVariables;
 		    DualVariables _dual_variables;
 		    typedef typename Digraph::template NodeMap<int> HeapCrossRef;
 		    HeapCrossRef *_heap_cross_ref;
 		    typedef BinHeap<int, HeapCrossRef> Heap;
 		    Heap *_heap;
 		  protected:
 		    MinCostArborescence() {}
 		  private:
 		    void createStructures() {
 		      if (!_pred) {
 		        local_pred = true;
 		        _pred = Traits::createPredMap(*_digraph);
+		      }
 		      if (!_arborescence) {
 		        local_arborescence = true;
 		        _arborescence = Traits::createArborescenceMap(*_digraph);
+		      }
 		      if (!_arc_order) {
 		        _arc_order = new ArcOrder(*_digraph);
+		      }
 		      if (!_node_order) {
 		        _node_order = new NodeOrder(*_digraph);
+		      }
 		      if (!_cost_arcs) {
 		        _cost_arcs = new CostArcMap(*_digraph);
+		      }
 		      if (!_heap_cross_ref) {
 		        _heap_cross_ref = new HeapCrossRef(*_digraph, -1);
+		      }
 		      if (!_heap) {
 		        _heap = new Heap(*_heap_cross_ref);
+		      }
+		    }
 		    void destroyStructures() {
 		      if (local_arborescence) {
 		        delete _arborescence;
+		      }
 		      if (local_pred) {
 		        delete _pred;
+		      }
 		      if (_arc_order) {
 		        delete _arc_order;
+		      }
 		      if (_node_order) {
 		        delete _node_order;
+		      }
 		      if (_cost_arcs) {
 		        delete _cost_arcs;
+		      }
 		      if (_heap) {
 		        delete _heap;
+		      }
 		      if (_heap_cross_ref) {
 		        delete _heap_cross_ref;
+		      }
+		    }
 		    Arc prepare(Node node) {
 		      std::vector<Node> nodes;
 		      (*_node_order)[node] = _dual_node_list.size();
 		      StackLevel level;
 		      level.node_level = _dual_node_list.size();
 		      _dual_node_list.push_back(node);
 		      for (InArcIt it(*_digraph, node); it != INVALID; ++it) {
 		        Arc arc = it;
 		        Node source = _digraph->source(arc);
 		        Value value = (*_cost)[it];
 		        if (source == node || (*_node_order)[source] == -3) continue;
 		        if ((*_cost_arcs)[source].arc == INVALID) {
 		          (*_cost_arcs)[source].arc = arc;
 		          (*_cost_arcs)[source].value = value;
 		          nodes.push_back(source);
 		        } else {
 		          if ((*_cost_arcs)[source].value > value) {
 		            (*_cost_arcs)[source].arc = arc;
 		            (*_cost_arcs)[source].value = value;
+		          }
+		        }
+		      }
 		      CostArc minimum = (*_cost_arcs)[nodes[0]];
 		      for (int i = 1; i < int(nodes.size()); ++i) {
 		        if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
 		          minimum = (*_cost_arcs)[nodes[i]];
+		        }
+		      }
-		      _arc_order->set(minimum.arc, _dual_variables.size());
+		      (*_arc_order)[minimum.arc] = _dual_variables.size();
 		      DualVariable var(_dual_node_list.size() - 1,
 		                       _dual_node_list.size(), minimum.value);
 		      _dual_variables.push_back(var);
 		      for (int i = 0; i < int(nodes.size()); ++i) {
 		        (*_cost_arcs)[nodes[i]].value -= minimum.value;
 		        level.arcs.push_back((*_cost_arcs)[nodes[i]]);
 		        (*_cost_arcs)[nodes[i]].arc = INVALID;
+		      }
 		      level_stack.push_back(level);
 		      return minimum.arc;
+		    }
 		    Arc contract(Node node) {
 		      int node_bottom = bottom(node);
 		      std::vector<Node> nodes;
 		      while (!level_stack.empty() &&
 		             level_stack.back().node_level >= node_bottom) {
 		        for (int i = 0; i < int(level_stack.back().arcs.size()); ++i) {
 		          Arc arc = level_stack.back().arcs[i].arc;
 		          Node source = _digraph->source(arc);
 		          Value value = level_stack.back().arcs[i].value;
 		          if ((*_node_order)[source] >= node_bottom) continue;
 		          if ((*_cost_arcs)[source].arc == INVALID) {
 		            (*_cost_arcs)[source].arc = arc;
 		            (*_cost_arcs)[source].value = value;
 		            nodes.push_back(source);
 		          } else {
 		            if ((*_cost_arcs)[source].value > value) {
 		              (*_cost_arcs)[source].arc = arc;
 		              (*_cost_arcs)[source].value = value;
+		            }
+		          }
+		        }
 		        level_stack.pop_back();
+		      }
 		      CostArc minimum = (*_cost_arcs)[nodes[0]];
 		      for (int i = 1; i < int(nodes.size()); ++i) {
 		        if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
 		          minimum = (*_cost_arcs)[nodes[i]];
+		        }
+		      }
-		      _arc_order->set(minimum.arc, _dual_variables.size());
+		      (*_arc_order)[minimum.arc] = _dual_variables.size();
 		      DualVariable var(node_bottom, _dual_node_list.size(), minimum.value);
 		      _dual_variables.push_back(var);
 		      StackLevel level;
 		      level.node_level = node_bottom;
 		      for (int i = 0; i < int(nodes.size()); ++i) {
 		        (*_cost_arcs)[nodes[i]].value -= minimum.value;
 		        level.arcs.push_back((*_cost_arcs)[nodes[i]]);
 		        (*_cost_arcs)[nodes[i]].arc = INVALID;
+		      }
 		      level_stack.push_back(level);
 		      return minimum.arc;
+		    }
 		    int bottom(Node node) {
 		      int k = level_stack.size() - 1;
 		      while (level_stack[k].node_level > (*_node_order)[node]) {
 		        --k;
+		      }
 		      return level_stack[k].node_level;
+		    }
 		    void finalize(Arc arc) {
 		      Node node = _digraph->target(arc);
 		      _heap->push(node, (*_arc_order)[arc]);
 		      _pred->set(node, arc);
 		      while (!_heap->empty()) {
 		        Node source = _heap->top();
 		        _heap->pop();
-		        _node_order->set(source, -1);
+		        (*_node_order)[source] = -1;
 		        for (OutArcIt it(*_digraph, source); it != INVALID; ++it) {
 		          if ((*_arc_order)[it] < 0) continue;
 		          Node target = _digraph->target(it);
 		          switch(_heap->state(target)) {
 		          case Heap::PRE_HEAP:
 		            _heap->push(target, (*_arc_order)[it]);
 		            _pred->set(target, it);
 		            break;
 		          case Heap::IN_HEAP:
 		            if ((*_arc_order)[it] < (*_heap)[target]) {
 		              _heap->decrease(target, (*_arc_order)[it]);
 		              _pred->set(target, it);
+		            }
 		            break;
 		          case Heap::POST_HEAP:
 		            break;
+		          }
+		        }
 		        _arborescence->set((*_pred)[source], true);
+		      }
+		    }
 		  public:
 		    /// \name Named template parameters
 		    /// @{
 		    template <class T>
 		    struct DefArborescenceMapTraits : public Traits {
 		      typedef T ArborescenceMap;
 		      static ArborescenceMap *createArborescenceMap(const Digraph &)
+		      {
 		        LEMON_ASSERT(false, "ArborescenceMap is not initialized");
 		        return 0; // ignore warnings
+		      }
 		    };
 		    /// \brief \ref named-templ-param "Named parameter" for
 		    /// setting ArborescenceMap type
 		    ///
 		    /// \ref named-templ-param "Named parameter" for setting
 		    /// ArborescenceMap type
 		    template <class T>
 		    struct DefArborescenceMap
 		      : public MinCostArborescence<Digraph, CostMap,
 		                                   DefArborescenceMapTraits<T> > {
 		    };
 		    template <class T>
 		    struct DefPredMapTraits : public Traits {
 		      typedef T PredMap;
 		      static PredMap *createPredMap(const Digraph &)
+		      {
 		        LEMON_ASSERT(false, "PredMap is not initialized");
+		      }
 		    };
 		    /// \brief \ref named-templ-param "Named parameter" for
 		    /// setting PredMap type
 		    ///
 		    /// \ref named-templ-param "Named parameter" for setting
 		    /// PredMap type
 		    template <class T>
 		    struct DefPredMap
 		      : public MinCostArborescence<Digraph, CostMap, DefPredMapTraits<T> > {
 		    };
 		    /// @}
 		    /// \brief Constructor.
 		    ///
 		    /// \param digraph The digraph the algorithm will run on.
 		    /// \param cost The cost map used by the algorithm.
 		    MinCostArborescence(const Digraph& digraph, const CostMap& cost)
 		      : _digraph(&digraph), _cost(&cost), _pred(0), local_pred(false),
 		        _arborescence(0), local_arborescence(false),
 		        _arc_order(0), _node_order(0), _cost_arcs(0),
 		        _heap_cross_ref(0), _heap(0) {}
 		    /// \brief Destructor.
 		    ~MinCostArborescence() {
 		      destroyStructures();
+		    }
 		    /// \brief Sets the arborescence map.
 		    ///
 		    /// Sets the arborescence map.
 		    /// \return <tt>(*this)</tt>
 		    MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
 		      if (local_arborescence) {
 		        delete _arborescence;
+		      }
 		      local_arborescence = false;
 		      _arborescence = &m;
 		      return *this;
+		    }
 		    /// \brief Sets the arborescence map.
 		    ///
 		    /// Sets the arborescence map.
 		    /// \return <tt>(*this)</tt>
 		    MinCostArborescence& predMap(PredMap& m) {
 		      if (local_pred) {
 		        delete _pred;
+		      }
 		      local_pred = false;
 		      _pred = &m;
 		      return *this;
+		    }
 		    /// \name Query Functions
 		    /// The result of the %MinCostArborescence algorithm can be obtained
 		    /// using these functions.\n
 		    /// Before the use of these functions,
 		    /// either run() or start() must be called.
 		    /// @{
 		    /// \brief Returns a reference to the arborescence map.
 		    ///
 		    /// Returns a reference to the arborescence map.
 		    const ArborescenceMap& arborescenceMap() const {
 		      return *_arborescence;
+		    }
 		    /// \brief Returns true if the arc is in the arborescence.
 		    ///
 		    /// Returns true if the arc is in the arborescence.
 		    /// \param arc The arc of the digraph.
 		    /// \pre \ref run() must be called before using this function.
 		    bool arborescence(Arc arc) const {
 		      return (*_pred)[_digraph->target(arc)] == arc;
+		    }
 		    /// \brief Returns a reference to the pred map.
 		    ///
 		    /// Returns a reference to the pred map.
 		    const PredMap& predMap() const {
 		      return *_pred;
+		    }
 		    /// \brief Returns the predecessor arc of the given node.
 		    ///
 		    /// Returns the predecessor arc of the given node.
 		    Arc pred(Node node) const {
 		      return (*_pred)[node];
+		    }
 		    /// \brief Returns the cost of the arborescence.
 		    ///
 		    /// Returns the cost of the arborescence.
 		    Value arborescenceValue() const {
 		      Value sum = 0;
 		      for (ArcIt it(*_digraph); it != INVALID; ++it) {
 		        if (arborescence(it)) {
 		          sum += (*_cost)[it];
+		        }
+		      }
 		      return sum;
+		    }
 		    /// \brief Indicates that a node is reachable from the sources.
 		    ///
 		    /// Indicates that a node is reachable from the sources.
 		    bool reached(Node node) const {
 		      return (*_node_order)[node] != -3;
+		    }
 		    /// \brief Indicates that a node is processed.
 		    ///
 		    /// Indicates that a node is processed. The arborescence path exists
 		    /// from the source to the given node.
 		    bool processed(Node node) const {
 		      return (*_node_order)[node] == -1;
+		    }
 		    /// \brief Returns the number of the dual variables in basis.
 		    ///
 		    /// Returns the number of the dual variables in basis.
 		    int dualNum() const {
 		      return _dual_variables.size();
+		    }
 		    /// \brief Returns the value of the dual solution.
 		    ///
 		    /// Returns the value of the dual solution. It should be
 		    /// equal to the arborescence value.
 		    Value dualValue() const {
 		      Value sum = 0;
 		      for (int i = 0; i < int(_dual_variables.size()); ++i) {
 		        sum += _dual_variables[i].value;
+		      }
 		      return sum;
+		    }
 		    /// \brief Returns the number of the nodes in the dual variable.
 		    ///
 		    /// Returns the number of the nodes in the dual variable.
 		    int dualSize(int k) const {
 		      return _dual_variables[k].end - _dual_variables[k].begin;
+		    }
 		    /// \brief Returns the value of the dual variable.
 		    ///
 		    /// Returns the the value of the dual variable.
 		    const Value& dualValue(int k) const {
 		      return _dual_variables[k].value;
+		    }
 		    /// \brief Lemon iterator for get a dual variable.
 		    ///
 		    /// Lemon iterator for get a dual variable. This class provides
 		    /// a common style lemon iterator which gives back a subset of
 		    /// the nodes.
 		    class DualIt {
 		    public:
 		      /// \brief Constructor.
 		      ///
 		      /// Constructor for get the nodeset of the variable.
 		      DualIt(const MinCostArborescence& algorithm, int variable)
 		        : _algorithm(&algorithm)
+		      {
 		        _index = _algorithm->_dual_variables[variable].begin;
 		        _last = _algorithm->_dual_variables[variable].end;
+		      }
 		      /// \brief Conversion to node.
 		      ///
 		      /// Conversion to node.
 		      operator Node() const {
 		        return _algorithm->_dual_node_list[_index];
+		      }
 		      /// \brief Increment operator.
 		      ///
 		      /// Increment operator.
 		      DualIt& operator++() {
 		        ++_index;
 		        return *this;
+		      }
 		      /// \brief Validity checking
 		      ///
 		      /// Checks whether the iterator is invalid.
 		      bool operator==(Invalid) const {
 		        return _index == _last;
+		      }
 		      /// \brief Validity checking
 		      ///
 		      /// Checks whether the iterator is valid.
 		      bool operator!=(Invalid) const {
 		        return _index != _last;
+		      }
 		    private:
 		      const MinCostArborescence* _algorithm;
 		      int _index, _last;
 		    };
 		    /// @}
 		    /// \name Execution control
 		    /// The simplest way to execute the algorithm is to use
 		    /// one of the member functions called \c run(...). \n
 		    /// If you need more control on the execution,
 		    /// first you must call \ref init(), then you can add several
 		    /// source nodes with \ref addSource().
 		    /// Finally \ref start() will perform the arborescence
 		    /// computation.
 		    ///@{
 		    /// \brief Initializes the internal data structures.
 		    ///
 		    /// Initializes the internal data structures.
 		    ///
 		    void init() {
 		      createStructures();
 		      _heap->clear();
 		      for (NodeIt it(*_digraph); it != INVALID; ++it) {
 		        (*_cost_arcs)[it].arc = INVALID;
 		        _node_order->set(it, -3);
 		        _heap_cross_ref->set(it, Heap::PRE_HEAP);
 		        (*_node_order)[it] = -3;
 		        (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
 		        _pred->set(it, INVALID);
+		      }
 		      for (ArcIt it(*_digraph); it != INVALID; ++it) {
 		        _arborescence->set(it, false);
-		        _arc_order->set(it, -1);
+		        (*_arc_order)[it] = -1;
+		      }
 		      _dual_node_list.clear();
 		      _dual_variables.clear();
+		    }
 		    /// \brief Adds a new source node.
 		    ///
 		    /// Adds a new source node to the algorithm.
 		    void addSource(Node source) {
 		      std::vector<Node> nodes;
 		      nodes.push_back(source);
 		      while (!nodes.empty()) {
 		        Node node = nodes.back();
 		        nodes.pop_back();
 		        for (OutArcIt it(*_digraph, node); it != INVALID; ++it) {
 		          Node target = _digraph->target(it);
 		          if ((*_node_order)[target] == -3) {
 		            (*_node_order)[target] = -2;
 		            nodes.push_back(target);
 		            queue.push_back(target);
+		          }
+		        }
+		      }
 		      (*_node_order)[source] = -1;
+		    }
 		    /// \brief Processes the next node in the priority queue.
 		    ///
 		    /// Processes the next node in the priority queue.
 		    ///
 		    /// \return The processed node.
 		    ///
 		    /// \warning The queue must not be empty!
 		    Node processNextNode() {
 		      Node node = queue.back();
 		      queue.pop_back();
 		      if ((*_node_order)[node] == -2) {
 		        Arc arc = prepare(node);
 		        Node source = _digraph->source(arc);
 		        while ((*_node_order)[source] != -1) {
 		          if ((*_node_order)[source] >= 0) {
 		            arc = contract(source);
 		          } else {
 		            arc = prepare(source);
+		          }
 		          source = _digraph->source(arc);
+		        }
 		        finalize(arc);
 		        level_stack.clear();
+		      }
 		      return node;
+		    }
 		    /// \brief Returns the number of the nodes to be processed.
 		    ///
 		    /// Returns the number of the nodes to be processed.
 		    int queueSize() const {
 		      return queue.size();
+		    }
 		    /// \brief Returns \c false if there are nodes to be processed.
 		    ///
 		    /// Returns \c false if there are nodes to be processed.
 		    bool emptyQueue() const {
 		      return queue.empty();
+		    }
 		    /// \brief Executes the algorithm.
 		    ///
 		    /// Executes the algorithm.
 		    ///
 		    /// \pre init() must be called and at least one node should be added
 		    /// with addSource() before using this function.
 		    ///
 		    ///\note mca.start() is just a shortcut of the following code.
 		    ///\code
 		    ///while (!mca.emptyQueue()) {
 		    ///  mca.processNextNode();
 		    ///}
 		    ///\endcode
 		    void start() {
 		      while (!emptyQueue()) {
 		        processNextNode();
+		      }
+		    }
 		    /// \brief Runs %MinCostArborescence algorithm from node \c s.
 		    ///
 		    /// This method runs the %MinCostArborescence algorithm from
 		    /// a root node \c s.
 		    ///
 		    /// \note mca.run(s) is just a shortcut of the following code.
 		    /// \code
 		    /// mca.init();
 		    /// mca.addSource(s);
 		    /// mca.start();
 		    /// \endcode
 		    void run(Node node) {
 		      init();
 		      addSource(node);
 		      start();
+		    }
 		    ///@}
 		  };
 		  /// \ingroup spantree
 		  ///
 		  /// \brief Function type interface for MinCostArborescence algorithm.
 		  ///
 		  /// Function type interface for MinCostArborescence algorithm.
 		  /// \param digraph The Digraph that the algorithm runs on.
 		  /// \param cost The CostMap of the arcs.
 		  /// \param source The source of the arborescence.
 		  /// \retval arborescence The bool ArcMap which stores the arborescence.
 		  /// \return The cost of the arborescence.
 		  ///
 		  /// \sa MinCostArborescence
 		  template <typename Digraph, typename CostMap, typename ArborescenceMap>
 		  typename CostMap::Value minCostArborescence(const Digraph& digraph,
 		                                              const CostMap& cost,
 		                                              typename Digraph::Node source,
 		                                              ArborescenceMap& arborescence) {
 		    typename MinCostArborescence<Digraph, CostMap>
 		      ::template DefArborescenceMap<ArborescenceMap>
 		      ::Create mca(digraph, cost);
 		    mca.arborescenceMap(arborescence);
 		    mca.run(source);
 		    return mca.arborescenceValue();
+		  }
+		}
 		#endif

lemon/preflow.h

Show inline comments

@@ -215,751 +215,751 @@
 		      typedef Preflow<Digraph, CapacityMap,
 		                      SetFlowMapTraits<T> > Create;
 		    };
 		    template <typename T>
 		    struct SetElevatorTraits : public Traits {
 		      typedef T Elevator;
 		      static Elevator *createElevator(const Digraph&, int) {
 		        LEMON_ASSERT(false, "Elevator is not initialized");
 		        return 0; // ignore warnings
+		      }
 		    };
 		    /// \brief \ref named-templ-param "Named parameter" for setting
 		    /// Elevator type
 		    ///
 		    /// \ref named-templ-param "Named parameter" for setting Elevator
 		    /// type. If this named parameter is used, then an external
 		    /// elevator object must be passed to the algorithm using the
 		    /// \ref elevator(Elevator&) "elevator()" function before calling
 		    /// \ref run() or \ref init().
 		    /// \sa SetStandardElevator
 		    template <typename T>
 		    struct SetElevator
 		      : public Preflow<Digraph, CapacityMap, SetElevatorTraits<T> > {
 		      typedef Preflow<Digraph, CapacityMap,
 		                      SetElevatorTraits<T> > Create;
 		    };
 		    template <typename T>
 		    struct SetStandardElevatorTraits : public Traits {
 		      typedef T Elevator;
 		      static Elevator *createElevator(const Digraph& digraph, int max_level) {
 		        return new Elevator(digraph, max_level);
+		      }
 		    };
 		    /// \brief \ref named-templ-param "Named parameter" for setting
 		    /// Elevator type with automatic allocation
 		    ///
 		    /// \ref named-templ-param "Named parameter" for setting Elevator
 		    /// type with automatic allocation.
 		    /// The Elevator should have standard constructor interface to be
 		    /// able to automatically created by the algorithm (i.e. the
 		    /// digraph and the maximum level should be passed to it).
 		    /// However an external elevator object could also be passed to the
 		    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
 		    /// before calling \ref run() or \ref init().
 		    /// \sa SetElevator
 		    template <typename T>
 		    struct SetStandardElevator
 		      : public Preflow<Digraph, CapacityMap,
 		                       SetStandardElevatorTraits<T> > {
 		      typedef Preflow<Digraph, CapacityMap,
 		                      SetStandardElevatorTraits<T> > Create;
 		    };
 		    /// @}
 		  protected:
 		    Preflow() {}
 		  public:
 		    /// \brief The constructor of the class.
 		    ///
 		    /// The constructor of the class.
 		    /// \param digraph The digraph the algorithm runs on.
 		    /// \param capacity The capacity of the arcs.
 		    /// \param source The source node.
 		    /// \param target The target node.
 		    Preflow(const Digraph& digraph, const CapacityMap& capacity,
 		            Node source, Node target)
 		      : _graph(digraph), _capacity(&capacity),
 		        _node_num(0), _source(source), _target(target),
 		        _flow(0), _local_flow(false),
 		        _level(0), _local_level(false),
 		        _excess(0), _tolerance(), _phase() {}
 		    /// \brief Destructor.
 		    ///
 		    /// Destructor.
 		    ~Preflow() {
 		      destroyStructures();
+		    }
 		    /// \brief Sets the capacity map.
 		    ///
 		    /// Sets the capacity map.
 		    /// \return <tt>(*this)</tt>
 		    Preflow& capacityMap(const CapacityMap& map) {
 		      _capacity = &map;
 		      return *this;
+		    }
 		    /// \brief Sets the flow map.
 		    ///
 		    /// Sets the flow map.
 		    /// If you don't use this function before calling \ref run() or
 		    /// \ref init(), an instance will be allocated automatically.
 		    /// The destructor deallocates this automatically allocated map,
 		    /// of course.
 		    /// \return <tt>(*this)</tt>
 		    Preflow& flowMap(FlowMap& map) {
 		      if (_local_flow) {
 		        delete _flow;
 		        _local_flow = false;
+		      }
 		      _flow = &map;
 		      return *this;
+		    }
 		    /// \brief Sets the source node.
 		    ///
 		    /// Sets the source node.
 		    /// \return <tt>(*this)</tt>
 		    Preflow& source(const Node& node) {
 		      _source = node;
 		      return *this;
+		    }
 		    /// \brief Sets the target node.
 		    ///
 		    /// Sets the target node.
 		    /// \return <tt>(*this)</tt>
 		    Preflow& target(const Node& node) {
 		      _target = node;
 		      return *this;
+		    }
 		    /// \brief Sets the elevator used by algorithm.
 		    ///
 		    /// Sets the elevator used by algorithm.
 		    /// If you don't use this function before calling \ref run() or
 		    /// \ref init(), an instance will be allocated automatically.
 		    /// The destructor deallocates this automatically allocated elevator,
 		    /// of course.
 		    /// \return <tt>(*this)</tt>
 		    Preflow& elevator(Elevator& elevator) {
 		      if (_local_level) {
 		        delete _level;
 		        _local_level = false;
+		      }
 		      _level = &elevator;
 		      return *this;
+		    }
 		    /// \brief Returns a const reference to the elevator.
 		    ///
 		    /// Returns a const reference to the elevator.
 		    ///
 		    /// \pre Either \ref run() or \ref init() must be called before
 		    /// using this function.
 		    const Elevator& elevator() const {
 		      return *_level;
+		    }
 		    /// \brief Sets the tolerance used by algorithm.
 		    ///
 		    /// Sets the tolerance used by algorithm.
 		    Preflow& tolerance(const Tolerance& tolerance) const {
 		      _tolerance = tolerance;
 		      return *this;
+		    }
 		    /// \brief Returns a const reference to the tolerance.
 		    ///
 		    /// Returns a const reference to the tolerance.
 		    const Tolerance& tolerance() const {
 		      return tolerance;
+		    }
 		    /// \name Execution Control
 		    /// The simplest way to execute the preflow algorithm is to use
 		    /// \ref run() or \ref runMinCut().\n
 		    /// If you need more control on the initial solution or the execution,
 		    /// first you have to call one of the \ref init() functions, then
 		    /// \ref startFirstPhase() and if you need it \ref startSecondPhase().
 		    ///@{
 		    /// \brief Initializes the internal data structures.
 		    ///
 		    /// Initializes the internal data structures and sets the initial
 		    /// flow to zero on each arc.
 		    void init() {
 		      createStructures();
 		      _phase = true;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
-		        _excess->set(n, 0);
+		        (*_excess)[n] = 0;
+		      }
 		      for (ArcIt e(_graph); e != INVALID; ++e) {
 		        _flow->set(e, 0);
+		      }
 		      typename Digraph::template NodeMap<bool> reached(_graph, false);
 		      _level->initStart();
 		      _level->initAddItem(_target);
 		      std::vector<Node> queue;
-		      reached.set(_source, true);
+		      reached[_source] = true;
 		      queue.push_back(_target);
-		      reached.set(_target, true);
+		      reached[_target] = true;
 		      while (!queue.empty()) {
 		        _level->initNewLevel();
 		        std::vector<Node> nqueue;
 		        for (int i = 0; i < int(queue.size()); ++i) {
 		          Node n = queue[i];
 		          for (InArcIt e(_graph, n); e != INVALID; ++e) {
 		            Node u = _graph.source(e);
 		            if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
-		              reached.set(u, true);
+		              reached[u] = true;
 		              _level->initAddItem(u);
 		              nqueue.push_back(u);
+		            }
+		          }
+		        }
 		        queue.swap(nqueue);
+		      }
 		      _level->initFinish();
 		      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
 		        if (_tolerance.positive((*_capacity)[e])) {
 		          Node u = _graph.target(e);
 		          if ((*_level)[u] == _level->maxLevel()) continue;
 		          _flow->set(e, (*_capacity)[e]);
-		          _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
+		          (*_excess)[u] += (*_capacity)[e];
 		          if (u != _target && !_level->active(u)) {
 		            _level->activate(u);
+		          }
+		        }
+		      }
+		    }
 		    /// \brief Initializes the internal data structures using the
 		    /// given flow map.
 		    ///
 		    /// Initializes the internal data structures and sets the initial
 		    /// flow to the given \c flowMap. The \c flowMap should contain a
 		    /// flow or at least a preflow, i.e. at each node excluding the
 		    /// source node the incoming flow should greater or equal to the
 		    /// outgoing flow.
 		    /// \return \c false if the given \c flowMap is not a preflow.
 		    template <typename FlowMap>
 		    bool init(const FlowMap& flowMap) {
 		      createStructures();
 		      for (ArcIt e(_graph); e != INVALID; ++e) {
 		        _flow->set(e, flowMap[e]);
+		      }
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        Value excess = 0;
 		        for (InArcIt e(_graph, n); e != INVALID; ++e) {
 		          excess += (*_flow)[e];
+		        }
 		        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		          excess -= (*_flow)[e];
+		        }
 		        if (excess < 0 && n != _source) return false;
-		        _excess->set(n, excess);
+		        (*_excess)[n] = excess;
+		      }
 		      typename Digraph::template NodeMap<bool> reached(_graph, false);
 		      _level->initStart();
 		      _level->initAddItem(_target);
 		      std::vector<Node> queue;
-		      reached.set(_source, true);
+		      reached[_source] = true;
 		      queue.push_back(_target);
-		      reached.set(_target, true);
+		      reached[_target] = true;
 		      while (!queue.empty()) {
 		        _level->initNewLevel();
 		        std::vector<Node> nqueue;
 		        for (int i = 0; i < int(queue.size()); ++i) {
 		          Node n = queue[i];
 		          for (InArcIt e(_graph, n); e != INVALID; ++e) {
 		            Node u = _graph.source(e);
 		            if (!reached[u] &&
 		                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
-		              reached.set(u, true);
+		              reached[u] = true;
 		              _level->initAddItem(u);
 		              nqueue.push_back(u);
+		            }
+		          }
 		          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		            Node v = _graph.target(e);
 		            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
-		              reached.set(v, true);
+		              reached[v] = true;
 		              _level->initAddItem(v);
 		              nqueue.push_back(v);
+		            }
+		          }
+		        }
 		        queue.swap(nqueue);
+		      }
 		      _level->initFinish();
 		      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
 		        Value rem = (*_capacity)[e] - (*_flow)[e];
 		        if (_tolerance.positive(rem)) {
 		          Node u = _graph.target(e);
 		          if ((*_level)[u] == _level->maxLevel()) continue;
 		          _flow->set(e, (*_capacity)[e]);
-		          _excess->set(u, (*_excess)[u] + rem);
+		          (*_excess)[u] += rem;
 		          if (u != _target && !_level->active(u)) {
 		            _level->activate(u);
+		          }
+		        }
+		      }
 		      for (InArcIt e(_graph, _source); e != INVALID; ++e) {
 		        Value rem = (*_flow)[e];
 		        if (_tolerance.positive(rem)) {
 		          Node v = _graph.source(e);
 		          if ((*_level)[v] == _level->maxLevel()) continue;
 		          _flow->set(e, 0);
-		          _excess->set(v, (*_excess)[v] + rem);
+		          (*_excess)[v] += rem;
 		          if (v != _target && !_level->active(v)) {
 		            _level->activate(v);
+		          }
+		        }
+		      }
 		      return true;
+		    }
 		    /// \brief Starts the first phase of the preflow algorithm.
 		    ///
 		    /// The preflow algorithm consists of two phases, this method runs
 		    /// the first phase. After the first phase the maximum flow value
 		    /// and a minimum value cut can already be computed, although a
 		    /// maximum flow is not yet obtained. So after calling this method
 		    /// \ref flowValue() returns the value of a maximum flow and \ref
 		    /// minCut() returns a minimum cut.
 		    /// \pre One of the \ref init() functions must be called before
 		    /// using this function.
 		    void startFirstPhase() {
 		      _phase = true;
 		      Node n = _level->highestActive();
 		      int level = _level->highestActiveLevel();
 		      while (n != INVALID) {
 		        int num = _node_num;
 		        while (num > 0 && n != INVALID) {
 		          Value excess = (*_excess)[n];
 		          int new_level = _level->maxLevel();
 		          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		            Value rem = (*_capacity)[e] - (*_flow)[e];
 		            if (!_tolerance.positive(rem)) continue;
 		            Node v = _graph.target(e);
 		            if ((*_level)[v] < level) {
 		              if (!_level->active(v) && v != _target) {
 		                _level->activate(v);
+		              }
 		              if (!_tolerance.less(rem, excess)) {
 		                _flow->set(e, (*_flow)[e] + excess);
-		                _excess->set(v, (*_excess)[v] + excess);
+		                (*_excess)[v] += excess;
 		                excess = 0;
 		                goto no_more_push_1;
 		              } else {
 		                excess -= rem;
-		                _excess->set(v, (*_excess)[v] + rem);
+		                (*_excess)[v] += rem;
 		                _flow->set(e, (*_capacity)[e]);
+		              }
 		            } else if (new_level > (*_level)[v]) {
 		              new_level = (*_level)[v];
+		            }
+		          }
 		          for (InArcIt e(_graph, n); e != INVALID; ++e) {
 		            Value rem = (*_flow)[e];
 		            if (!_tolerance.positive(rem)) continue;
 		            Node v = _graph.source(e);
 		            if ((*_level)[v] < level) {
 		              if (!_level->active(v) && v != _target) {
 		                _level->activate(v);
+		              }
 		              if (!_tolerance.less(rem, excess)) {
 		                _flow->set(e, (*_flow)[e] - excess);
-		                _excess->set(v, (*_excess)[v] + excess);
+		                (*_excess)[v] += excess;
 		                excess = 0;
 		                goto no_more_push_1;
 		              } else {
 		                excess -= rem;
-		                _excess->set(v, (*_excess)[v] + rem);
+		                (*_excess)[v] += rem;
 		                _flow->set(e, 0);
+		              }
 		            } else if (new_level > (*_level)[v]) {
 		              new_level = (*_level)[v];
+		            }
+		          }
 		        no_more_push_1:
-		          _excess->set(n, excess);
+		          (*_excess)[n] = excess;
 		          if (excess != 0) {
 		            if (new_level + 1 < _level->maxLevel()) {
 		              _level->liftHighestActive(new_level + 1);
 		            } else {
 		              _level->liftHighestActiveToTop();
+		            }
 		            if (_level->emptyLevel(level)) {
 		              _level->liftToTop(level);
+		            }
 		          } else {
 		            _level->deactivate(n);
+		          }
 		          n = _level->highestActive();
 		          level = _level->highestActiveLevel();
 		          --num;
+		        }
 		        num = _node_num * 20;
 		        while (num > 0 && n != INVALID) {
 		          Value excess = (*_excess)[n];
 		          int new_level = _level->maxLevel();
 		          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		            Value rem = (*_capacity)[e] - (*_flow)[e];
 		            if (!_tolerance.positive(rem)) continue;
 		            Node v = _graph.target(e);
 		            if ((*_level)[v] < level) {
 		              if (!_level->active(v) && v != _target) {
 		                _level->activate(v);
+		              }
 		              if (!_tolerance.less(rem, excess)) {
 		                _flow->set(e, (*_flow)[e] + excess);
-		                _excess->set(v, (*_excess)[v] + excess);
+		                (*_excess)[v] += excess;
 		                excess = 0;
 		                goto no_more_push_2;
 		              } else {
 		                excess -= rem;
-		                _excess->set(v, (*_excess)[v] + rem);
+		                (*_excess)[v] += rem;
 		                _flow->set(e, (*_capacity)[e]);
+		              }
 		            } else if (new_level > (*_level)[v]) {
 		              new_level = (*_level)[v];
+		            }
+		          }
 		          for (InArcIt e(_graph, n); e != INVALID; ++e) {
 		            Value rem = (*_flow)[e];
 		            if (!_tolerance.positive(rem)) continue;
 		            Node v = _graph.source(e);
 		            if ((*_level)[v] < level) {
 		              if (!_level->active(v) && v != _target) {
 		                _level->activate(v);
+		              }
 		              if (!_tolerance.less(rem, excess)) {
 		                _flow->set(e, (*_flow)[e] - excess);
-		                _excess->set(v, (*_excess)[v] + excess);
+		                (*_excess)[v] += excess;
 		                excess = 0;
 		                goto no_more_push_2;
 		              } else {
 		                excess -= rem;
-		                _excess->set(v, (*_excess)[v] + rem);
+		                (*_excess)[v] += rem;
 		                _flow->set(e, 0);
+		              }
 		            } else if (new_level > (*_level)[v]) {
 		              new_level = (*_level)[v];
+		            }
+		          }
 		        no_more_push_2:
-		          _excess->set(n, excess);
+		          (*_excess)[n] = excess;
 		          if (excess != 0) {
 		            if (new_level + 1 < _level->maxLevel()) {
 		              _level->liftActiveOn(level, new_level + 1);
 		            } else {
 		              _level->liftActiveToTop(level);
+		            }
 		            if (_level->emptyLevel(level)) {
 		              _level->liftToTop(level);
+		            }
 		          } else {
 		            _level->deactivate(n);
+		          }
 		          while (level >= 0 && _level->activeFree(level)) {
 		            --level;
+		          }
 		          if (level == -1) {
 		            n = _level->highestActive();
 		            level = _level->highestActiveLevel();
 		          } else {
 		            n = _level->activeOn(level);
+		          }
 		          --num;
+		        }
+		      }
+		    }
 		    /// \brief Starts the second phase of the preflow algorithm.
 		    ///
 		    /// The preflow algorithm consists of two phases, this method runs
 		    /// the second phase. After calling one of the \ref init() functions
 		    /// and \ref startFirstPhase() and then \ref startSecondPhase(),
 		    /// \ref flowMap() returns a maximum flow, \ref flowValue() returns the
 		    /// value of a maximum flow, \ref minCut() returns a minimum cut
 		    /// \pre One of the \ref init() functions and \ref startFirstPhase()
 		    /// must be called before using this function.
 		    void startSecondPhase() {
 		      _phase = false;
 		      typename Digraph::template NodeMap<bool> reached(_graph);
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
-		        reached.set(n, (*_level)[n] < _level->maxLevel());
+		        reached[n] = (*_level)[n] < _level->maxLevel();
+		      }
 		      _level->initStart();
 		      _level->initAddItem(_source);
 		      std::vector<Node> queue;
 		      queue.push_back(_source);
-		      reached.set(_source, true);
+		      reached[_source] = true;
 		      while (!queue.empty()) {
 		        _level->initNewLevel();
 		        std::vector<Node> nqueue;
 		        for (int i = 0; i < int(queue.size()); ++i) {
 		          Node n = queue[i];
 		          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		            Node v = _graph.target(e);
 		            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
-		              reached.set(v, true);
+		              reached[v] = true;
 		              _level->initAddItem(v);
 		              nqueue.push_back(v);
+		            }
+		          }
 		          for (InArcIt e(_graph, n); e != INVALID; ++e) {
 		            Node u = _graph.source(e);
 		            if (!reached[u] &&
 		                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
-		              reached.set(u, true);
+		              reached[u] = true;
 		              _level->initAddItem(u);
 		              nqueue.push_back(u);
+		            }
+		          }
+		        }
 		        queue.swap(nqueue);
+		      }
 		      _level->initFinish();
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if (!reached[n]) {
 		          _level->dirtyTopButOne(n);
 		        } else if ((*_excess)[n] > 0 && _target != n) {
 		          _level->activate(n);
+		        }
+		      }
 		      Node n;
 		      while ((n = _level->highestActive()) != INVALID) {
 		        Value excess = (*_excess)[n];
 		        int level = _level->highestActiveLevel();
 		        int new_level = _level->maxLevel();
 		        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		          Value rem = (*_capacity)[e] - (*_flow)[e];
 		          if (!_tolerance.positive(rem)) continue;
 		          Node v = _graph.target(e);
 		          if ((*_level)[v] < level) {
 		            if (!_level->active(v) && v != _source) {
 		              _level->activate(v);
+		            }
 		            if (!_tolerance.less(rem, excess)) {
 		              _flow->set(e, (*_flow)[e] + excess);
-		              _excess->set(v, (*_excess)[v] + excess);
+		              (*_excess)[v] += excess;
 		              excess = 0;
 		              goto no_more_push;
 		            } else {
 		              excess -= rem;
-		              _excess->set(v, (*_excess)[v] + rem);
+		              (*_excess)[v] += rem;
 		              _flow->set(e, (*_capacity)[e]);
+		            }
 		          } else if (new_level > (*_level)[v]) {
 		            new_level = (*_level)[v];
+		          }
+		        }
 		        for (InArcIt e(_graph, n); e != INVALID; ++e) {
 		          Value rem = (*_flow)[e];
 		          if (!_tolerance.positive(rem)) continue;
 		          Node v = _graph.source(e);
 		          if ((*_level)[v] < level) {
 		            if (!_level->active(v) && v != _source) {
 		              _level->activate(v);
+		            }
 		            if (!_tolerance.less(rem, excess)) {
 		              _flow->set(e, (*_flow)[e] - excess);
-		              _excess->set(v, (*_excess)[v] + excess);
+		              (*_excess)[v] += excess;
 		              excess = 0;
 		              goto no_more_push;
 		            } else {
 		              excess -= rem;
-		              _excess->set(v, (*_excess)[v] + rem);
+		              (*_excess)[v] += rem;
 		              _flow->set(e, 0);
+		            }
 		          } else if (new_level > (*_level)[v]) {
 		            new_level = (*_level)[v];
+		          }
+		        }
 		      no_more_push:
-		        _excess->set(n, excess);
+		        (*_excess)[n] = excess;
 		        if (excess != 0) {
 		          if (new_level + 1 < _level->maxLevel()) {
 		            _level->liftHighestActive(new_level + 1);
 		          } else {
 		            // Calculation error
 		            _level->liftHighestActiveToTop();
+		          }
 		          if (_level->emptyLevel(level)) {
 		            // Calculation error
 		            _level->liftToTop(level);
+		          }
 		        } else {
 		          _level->deactivate(n);
+		        }
+		      }
+		    }
 		    /// \brief Runs the preflow algorithm.
 		    ///
 		    /// Runs the preflow algorithm.
 		    /// \note pf.run() is just a shortcut of the following code.
 		    /// \code
 		    ///   pf.init();
 		    ///   pf.startFirstPhase();
 		    ///   pf.startSecondPhase();
 		    /// \endcode
 		    void run() {
 		      init();
 		      startFirstPhase();
 		      startSecondPhase();
+		    }
 		    /// \brief Runs the preflow algorithm to compute the minimum cut.
 		    ///
 		    /// Runs the preflow algorithm to compute the minimum cut.
 		    /// \note pf.runMinCut() is just a shortcut of the following code.
 		    /// \code
 		    ///   pf.init();
 		    ///   pf.startFirstPhase();
 		    /// \endcode
 		    void runMinCut() {
 		      init();
 		      startFirstPhase();
+		    }
 		    /// @}
 		    /// \name Query Functions
 		    /// The results of the preflow algorithm can be obtained using these
 		    /// functions.\n
 		    /// Either one of the \ref run() "run*()" functions or one of the
 		    /// \ref startFirstPhase() "start*()" functions should be called
 		    /// before using them.
 		    ///@{
 		    /// \brief Returns the value of the maximum flow.
 		    ///
 		    /// Returns the value of the maximum flow by returning the excess
 		    /// of the target node. This value equals to the value of
 		    /// the maximum flow already after the first phase of the algorithm.
 		    ///
 		    /// \pre Either \ref run() or \ref init() must be called before
 		    /// using this function.
 		    Value flowValue() const {
 		      return (*_excess)[_target];
+		    }
 		    /// \brief Returns the flow on the given arc.
 		    ///
 		    /// Returns the flow on the given arc. This method can
 		    /// be called after the second phase of the algorithm.
 		    ///
 		    /// \pre Either \ref run() or \ref init() must be called before
 		    /// using this function.
 		    Value flow(const Arc& arc) const {
 		      return (*_flow)[arc];
+		    }
 		    /// \brief Returns a const reference to the flow map.
 		    ///
 		    /// Returns a const reference to the arc map storing the found flow.
 		    /// This method can be called after the second phase of the algorithm.
 		    ///
 		    /// \pre Either \ref run() or \ref init() must be called before
 		    /// using this function.
 		    const FlowMap& flowMap() const {
 		      return *_flow;
+		    }
 		    /// \brief Returns \c true when the node is on the source side of the
 		    /// minimum cut.
 		    ///
 		    /// Returns true when the node is on the source side of the found
 		    /// minimum cut. This method can be called both after running \ref
 		    /// startFirstPhase() and \ref startSecondPhase().
 		    ///
 		    /// \pre Either \ref run() or \ref init() must be called before
 		    /// using this function.
 		    bool minCut(const Node& node) const {
 		      return ((*_level)[node] == _level->maxLevel()) == _phase;
+		    }
 		    /// \brief Gives back a minimum value cut.
 		    ///
 		    /// Sets \c cutMap to the characteristic vector of a minimum value
 		    /// cut. \c cutMap should be a \ref concepts::WriteMap "writable"
 		    /// node map with \c bool (or convertible) value type.
 		    ///
 		    /// This method can be called both after running \ref startFirstPhase()
 		    /// and \ref startSecondPhase(). The result after the second phase
 		    /// could be slightly different if inexact computation is used.
 		    ///
 		    /// \note This function calls \ref minCut() for each node, so it runs in
 		    /// O(n) time.
 		    ///
 		    /// \pre Either \ref run() or \ref init() must be called before
 		    /// using this function.
 		    template <typename CutMap>
 		    void minCutMap(CutMap& cutMap) const {
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        cutMap.set(n, minCut(n));
+		      }
+		    }
 		    /// @}
 		  };
+		}
 		#endif

test/kruskal_test.cc

Show inline comments

 		/* -*- mode: C++; indent-tabs-mode: nil; -*-
+		 *
 		 * This file is a part of LEMON, a generic C++ optimization library.
+		 *
 		 * Copyright (C) 2003-2009
 		 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
 		 * (Egervary Research Group on Combinatorial Optimization, EGRES).
+		 *
 		 * Permission to use, modify and distribute this software is granted
 		 * provided that this copyright notice appears in all copies. For
 		 * precise terms see the accompanying LICENSE file.
+		 *
 		 * This software is provided "AS IS" with no warranty of any kind,
 		 * express or implied, and with no claim as to its suitability for any
 		 * purpose.
+		 *
 		 */
 		#include <iostream>
 		#include <vector>
 		#include "test_tools.h"
 		#include <lemon/maps.h>
 		#include <lemon/kruskal.h>
 		#include <lemon/list_graph.h>
 		#include <lemon/concepts/maps.h>
 		#include <lemon/concepts/digraph.h>
 		#include <lemon/concepts/graph.h>
 		using namespace std;
 		using namespace lemon;
 		void checkCompileKruskal()
+		{
 		  concepts::WriteMap<concepts::Digraph::Arc,bool> w;
 		  concepts::WriteMap<concepts::Graph::Edge,bool> uw;
 		  concepts::ReadMap<concepts::Digraph::Arc,int> r;
 		  concepts::ReadMap<concepts::Graph::Edge,int> ur;
 		  concepts::Digraph g;
 		  concepts::Graph ug;
 		  kruskal(g, r, w);
 		  kruskal(ug, ur, uw);
 		  std::vector<std::pair<concepts::Digraph::Arc, int> > rs;
 		  std::vector<std::pair<concepts::Graph::Edge, int> > urs;
 		  kruskal(g, rs, w);
 		  kruskal(ug, urs, uw);
 		  std::vector<concepts::Digraph::Arc> ws;
 		  std::vector<concepts::Graph::Edge> uws;
 		  kruskal(g, r, ws.begin());
 		  kruskal(ug, ur, uws.begin());
+		}
 		int main() {
 		  typedef ListGraph::Node Node;
 		  typedef ListGraph::Edge Edge;
 		  typedef ListGraph::NodeIt NodeIt;
 		  typedef ListGraph::ArcIt ArcIt;
 		  ListGraph G;
 		  Node s=G.addNode();
 		  Node v1=G.addNode();
 		  Node v2=G.addNode();
 		  Node v3=G.addNode();
 		  Node v4=G.addNode();
 		  Node t=G.addNode();
 		  Edge e1 = G.addEdge(s, v1);
 		  Edge e2 = G.addEdge(s, v2);
 		  Edge e3 = G.addEdge(v1, v2);
 		  Edge e4 = G.addEdge(v2, v1);
 		  Edge e5 = G.addEdge(v1, v3);
 		  Edge e6 = G.addEdge(v3, v2);
 		  Edge e7 = G.addEdge(v2, v4);
 		  Edge e8 = G.addEdge(v4, v3);
 		  Edge e9 = G.addEdge(v3, t);
 		  Edge e10 = G.addEdge(v4, t);
 		  typedef ListGraph::EdgeMap<int> ECostMap;
 		  typedef ListGraph::EdgeMap<bool> EBoolMap;
 		  ECostMap edge_cost_map(G, 2);
 		  EBoolMap tree_map(G);
 		  //Test with const map.
 		  check(kruskal(G, ConstMap<ListGraph::Edge,int>(2), tree_map)==10,
 		        "Total cost should be 10");
 		  //Test with an edge map (filled with uniform costs).
 		  check(kruskal(G, edge_cost_map, tree_map)==10,
 		        "Total cost should be 10");
 		  edge_cost_map.set(e1, -10);
 		  edge_cost_map.set(e2, -9);
 		  edge_cost_map.set(e3, -8);
 		  edge_cost_map.set(e4, -7);
 		  edge_cost_map.set(e5, -6);
 		  edge_cost_map.set(e6, -5);
 		  edge_cost_map.set(e7, -4);
 		  edge_cost_map.set(e8, -3);
 		  edge_cost_map.set(e9, -2);
 		  edge_cost_map.set(e10, -1);
 		  edge_cost_map[e1] = -10;
 		  edge_cost_map[e2] = -9;
 		  edge_cost_map[e3] = -8;
 		  edge_cost_map[e4] = -7;
 		  edge_cost_map[e5] = -6;
 		  edge_cost_map[e6] = -5;
 		  edge_cost_map[e7] = -4;
 		  edge_cost_map[e8] = -3;
 		  edge_cost_map[e9] = -2;
 		  edge_cost_map[e10] = -1;
 		  vector<Edge> tree_edge_vec(5);
 		  //Test with a edge map and inserter.
 		  check(kruskal(G, edge_cost_map,
 		                 tree_edge_vec.begin())
 		        ==-31,
 		        "Total cost should be -31.");
 		  tree_edge_vec.clear();
 		  check(kruskal(G, edge_cost_map,
 		                back_inserter(tree_edge_vec))
 		        ==-31,
 		        "Total cost should be -31.");
 		//   tree_edge_vec.clear();
 		//   //The above test could also be coded like this:
 		//   check(kruskal(G,
 		//                 makeKruskalMapInput(G, edge_cost_map),
 		//                 makeKruskalSequenceOutput(back_inserter(tree_edge_vec)))
 		//         ==-31,
 		//         "Total cost should be -31.");
 		  check(tree_edge_vec.size()==5,"The tree should have 5 edges.");
 		  check(tree_edge_vec[0]==e1 &&
 		        tree_edge_vec[1]==e2 &&
 		        tree_edge_vec[2]==e5 &&
 		        tree_edge_vec[3]==e7 &&
 		        tree_edge_vec[4]==e9,
 		        "Wrong tree.");
 		  return 0;
+		}

0 comments (0 inline)

RhodeCode

Login to your account