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:

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

Changeset was too big and was cut off... Show full diff

↑ Collapse diff ↑

lemon/circulation.h

Show inline comments

@@ -72,682 +72,682 @@
 		    /// The type of the map that stores the flow values.
 		    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
 		    typedef typename Digraph::template ArcMap<Value> FlowMap;
 		    /// \brief Instantiates a FlowMap.
 		    ///
 		    /// This function instantiates a \ref FlowMap.
 		    /// \param digraph The digraph, to which we would like to define
 		    /// the flow map.
 		    static FlowMap* createFlowMap(const Digraph& digraph) {
 		      return new FlowMap(digraph);
+		    }
 		    /// \brief The elevator type used by the algorithm.
 		    ///
 		    /// The elevator type used by the algorithm.
 		    ///
 		    /// \sa Elevator
 		    /// \sa LinkedElevator
 		    typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
 		    /// \brief Instantiates an Elevator.
 		    ///
 		    /// This function instantiates an \ref Elevator.
 		    /// \param digraph The digraph, to which we would like to define
 		    /// the elevator.
 		    /// \param max_level The maximum level of the elevator.
 		    static Elevator* createElevator(const Digraph& digraph, int max_level) {
 		      return new Elevator(digraph, max_level);
+		    }
 		    /// \brief The tolerance used by the algorithm
 		    ///
 		    /// The tolerance used by the algorithm to handle inexact computation.
 		    typedef lemon::Tolerance<Value> Tolerance;
 		  };
 		  /**
 		     \brief Push-relabel algorithm for the network circulation problem.
 		     \ingroup max_flow
 		     This class implements a push-relabel algorithm for the network
 		     circulation problem.
 		     It is to find a feasible circulation when lower and upper bounds
 		     are given for the flow values on the arcs and lower bounds
 		     are given for the supply values of the nodes.
 		     The exact formulation of this problem is the following.
 		     Let \f$G=(V,A)\f$ be a digraph,
 		     \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$,
 		     \f$delta: V\rightarrow\mathbf{R}\f$. Find a feasible circulation
 		     \f$f: A\rightarrow\mathbf{R}^+_0\f$ so that
 		     \f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a)
 		     \geq delta(v) \quad \forall v\in V, \f]
 		     \f[ lower(a)\leq f(a) \leq upper(a) \quad \forall a\in A. \f]
 		     \note \f$delta(v)\f$ specifies a lower bound for the supply of node
 		     \f$v\f$. It can be either positive or negative, however note that
 		     \f$\sum_{v\in V}delta(v)\f$ should be zero or negative in order to
 		     have a feasible solution.
 		     \note A special case of this problem is when
 		     \f$\sum_{v\in V}delta(v) = 0\f$. Then the supply of each node \f$v\f$
 		     will be \e equal \e to \f$delta(v)\f$, if a circulation can be found.
 		     Thus a feasible solution for the
 		     \ref min_cost_flow "minimum cost flow" problem can be calculated
 		     in this way.
 		     \tparam GR The type of the digraph the algorithm runs on.
 		     \tparam LM The type of the lower bound capacity map. The default
 		     map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
 		     \tparam UM The type of the upper bound capacity map. The default
 		     map type is \c LM.
 		     \tparam DM The type of the map that stores the lower bound
 		     for the supply of the nodes. The default map type is
 		     \ref concepts::Digraph::NodeMap "GR::NodeMap<UM::Value>".
 		  */
 		#ifdef DOXYGEN
 		template< typename GR,
 		          typename LM,
 		          typename UM,
 		          typename DM,
 		          typename TR >
 		#else
 		template< typename GR,
 		          typename LM = typename GR::template ArcMap<int>,
 		          typename UM = LM,
 		          typename DM = typename GR::template NodeMap<typename UM::Value>,
 		          typename TR = CirculationDefaultTraits<GR, LM, UM, DM> >
 		#endif
 		  class Circulation {
 		  public:
 		    ///The \ref CirculationDefaultTraits "traits class" of the algorithm.
 		    typedef TR Traits;
 		    ///The type of the digraph the algorithm runs on.
 		    typedef typename Traits::Digraph Digraph;
 		    ///The type of the flow values.
 		    typedef typename Traits::Value Value;
 		    /// The type of the lower bound capacity map.
 		    typedef typename Traits::LCapMap LCapMap;
 		    /// The type of the upper bound capacity map.
 		    typedef typename Traits::UCapMap UCapMap;
 		    /// \brief The type of the map that stores the lower bound for
 		    /// the supply of the nodes.
 		    typedef typename Traits::DeltaMap DeltaMap;
 		    ///The type of the flow map.
 		    typedef typename Traits::FlowMap FlowMap;
 		    ///The type of the elevator.
 		    typedef typename Traits::Elevator Elevator;
 		    ///The type of the tolerance.
 		    typedef typename Traits::Tolerance Tolerance;
 		  private:
 		    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
 		    const Digraph &_g;
 		    int _node_num;
 		    const LCapMap *_lo;
 		    const UCapMap *_up;
 		    const DeltaMap *_delta;
 		    FlowMap *_flow;
 		    bool _local_flow;
 		    Elevator* _level;
 		    bool _local_level;
 		    typedef typename Digraph::template NodeMap<Value> ExcessMap;
 		    ExcessMap* _excess;
 		    Tolerance _tol;
 		    int _el;
 		  public:
 		    typedef Circulation Create;
 		    ///\name Named Template Parameters
 		    ///@{
 		    template <typename T>
 		    struct SetFlowMapTraits : public Traits {
 		      typedef T FlowMap;
 		      static FlowMap *createFlowMap(const Digraph&) {
 		        LEMON_ASSERT(false, "FlowMap is not initialized");
 		        return 0; // ignore warnings
+		      }
 		    };
 		    /// \brief \ref named-templ-param "Named parameter" for setting
 		    /// FlowMap type
 		    ///
 		    /// \ref named-templ-param "Named parameter" for setting FlowMap
 		    /// type.
 		    template <typename T>
 		    struct SetFlowMap
 		      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
 		                           SetFlowMapTraits<T> > {
 		      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
 		                          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 Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
 		                           SetElevatorTraits<T> > {
 		      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
 		                          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 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

@@ -934,913 +934,913 @@
 		    std::vector<_core_bits::MapCopyBase<From, Arc, ArcRefMap>* >
 		      _arc_maps;
 		    std::vector<_core_bits::MapCopyBase<From, Edge, EdgeRefMap>* >
 		      _edge_maps;
 		  };
 		  /// \brief Copy a graph to another graph.
 		  ///
 		  /// This function copies a graph to another graph.
 		  /// The complete usage of it is detailed in the GraphCopy class,
 		  /// but a short example shows a basic work:
 		  ///\code
 		  /// graphCopy(src, trg).nodeRef(nr).edgeCrossRef(ecr).run();
 		  ///\endcode
 		  ///
 		  /// After the copy the \c nr map will contain the mapping from the
 		  /// nodes of the \c from graph to the nodes of the \c to graph and
 		  /// \c ecr will contain the mapping from the edges of the \c to graph
 		  /// to the edges of the \c from graph.
 		  ///
 		  /// \see GraphCopy
 		  template <typename From, typename To>
 		  GraphCopy<From, To>
 		  graphCopy(const From& from, To& to) {
 		    return GraphCopy<From, To>(from, to);
+		  }
 		  namespace _core_bits {
 		    template <typename Graph, typename Enable = void>
 		    struct FindArcSelector {
 		      typedef typename Graph::Node Node;
 		      typedef typename Graph::Arc Arc;
 		      static Arc find(const Graph &g, Node u, Node v, Arc e) {
 		        if (e == INVALID) {
 		          g.firstOut(e, u);
 		        } else {
 		          g.nextOut(e);
+		        }
 		        while (e != INVALID && g.target(e) != v) {
 		          g.nextOut(e);
+		        }
 		        return e;
+		      }
 		    };
 		    template <typename Graph>
 		    struct FindArcSelector<
 		      Graph,
 		      typename enable_if<typename Graph::FindArcTag, void>::type>
+		    {
 		      typedef typename Graph::Node Node;
 		      typedef typename Graph::Arc Arc;
 		      static Arc find(const Graph &g, Node u, Node v, Arc prev) {
 		        return g.findArc(u, v, prev);
+		      }
 		    };
+		  }
 		  /// \brief Find an arc between two nodes of a digraph.
 		  ///
 		  /// This function finds an arc from node \c u to node \c v in the
 		  /// digraph \c g.
 		  ///
 		  /// If \c prev is \ref INVALID (this is the default value), then
 		  /// it finds the first arc from \c u to \c v. Otherwise it looks for
 		  /// the next arc from \c u to \c v after \c prev.
 		  /// \return The found arc or \ref INVALID if there is no such an arc.
 		  ///
 		  /// Thus you can iterate through each arc from \c u to \c v as it follows.
 		  ///\code
 		  /// for(Arc e = findArc(g,u,v); e != INVALID; e = findArc(g,u,v,e)) {
 		  ///   ...
 		  /// }
 		  ///\endcode
 		  ///
 		  /// \note \ref ConArcIt provides iterator interface for the same
 		  /// functionality.
 		  ///
 		  ///\sa ConArcIt
 		  ///\sa ArcLookUp, AllArcLookUp, DynArcLookUp
 		  template <typename Graph>
 		  inline typename Graph::Arc
 		  findArc(const Graph &g, typename Graph::Node u, typename Graph::Node v,
 		          typename Graph::Arc prev = INVALID) {
 		    return _core_bits::FindArcSelector<Graph>::find(g, u, v, prev);
+		  }
 		  /// \brief Iterator for iterating on parallel arcs connecting the same nodes.
 		  ///
 		  /// Iterator for iterating on parallel arcs connecting the same nodes. It is
 		  /// a higher level interface for the \ref findArc() function. You can
 		  /// use it the following way:
 		  ///\code
 		  /// for (ConArcIt<Graph> it(g, src, trg); it != INVALID; ++it) {
 		  ///   ...
 		  /// }
 		  ///\endcode
 		  ///
 		  ///\sa findArc()
 		  ///\sa ArcLookUp, AllArcLookUp, DynArcLookUp
 		  template <typename GR>
 		  class ConArcIt : public GR::Arc {
 		  public:
 		    typedef GR Graph;
 		    typedef typename Graph::Arc Parent;
 		    typedef typename Graph::Arc Arc;
 		    typedef typename Graph::Node Node;
 		    /// \brief Constructor.
 		    ///
 		    /// Construct a new ConArcIt iterating on the arcs that
 		    /// connects nodes \c u and \c v.
 		    ConArcIt(const Graph& g, Node u, Node v) : _graph(g) {
 		      Parent::operator=(findArc(_graph, u, v));
+		    }
 		    /// \brief Constructor.
 		    ///
 		    /// Construct a new ConArcIt that continues the iterating from arc \c a.
 		    ConArcIt(const Graph& g, Arc a) : Parent(a), _graph(g) {}
 		    /// \brief Increment operator.
 		    ///
 		    /// It increments the iterator and gives back the next arc.
 		    ConArcIt& operator++() {
 		      Parent::operator=(findArc(_graph, _graph.source(*this),
 		                                _graph.target(*this), *this));
 		      return *this;
+		    }
 		  private:
 		    const Graph& _graph;
 		  };
 		  namespace _core_bits {
 		    template <typename Graph, typename Enable = void>
 		    struct FindEdgeSelector {
 		      typedef typename Graph::Node Node;
 		      typedef typename Graph::Edge Edge;
 		      static Edge find(const Graph &g, Node u, Node v, Edge e) {
 		        bool b;
 		        if (u != v) {
 		          if (e == INVALID) {
 		            g.firstInc(e, b, u);
 		          } else {
 		            b = g.u(e) == u;
 		            g.nextInc(e, b);
+		          }
 		          while (e != INVALID && (b ? g.v(e) : g.u(e)) != v) {
 		            g.nextInc(e, b);
+		          }
 		        } else {
 		          if (e == INVALID) {
 		            g.firstInc(e, b, u);
 		          } else {
 		            b = true;
 		            g.nextInc(e, b);
+		          }
 		          while (e != INVALID && (!b || g.v(e) != v)) {
 		            g.nextInc(e, b);
+		          }
+		        }
 		        return e;
+		      }
 		    };
 		    template <typename Graph>
 		    struct FindEdgeSelector<
 		      Graph,
 		      typename enable_if<typename Graph::FindEdgeTag, void>::type>
+		    {
 		      typedef typename Graph::Node Node;
 		      typedef typename Graph::Edge Edge;
 		      static Edge find(const Graph &g, Node u, Node v, Edge prev) {
 		        return g.findEdge(u, v, prev);
+		      }
 		    };
+		  }
 		  /// \brief Find an edge between two nodes of a graph.
 		  ///
 		  /// This function finds an edge from node \c u to node \c v in graph \c g.
 		  /// If node \c u and node \c v is equal then each loop edge
 		  /// will be enumerated once.
 		  ///
 		  /// If \c prev is \ref INVALID (this is the default value), then
 		  /// it finds the first edge from \c u to \c v. Otherwise it looks for
 		  /// 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);
+		      }
 		      else _head[n]=INVALID;
+		    }
 		    ///Refresh the full data structure.
 		    ///Build up the full search database. In fact, it simply calls
 		    ///\ref refresh(Node) "refresh(n)" for each node \c n.
 		    ///
 		    ///It runs in time <em>O</em>(<em>m</em> log<em>D</em>), where <em>m</em> is
 		    ///the number of the arcs in the digraph and <em>D</em> is the maximum
 		    ///out-degree of the digraph.
 		    void refresh()
+		    {
 		      for(NodeIt n(_g);n!=INVALID;++n) refresh(n);
+		    }
 		    ///Find an arc between two nodes.
 		    ///Find an arc between two nodes in time <em>O</em>(log<em>d</em>),
 		    ///where <em>d</em> is the number of outgoing arcs of \c s.
 		    ///\param s The source node.
 		    ///\param t The target node.
 		    ///\return An arc from \c s to \c t if there exists,
 		    ///\ref INVALID otherwise.
 		    ///
 		    ///\warning If you change the digraph, refresh() must be called before using
 		    ///this operator. If you change the outgoing arcs of
 		    ///a single node \c n, then \ref refresh(Node) "refresh(n)" is enough.
 		    Arc operator()(Node s, Node t) const
+		    {
 		      Arc e;
 		      for(e=_head[s];
 		          e!=INVALID&&_g.target(e)!=t;
 		          e = t < _g.target(e)?_left[e]:_right[e]) ;
 		      return e;
+		    }
 		  };
 		  ///Fast look-up of all arcs between given endpoints.
 		  ///This class is the same as \ref ArcLookUp, with the addition
 		  ///that it makes it possible to find all parallel arcs between given
 		  ///endpoints.
 		  ///
 		  ///\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 ArcLookUp
 		  template<class GR>
 		  class AllArcLookUp : public ArcLookUp<GR>
+		  {
 		    using ArcLookUp<GR>::_g;
 		    using ArcLookUp<GR>::_right;
 		    using ArcLookUp<GR>::_left;
 		    using ArcLookUp<GR>::_head;
 		    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
 		    typedef GR Digraph;
 		    typename Digraph::template ArcMap<Arc> _next;
 		    Arc refreshNext(Arc head,Arc next=INVALID)
+		    {
 		      if(head==INVALID) return next;
 		      else {
 		        next=refreshNext(_right[head],next);
 		        _next[head]=( next!=INVALID && _g.target(next)==_g.target(head))
 		          ? next : INVALID;
 		        return refreshNext(_left[head],head);
+		      }
+		    }
 		    void refreshNext()
+		    {
 		      for(NodeIt n(_g);n!=INVALID;++n) refreshNext(_head[n]);
+		    }
 		  public:
 		    ///Constructor
 		    ///Constructor.
 		    ///
 		    ///It builds up the search database, which remains valid until the digraph
 		    ///changes.
 		    AllArcLookUp(const Digraph &g) : ArcLookUp<GR>(g), _next(g) {refreshNext();}
 		    ///Refresh the 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)
+		    {
 		      ArcLookUp<GR>::refresh(n);
 		      refreshNext(_head[n]);
+		    }
 		    ///Refresh the full data structure.
 		    ///Build up the full search database. In fact, it simply calls
 		    ///\ref refresh(Node) "refresh(n)" for each node \c n.
 		    ///
 		    ///It runs in time <em>O</em>(<em>m</em> log<em>D</em>), where <em>m</em> is
 		    ///the number of the arcs in the digraph and <em>D</em> is the maximum
 		    ///out-degree of the digraph.
 		    void refresh()
+		    {
 		      for(NodeIt n(_g);n!=INVALID;++n) refresh(_head[n]);
+		    }
 		    ///Find an arc between two nodes.
 		    ///Find an arc between two nodes.
 		    ///\param s The source node.
 		    ///\param t The target node.
 		    ///\param prev 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 prev 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
 		    ///AllArcLookUp<ListDigraph> ae(g);
 		    ///...
 		    ///int n = 0;
 		    ///for(Arc a = ae(u,v); a != INVALID; a=ae(u,v,a)) n++;
 		    ///\endcode
 		    ///
 		    ///Finding the first arc take <em>O</em>(log<em>d</em>) time,
 		    ///where <em>d</em> is the number of outgoing arcs of \c s. Then the
 		    ///consecutive arcs are found in constant time.
 		    ///
 		    ///\warning If you change the digraph, refresh() must be called before using
 		    ///this operator. If you change the outgoing arcs of
 		    ///a single node \c n, then \ref refresh(Node) "refresh(n)" is enough.
 		    ///
 		#ifdef DOXYGEN
 		    Arc operator()(Node s, Node t, Arc prev=INVALID) const {}
 		#else
 		    using ArcLookUp<GR>::operator() ;
 		    Arc operator()(Node s, Node t, Arc prev) const
+		    {
 		      return prev==INVALID?(*this)(s,t):_next[prev];
+		    }
 		#endif
 		  };
 		  /// @}
 		} //namespace lemon
 		#endif

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
 		      /// Conversion to \c Arc.
 		      ///
 		      operator typename Graph::Arc() const
+		      {
 		        return _arc_it;
+		      }
 		      /// Conversion to \c Edge
 		      /// Conversion to \c Edge.
 		      ///
 		      operator typename Graph::Edge() const
+		      {
 		        return _arc_it;
+		      }
 		      bool operator==(Invalid) { return _node_it==INVALID; }
 		      bool operator!=(Invalid) { return _node_it!=INVALID; }
 		      /// Next edge
 		      /// Next edge.
 		      ///
 		      MinCutEdgeIt &operator++()
+		      {
 		        step();
 		        while(_arc_it!=INVALID && _cut[_graph.target(_arc_it)]) step();
 		        return *this;
+		      }
 		      /// Postfix incrementation
 		      /// Postfix incrementation.
 		      ///
 		      /// \warning This incrementation
 		      /// returns an \c Arc, not a \c MinCutEdgeIt, as one may expect.
 		      typename Graph::Arc operator++(int)
+		      {
 		        typename Graph::Arc e=*this;
 		        ++(*this);
 		        return e;
+		      }
 		    };
 		  };
+		}
 		#endif

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

 		/* -*- 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_MAX_MATCHING_H
 		#define LEMON_MAX_MATCHING_H
 		#include <vector>
 		#include <queue>
 		#include <set>
 		#include <limits>
 		#include <lemon/core.h>
 		#include <lemon/unionfind.h>
 		#include <lemon/bin_heap.h>
 		#include <lemon/maps.h>
 		///\ingroup matching
 		///\file
 		///\brief Maximum matching algorithms in general graphs.
 		namespace lemon {
 		  /// \ingroup matching
 		  ///
 		  /// \brief Edmonds' alternating forest maximum matching algorithm.
 		  ///
 		  /// This class implements Edmonds' alternating forest matching
 		  /// algorithm. The algorithm can be started from an arbitrary initial
 		  /// matching (the default is the empty one)
 		  ///
 		  /// The dual solution of the problem is a map of the nodes to
 		  /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c
 		  /// MATCHED/C showing the Gallai-Edmonds decomposition of the
 		  /// graph. The nodes in \c EVEN/D induce a graph with
 		  /// factor-critical components, the nodes in \c ODD/A form the
 		  /// barrier, and the nodes in \c MATCHED/C induce a graph having a
 		  /// perfect matching. The number of the factor-critical components
 		  /// minus the number of barrier nodes is a lower bound on the
 		  /// unmatched nodes, and the matching is optimal if and only if this bound is
 		  /// tight. This decomposition can be attained by calling \c
 		  /// decomposition() after running the algorithm.
 		  ///
 		  /// \param GR The graph type the algorithm runs on.
 		  template <typename GR>
 		  class MaxMatching {
 		  public:
 		    typedef GR Graph;
 		    typedef typename Graph::template NodeMap<typename Graph::Arc>
 		    MatchingMap;
 		    ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
 		    ///
 		    ///Indicates the Gallai-Edmonds decomposition of the graph. The
 		    ///nodes with Status \c EVEN/D induce a graph with factor-critical
 		    ///components, the nodes in \c ODD/A form the canonical barrier,
 		    ///and the nodes in \c MATCHED/C induce a graph having a perfect
 		    ///matching.
 		    enum Status {
 		      EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2
 		    };
 		    typedef typename Graph::template NodeMap<Status> StatusMap;
 		  private:
 		    TEMPLATE_GRAPH_TYPEDEFS(Graph);
 		    typedef UnionFindEnum<IntNodeMap> BlossomSet;
 		    typedef ExtendFindEnum<IntNodeMap> TreeSet;
 		    typedef RangeMap<Node> NodeIntMap;
 		    typedef MatchingMap EarMap;
 		    typedef std::vector<Node> NodeQueue;
 		    const Graph& _graph;
 		    MatchingMap* _matching;
 		    StatusMap* _status;
 		    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
 		    };
 		    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
 		    struct BlossomData {
 		      int tree;
 		      Status status;
 		      Arc pred, next;
 		      Value pot, offset;
 		      Node base;
 		    };
 		    IntNodeMap *_blossom_index;
 		    BlossomSet *_blossom_set;
 		    RangeMap<BlossomData>* _blossom_data;
 		    IntNodeMap *_node_index;
 		    IntArcMap *_node_heap_index;
 		    struct NodeData {
 		      NodeData(IntArcMap& node_heap_index)
 		        : heap(node_heap_index) {}
 		      int blossom;
 		      Value pot;
 		      BinHeap<Value, IntArcMap> heap;
 		      std::map<int, Arc> heap_index;
 		      int tree;
 		    };
 		    RangeMap<NodeData>* _node_data;
 		    typedef ExtendFindEnum<IntIntMap> TreeSet;
 		    IntIntMap *_tree_set_index;
 		    TreeSet *_tree_set;
 		    IntNodeMap *_delta1_index;
 		    BinHeap<Value, IntNodeMap> *_delta1;
 		    IntIntMap *_delta2_index;
 		    BinHeap<Value, IntIntMap> *_delta2;
 		    IntEdgeMap *_delta3_index;
 		    BinHeap<Value, IntEdgeMap> *_delta3;
 		    IntIntMap *_delta4_index;
 		    BinHeap<Value, IntIntMap> *_delta4;
 		    Value _delta_sum;
 		    void createStructures() {
 		      _node_num = countNodes(_graph);
 		      _blossom_num = _node_num * 3 / 2;
 		      if (!_matching) {
 		        _matching = new MatchingMap(_graph);
+		      }
 		      if (!_node_potential) {
 		        _node_potential = new NodePotential(_graph);
+		      }
 		      if (!_blossom_set) {
 		        _blossom_index = new IntNodeMap(_graph);
 		        _blossom_set = new BlossomSet(*_blossom_index);
 		        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
+		      }
 		      if (!_node_index) {
 		        _node_index = new IntNodeMap(_graph);
 		        _node_heap_index = new IntArcMap(_graph);
 		        _node_data = new RangeMap<NodeData>(_node_num,
 		                                              NodeData(*_node_heap_index));
+		      }
 		      if (!_tree_set) {
 		        _tree_set_index = new IntIntMap(_blossom_num);
 		        _tree_set = new TreeSet(*_tree_set_index);
+		      }
 		      if (!_delta1) {
 		        _delta1_index = new IntNodeMap(_graph);
 		        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
+		      }
 		      if (!_delta2) {
 		        _delta2_index = new IntIntMap(_blossom_num);
 		        _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
+		      }
 		      if (!_delta3) {
 		        _delta3_index = new IntEdgeMap(_graph);
 		        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
+		      }
 		      if (!_delta4) {
 		        _delta4_index = new IntIntMap(_blossom_num);
 		        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
+		      }
+		    }
 		    void destroyStructures() {
 		      _node_num = countNodes(_graph);
 		      _blossom_num = _node_num * 3 / 2;
 		      if (_matching) {
 		        delete _matching;
+		      }
 		      if (_node_potential) {
 		        delete _node_potential;
+		      }
 		      if (_blossom_set) {
 		        delete _blossom_index;
 		        delete _blossom_set;
 		        delete _blossom_data;
+		      }
 		      if (_node_index) {
 		        delete _node_index;
 		        delete _node_heap_index;
 		        delete _node_data;
+		      }
 		      if (_tree_set) {
 		        delete _tree_set_index;
 		        delete _tree_set;
+		      }
 		      if (_delta1) {
 		        delete _delta1_index;
 		        delete _delta1;
+		      }
 		      if (_delta2) {
 		        delete _delta2_index;
 		        delete _delta2;
+		      }
 		      if (_delta3) {
 		        delete _delta3_index;
 		        delete _delta3;
+		      }
 		      if (_delta4) {
 		        delete _delta4_index;
 		        delete _delta4;
+		      }
+		    }
 		    void matchedToEven(int blossom, int tree) {
 		      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
 		        _delta2->erase(blossom);
+		      }
 		      if (!_blossom_set->trivial(blossom)) {
 		        (*_blossom_data)[blossom].pot -=
 * (_delta_sum - (*_blossom_data)[blossom].offset);
+		      }
 		      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
 		           n != INVALID; ++n) {
 		        _blossom_set->increase(n, std::numeric_limits<Value>::max());
 		        int ni = (*_node_index)[n];
 		        (*_node_data)[ni].heap.clear();
 		        (*_node_data)[ni].heap_index.clear();
 		        (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
 		        _delta1->push(n, (*_node_data)[ni].pot);
 		        for (InArcIt e(_graph, n); e != INVALID; ++e) {
 		          Node v = _graph.source(e);
 		          int vb = _blossom_set->find(v);
 		          int vi = (*_node_index)[v];
 		          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
 		            dualScale * _weight[e];
 		          if ((*_blossom_data)[vb].status == EVEN) {
 		            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
 		              _delta3->push(e, rw / 2);
+		            }
 		          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
 		            if (_delta3->state(e) != _delta3->IN_HEAP) {
 		              _delta3->push(e, rw);
+		            }
 		          } else {
 		            typename std::map<int, Arc>::iterator it =
 		              (*_node_data)[vi].heap_index.find(tree);
 		            if (it != (*_node_data)[vi].heap_index.end()) {
 		              if ((*_node_data)[vi].heap[it->second] > rw) {
 		                (*_node_data)[vi].heap.replace(it->second, e);
 		                (*_node_data)[vi].heap.decrease(e, rw);
 		                it->second = e;
+		              }
 		            } else {
 		              (*_node_data)[vi].heap.push(e, rw);
 		              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
@@ -1167,887 +1167,887 @@
 		          } else if ((*_blossom_data)[vb].status == EVEN) {
 		            if (_delta3->state(e) == _delta3->IN_HEAP) {
 		              _delta3->erase(e);
+		            }
 		            int vt = _tree_set->find(vb);
 		            Arc r = _graph.oppositeArc(e);
 		            typename std::map<int, Arc>::iterator it =
 		              (*_node_data)[ni].heap_index.find(vt);
 		            if (it != (*_node_data)[ni].heap_index.end()) {
 		              if ((*_node_data)[ni].heap[it->second] > rw) {
 		                (*_node_data)[ni].heap.replace(it->second, r);
 		                (*_node_data)[ni].heap.decrease(r, rw);
 		                it->second = r;
+		              }
 		            } else {
 		              (*_node_data)[ni].heap.push(r, rw);
 		              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
+		            }
 		            if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
 		              _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
 		              if (_delta2->state(blossom) != _delta2->IN_HEAP) {
 		                _delta2->push(blossom, _blossom_set->classPrio(blossom) -
 		                             (*_blossom_data)[blossom].offset);
 		              } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
 		                         (*_blossom_data)[blossom].offset){
 		                _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
 		                                 (*_blossom_data)[blossom].offset);
+		              }
+		            }
 		          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
 		            if (_delta3->state(e) == _delta3->IN_HEAP) {
 		              _delta3->erase(e);
+		            }
+		          }
+		        }
+		      }
+		    }
 		    void alternatePath(int even, int tree) {
 		      int odd;
 		      evenToMatched(even, tree);
 		      (*_blossom_data)[even].status = MATCHED;
 		      while ((*_blossom_data)[even].pred != INVALID) {
 		        odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
 		        (*_blossom_data)[odd].status = MATCHED;
 		        oddToMatched(odd);
 		        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
 		        even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
 		        (*_blossom_data)[even].status = MATCHED;
 		        evenToMatched(even, tree);
 		        (*_blossom_data)[even].next =
 		          _graph.oppositeArc((*_blossom_data)[odd].pred);
+		      }
+		    }
 		    void destroyTree(int tree) {
 		      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
 		        if ((*_blossom_data)[b].status == EVEN) {
 		          (*_blossom_data)[b].status = MATCHED;
 		          evenToMatched(b, tree);
 		        } else if ((*_blossom_data)[b].status == ODD) {
 		          (*_blossom_data)[b].status = MATCHED;
 		          oddToMatched(b);
+		        }
+		      }
 		      _tree_set->eraseClass(tree);
+		    }
 		    void unmatchNode(const Node& node) {
 		      int blossom = _blossom_set->find(node);
 		      int tree = _tree_set->find(blossom);
 		      alternatePath(blossom, tree);
 		      destroyTree(tree);
 		      (*_blossom_data)[blossom].status = UNMATCHED;
 		      (*_blossom_data)[blossom].base = node;
 		      matchedToUnmatched(blossom);
+		    }
 		    void augmentOnEdge(const Edge& edge) {
 		      int left = _blossom_set->find(_graph.u(edge));
 		      int right = _blossom_set->find(_graph.v(edge));
 		      if ((*_blossom_data)[left].status == EVEN) {
 		        int left_tree = _tree_set->find(left);
 		        alternatePath(left, left_tree);
 		        destroyTree(left_tree);
 		      } else {
 		        (*_blossom_data)[left].status = MATCHED;
 		        unmatchedToMatched(left);
+		      }
 		      if ((*_blossom_data)[right].status == EVEN) {
 		        int right_tree = _tree_set->find(right);
 		        alternatePath(right, right_tree);
 		        destroyTree(right_tree);
 		      } else {
 		        (*_blossom_data)[right].status = MATCHED;
 		        unmatchedToMatched(right);
+		      }
 		      (*_blossom_data)[left].next = _graph.direct(edge, true);
 		      (*_blossom_data)[right].next = _graph.direct(edge, false);
+		    }
 		    void extendOnArc(const Arc& arc) {
 		      int base = _blossom_set->find(_graph.target(arc));
 		      int tree = _tree_set->find(base);
 		      int odd = _blossom_set->find(_graph.source(arc));
 		      _tree_set->insert(odd, tree);
 		      (*_blossom_data)[odd].status = ODD;
 		      matchedToOdd(odd);
 		      (*_blossom_data)[odd].pred = arc;
 		      int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
 		      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
 		      _tree_set->insert(even, tree);
 		      (*_blossom_data)[even].status = EVEN;
 		      matchedToEven(even, tree);
+		    }
 		    void shrinkOnEdge(const Edge& edge, int tree) {
 		      int nca = -1;
 		      std::vector<int> left_path, right_path;
+		      {
 		        std::set<int> left_set, right_set;
 		        int left = _blossom_set->find(_graph.u(edge));
 		        left_path.push_back(left);
 		        left_set.insert(left);
 		        int right = _blossom_set->find(_graph.v(edge));
 		        right_path.push_back(right);
 		        right_set.insert(right);
 		        while (true) {
 		          if ((*_blossom_data)[left].pred == INVALID) break;
 		          left =
 		            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
 		          left_path.push_back(left);
 		          left =
 		            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
 		          left_path.push_back(left);
 		          left_set.insert(left);
 		          if (right_set.find(left) != right_set.end()) {
 		            nca = left;
 		            break;
+		          }
 		          if ((*_blossom_data)[right].pred == INVALID) break;
 		          right =
 		            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
 		          right_path.push_back(right);
 		          right =
 		            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
 		          right_path.push_back(right);
 		          right_set.insert(right);
 		          if (left_set.find(right) != left_set.end()) {
 		            nca = right;
 		            break;
+		          }
+		        }
 		        if (nca == -1) {
 		          if ((*_blossom_data)[left].pred == INVALID) {
 		            nca = right;
 		            while (left_set.find(nca) == left_set.end()) {
 		              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 {
 		      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.
 		    /// \see BlossomIt
 		    Value blossomValue(int k) const {
 		      return _blossom_potential[k].value;
+		    }
 		    /// \brief Iterator for obtaining the nodes of the blossom.
 		    ///
 		    /// Iterator for obtaining the nodes of the blossom. This class
 		    /// provides a common lemon style iterator for listing a
 		    /// subset of the nodes.
 		    class BlossomIt {
 		    public:
 		      /// \brief Constructor.
 		      ///
 		      /// Constructor to get the nodes of the variable.
 		      BlossomIt(const MaxWeightedMatching& algorithm, int variable)
 		        : _algorithm(&algorithm)
+		      {
 		        _index = _algorithm->_blossom_potential[variable].begin;
 		        _last = _algorithm->_blossom_potential[variable].end;
+		      }
 		      /// \brief Conversion to node.
 		      ///
 		      /// Conversion to node.
 		      operator Node() const {
 		        return _algorithm->_blossom_node_list[_index];
+		      }
 		      /// \brief Increment operator.
 		      ///
 		      /// Increment operator.
 		      BlossomIt& 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 MaxWeightedMatching* _algorithm;
 		      int _last;
 		      int _index;
 		    };
 		    /// @}
 		  };
 		  /// \ingroup matching
 		  ///
 		  /// \brief Weighted perfect matching in general graphs
 		  ///
 		  /// This class provides an efficient implementation of Edmond's
 		  /// maximum weighted perfect 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 and covers all nodes. The problem can be
 		  /// formulated with the following linear program.
 		  /// \f[ \sum_{e \in \delta(u)}x_e = 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[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 MaxWeightedPerfectMatching {
 		  public:
 		    typedef GR Graph;
 		    typedef WM WeightMap;
 		    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
 		    };
 		    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
 		    struct BlossomData {
@@ -2360,748 +2360,748 @@
 		        _delta4->erase(blossom);
+		      }
+		    }
 		    void oddToEven(int blossom, int tree) {
 		      if (!_blossom_set->trivial(blossom)) {
 		        _delta4->erase(blossom);
 		        (*_blossom_data)[blossom].pot -=
 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
+		      }
 		      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
 		           n != INVALID; ++n) {
 		        int ni = (*_node_index)[n];
 		        _blossom_set->increase(n, std::numeric_limits<Value>::max());
 		        (*_node_data)[ni].heap.clear();
 		        (*_node_data)[ni].heap_index.clear();
 		        (*_node_data)[ni].pot +=
 * _delta_sum - (*_blossom_data)[blossom].offset;
 		        for (InArcIt e(_graph, n); e != INVALID; ++e) {
 		          Node v = _graph.source(e);
 		          int vb = _blossom_set->find(v);
 		          int vi = (*_node_index)[v];
 		          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
 		            dualScale * _weight[e];
 		          if ((*_blossom_data)[vb].status == EVEN) {
 		            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
 		              _delta3->push(e, rw / 2);
+		            }
 		          } else {
 		            typename std::map<int, Arc>::iterator it =
 		              (*_node_data)[vi].heap_index.find(tree);
 		            if (it != (*_node_data)[vi].heap_index.end()) {
 		              if ((*_node_data)[vi].heap[it->second] > rw) {
 		                (*_node_data)[vi].heap.replace(it->second, e);
 		                (*_node_data)[vi].heap.decrease(e, rw);
 		                it->second = e;
+		              }
 		            } else {
 		              (*_node_data)[vi].heap.push(e, rw);
 		              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
+		            }
 		            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
 		              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
 		              if ((*_blossom_data)[vb].status == MATCHED) {
 		                if (_delta2->state(vb) != _delta2->IN_HEAP) {
 		                  _delta2->push(vb, _blossom_set->classPrio(vb) -
 		                               (*_blossom_data)[vb].offset);
 		                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
 		                           (*_blossom_data)[vb].offset) {
 		                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
 		                                   (*_blossom_data)[vb].offset);
+		                }
+		              }
+		            }
+		          }
+		        }
+		      }
 		      (*_blossom_data)[blossom].offset = 0;
+		    }
 		    void alternatePath(int even, int tree) {
 		      int odd;
 		      evenToMatched(even, tree);
 		      (*_blossom_data)[even].status = MATCHED;
 		      while ((*_blossom_data)[even].pred != INVALID) {
 		        odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
 		        (*_blossom_data)[odd].status = MATCHED;
 		        oddToMatched(odd);
 		        (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
 		        even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
 		        (*_blossom_data)[even].status = MATCHED;
 		        evenToMatched(even, tree);
 		        (*_blossom_data)[even].next =
 		          _graph.oppositeArc((*_blossom_data)[odd].pred);
+		      }
+		    }
 		    void destroyTree(int tree) {
 		      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
 		        if ((*_blossom_data)[b].status == EVEN) {
 		          (*_blossom_data)[b].status = MATCHED;
 		          evenToMatched(b, tree);
 		        } else if ((*_blossom_data)[b].status == ODD) {
 		          (*_blossom_data)[b].status = MATCHED;
 		          oddToMatched(b);
+		        }
+		      }
 		      _tree_set->eraseClass(tree);
+		    }
 		    void augmentOnEdge(const Edge& edge) {
 		      int left = _blossom_set->find(_graph.u(edge));
 		      int right = _blossom_set->find(_graph.v(edge));
 		      int left_tree = _tree_set->find(left);
 		      alternatePath(left, left_tree);
 		      destroyTree(left_tree);
 		      int right_tree = _tree_set->find(right);
 		      alternatePath(right, right_tree);
 		      destroyTree(right_tree);
 		      (*_blossom_data)[left].next = _graph.direct(edge, true);
 		      (*_blossom_data)[right].next = _graph.direct(edge, false);
+		    }
 		    void extendOnArc(const Arc& arc) {
 		      int base = _blossom_set->find(_graph.target(arc));
 		      int tree = _tree_set->find(base);
 		      int odd = _blossom_set->find(_graph.source(arc));
 		      _tree_set->insert(odd, tree);
 		      (*_blossom_data)[odd].status = ODD;
 		      matchedToOdd(odd);
 		      (*_blossom_data)[odd].pred = arc;
 		      int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
 		      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
 		      _tree_set->insert(even, tree);
 		      (*_blossom_data)[even].status = EVEN;
 		      matchedToEven(even, tree);
+		    }
 		    void shrinkOnEdge(const Edge& edge, int tree) {
 		      int nca = -1;
 		      std::vector<int> left_path, right_path;
+		      {
 		        std::set<int> left_set, right_set;
 		        int left = _blossom_set->find(_graph.u(edge));
 		        left_path.push_back(left);
 		        left_set.insert(left);
 		        int right = _blossom_set->find(_graph.v(edge));
 		        right_path.push_back(right);
 		        right_set.insert(right);
 		        while (true) {
 		          if ((*_blossom_data)[left].pred == INVALID) break;
 		          left =
 		            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
 		          left_path.push_back(left);
 		          left =
 		            _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
 		          left_path.push_back(left);
 		          left_set.insert(left);
 		          if (right_set.find(left) != right_set.end()) {
 		            nca = left;
 		            break;
+		          }
 		          if ((*_blossom_data)[right].pred == INVALID) break;
 		          right =
 		            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
 		          right_path.push_back(right);
 		          right =
 		            _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
 		          right_path.push_back(right);
 		          right_set.insert(right);
 		          if (left_set.find(right) != left_set.end()) {
 		            nca = right;
 		            break;
+		          }
+		        }
 		        if (nca == -1) {
 		          if ((*_blossom_data)[left].pred == INVALID) {
 		            nca = right;
 		            while (left_set.find(nca) == left_set.end()) {
 		              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.
 		    /// \see BlossomIt
 		    Value blossomValue(int k) const {
 		      return _blossom_potential[k].value;
+		    }
 		    /// \brief Iterator for obtaining the nodes of the blossom.
 		    ///
 		    /// Iterator for obtaining the nodes of the blossom. This class
 		    /// provides a common lemon style iterator for listing a
 		    /// subset of the nodes.
 		    class BlossomIt {
 		    public:
 		      /// \brief Constructor.
 		      ///
 		      /// Constructor to get the nodes of the variable.
 		      BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
 		        : _algorithm(&algorithm)
+		      {
 		        _index = _algorithm->_blossom_potential[variable].begin;
 		        _last = _algorithm->_blossom_potential[variable].end;
+		      }
 		      /// \brief Conversion to node.
 		      ///
 		      /// Conversion to node.
 		      operator Node() const {
 		        return _algorithm->_blossom_node_list[_index];
+		      }
 		      /// \brief Increment operator.
 		      ///
 		      /// Increment operator.
 		      BlossomIt& 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 MaxWeightedPerfectMatching* _algorithm;
 		      int _last;
 		      int _index;
 		    };
 		    /// @}
 		  };
 		} //END OF NAMESPACE LEMON
 		#endif //LEMON_MAX_MATCHING_H

lemon/min_cost_arborescence.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-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_MIN_COST_ARBORESCENCE_H
 		#define LEMON_MIN_COST_ARBORESCENCE_H
 		///\ingroup spantree
 		///\file
 		///\brief Minimum Cost Arborescence algorithm.
 		#include <vector>
 		#include <lemon/list_graph.h>
 		#include <lemon/bin_heap.h>
 		#include <lemon/assert.h>
 		namespace lemon {
 		  /// \brief Default traits class for MinCostArborescence class.
 		  ///
 		  /// Default traits class for MinCostArborescence class.
 		  /// \param GR Digraph type.
 		  /// \param CM Type of cost map.
 		  template <class GR, class CM>
 		  struct MinCostArborescenceDefaultTraits{
 		    /// \brief The digraph type the algorithm runs on.
 		    typedef GR Digraph;
 		    /// \brief The type of the map that stores the arc costs.
 		    ///
 		    /// The type of the map that stores the arc costs.
 		    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
 		    typedef CM CostMap;
 		    /// \brief The value type of the costs.
 		    ///
 		    /// The value type of the costs.
 		    typedef typename CostMap::Value Value;
 		    /// \brief The type of the map that stores which arcs are in the
 		    /// arborescence.
 		    ///
 		    /// The type of the map that stores which arcs are in the
 		    /// arborescence.  It must meet the \ref concepts::WriteMap
 		    /// "WriteMap" concept.  Initially it will be set to false on each
 		    /// arc. After it will set all arborescence arcs once.
 		    typedef typename Digraph::template ArcMap<bool> ArborescenceMap;
 		    /// \brief Instantiates a \c ArborescenceMap.
 		    ///
 		    /// This function instantiates a \c ArborescenceMap.
 		    /// \param digraph is the graph, to which we would like to
 		    /// calculate the \c ArborescenceMap.
 		    static ArborescenceMap *createArborescenceMap(const Digraph &digraph){
 		      return new ArborescenceMap(digraph);
+		    }
 		    /// \brief The type of the \c PredMap
 		    ///
 		    /// The type of the \c PredMap. It is a node map with an arc value type.
 		    typedef typename Digraph::template NodeMap<typename Digraph::Arc> PredMap;
 		    /// \brief Instantiates a \c PredMap.
 		    ///
 		    /// This function instantiates a \c PredMap.
 		    /// \param digraph The digraph to which we would like to define the
 		    /// \c PredMap.
 		    static PredMap *createPredMap(const Digraph &digraph){
 		      return new PredMap(digraph);
+		    }
 		  };
 		  /// \ingroup spantree
 		  ///
 		  /// \brief %MinCostArborescence algorithm class.
 		  ///
 		  /// This class provides an efficient implementation of
 		  /// %MinCostArborescence algorithm. The arborescence is a tree
 		  /// which is directed from a given source node of the digraph. One or
 		  /// more sources should be given for the algorithm and it will calculate
 		  /// the minimum cost subgraph which are union of arborescences with the
 		  /// given sources and spans all the nodes which are reachable from the
 		  /// sources. The time complexity of the algorithm is O(n<sup>2</sup>+e).
 		  ///
 		  /// The algorithm provides also an optimal dual solution, therefore
 		  /// 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

Changeset was too big and was cut off... Show full diff

0 comments (0 inline)

RhodeCode

Login to your account

Changeset was too big and was cut off... Show full diff

Changeset was too big and was cut off... Show full diff