RhodeCode

#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.

  ///

  ///

  /// 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

  /// minus the number of barrier nodes is a lower bound on the

  /// tight. This decomposition can be attained by calling \c

  /// decomposition() after running the algorithm.

  ///

  /// \param _Graph The graph type the algorithm runs on.

  template <typename _Graph>

  class MaxMatching {

  public:

    typedef _Graph Graph;

    typedef typename Graph::template NodeMap<typename Graph::Arc>

    MatchingMap;

    ///\brief Indicates the Gallai-Edmonds decomposition of the graph.

    ///

    ///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:

    ///

    /// 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

    /// \c run() member function.

    /// \n

    /// \ref init(), \ref greedyInit() or \ref matchingInit()

    /// startParse() 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);

      }

    }

    ///

    void greedyInit() {

      createStructures();

      for (NodeIt n(_graph); n != INVALID; ++n) {

        _matching->set(n, INVALID);

        _status->set(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);

              break;

            }

          }

        }

      }

    }

    ///

    /// 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 %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);

        if ((*_status)[n] == UNMATCHED) {

          (*_blossom_rep)[_blossom_set->insert(n)] = n;

          _tree_set->insert(n);

          _status->set(n, EVEN);

          processSparse(n);

        }

      }

    }

    /// \brief Starts Edmonds' algorithm.

    ///

    /// It runs Edmonds' algorithm with a heuristic of postponing

    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);

          processDense(n);

        }

      }

    }

      if (countEdges(_graph) < 2 * countNodes(_graph)) {

        greedyInit();

        startSparse();

      } else {

        init();

        startDense();

      }

    }

    /// @}

    /// \name Primal solution

    /// @{

    ///

    ///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.

    ///\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

    /// @{

    /// \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];

    }

  /** \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

  /// of a blossom. If the value type is integral then the dual

  /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".

  template <typename _Graph,

            typename _WeightMap = typename _Graph::template EdgeMap<int> >

  class MaxWeightedMatching {

  public:

    typedef _Graph Graph;

    typedef _WeightMap WeightMap;

    typedef typename WeightMap::Value Value;

    /// \brief Scaling factor for dual solution

        _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

    /// \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);

    /// \code

    ///   mwm.init();

    ///   mwm.start();

    /// \endcode

    void run() {

      init();

      start();

    }

    /// @}

    /// \name Primal solution

    /// @{

    ///

    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.

    ///

    /// \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

    /// @{

    /// \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);

      }

    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;

    }

    ///

    /// subset of the nodes.

    class BlossomIt {

    public:

      /// \brief Constructor.

      ///

      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];

        _delta2_index(0), _delta2(0),

        _delta3_index(0), _delta3(0),

        _delta4_index(0), _delta4(0),

        _delta_sum() {}

    ~MaxWeightedPerfectMatching() {

      destroyStructures();

    }

    /// \name Execution control

    /// \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);

    /// \code

    ///   mwm.init();

    ///   mwm.start();

    /// \endcode

    bool run() {

      init();

      return start();

    }

    /// @}

    /// \name Primal solution

    /// @{

    /// \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]];

        }

    }

    /// \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

    /// @{

    /// \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);

      }

    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;

    }

    ///

    /// subset of the nodes.

    class BlossomIt {

    public:

      /// \brief Constructor.

      ///

      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];