RhodeCode

This group contains the algorithms for discovering the graph properties

like connectivity, bipartiteness, euler property, simplicity etc.

\image html edge_biconnected_components.png

\image latex edge_biconnected_components.eps "bi-edge-connected components" width=\textwidth

*/

/**

@defgroup planar Planarity Embedding and Drawing

@ingroup algs

\brief Algorithms for planarity checking, embedding and drawing

This group contains the algorithms for planarity checking,

embedding and drawing.

\image html planar.png

\image latex planar.eps "Plane graph" width=\textwidth

*/

/**

@defgroup matching Matching Algorithms

@ingroup algs

\brief Algorithms for finding matchings in graphs and bipartite graphs.

matchings in graphs and bipartite graphs. The general matching problem is

There are several different algorithms for calculate matchings in

graphs.  The matching problems in bipartite graphs are generally

easier than in general graphs. The goal of the matching optimization

can be finding maximum cardinality, maximum weight or minimum cost

matching. The search can be constrained to find perfect or

maximum cardinality matching.

The matching algorithms implemented in LEMON:

- \ref MaxBipartiteMatching Hopcroft-Karp augmenting path algorithm

  for calculating maximum cardinality matching in bipartite graphs.

- \ref PrBipartiteMatching Push-relabel algorithm

  for calculating maximum cardinality matching in bipartite graphs.

- \ref MaxWeightedBipartiteMatching

  Successive shortest path algorithm for calculating maximum weighted

  matching and maximum weighted bipartite matching in bipartite graphs.

- \ref MinCostMaxBipartiteMatching

  Successive shortest path algorithm for calculating minimum cost maximum

  matching in bipartite graphs.

- \ref MaxMatching Edmond's blossom shrinking algorithm for calculating

  maximum cardinality matching in general graphs.

- \ref MaxWeightedMatching Edmond's blossom shrinking algorithm for calculating

  maximum weighted matching in general graphs.

- \ref MaxWeightedPerfectMatching

 *

 */

#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

  ///

  ///

  ///

  /// The dual solution of the problem is a map of the nodes to

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

  /// decomposition() after running the algorithm.

  ///

  template <typename GR>

  class MaxMatching {

  public:

    typedef GR Graph;

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

    MatchingMap;

    ///

    enum Status {

    };

    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;

          Node n = node;

          while (n != base) {

            n = _graph.target((*_matching)[n]);

            Arc a = (*_ear)[n];

            n = _graph.target(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)[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)[_blossom_set->find(nca)] = nca;

    }

    void extendOnArc(const Arc& a) {

      Node base = _graph.source(a);

      Node odd = _graph.target(a);

      (*_ear)[odd] = _graph.oppositeArc(a);

      Node even = _graph.target((*_matching)[odd]);

      (*_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)[odd] = _graph.oppositeArc(a);

      (*_status)[odd] = MATCHED;

               jt != INVALID; ++jt) {

            (*_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();

    }

    /// The simplest way to execute the algorithm is to use the

    ///@{

    ///

    void init() {

      createStructures();

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

        (*_matching)[n] = INVALID;

        (*_status)[n] = UNMATCHED;

      }

    }

    ///

    void greedyInit() {

      createStructures();

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

        (*_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)[n] = a;

              (*_status)[n] = MATCHED;

              (*_matching)[v] = _graph.oppositeArc(a);

              (*_status)[v] = MATCHED;

              break;

            }

          }

        }

      }

    }

    ///

    /// \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)[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)[u] = _graph.direct(e, true);

          (*_status)[u] = MATCHED;

          Node v = _graph.v(e);

          if ((*_matching)[v] != INVALID) return false;

          (*_matching)[v] = _graph.direct(e, false);

          (*_status)[v] = MATCHED;

        }

      }

      return true;

    }

    ///

    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)[n] = EVEN;

          processSparse(n);

        }

      }

    }

    ///

    /// 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)[n] = EVEN;

          processDense(n);

        }

      }

    }

    ///

    void run() {

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

        greedyInit();

        startSparse();

      } else {

        init();

        startDense();

      }

    }

    /// @}

    /// @{

    ///

    int matchingSize() const {

      int size = 0;

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

        if ((*_matching)[n] != INVALID) {

          ++size;

        }

      }

      return size / 2;

    }

    ///

    bool matching(const Edge& edge) const {

      return edge == (*_matching)[_graph.u(edge)];

    }

    ///

    Arc matching(const Node& n) const {

      return (*_matching)[n];

    }

    ///

    Node mate(const Node& n) const {

      return (*_matching)[n] != INVALID ?

        _graph.target((*_matching)[n]) : INVALID;

    }

    /// @}

    /// @{

    /// decomposition.

    ///

    Status decomposition(const Node& n) const {

      return (*_status)[n];

    }

    ///

    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.

  ///

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

  ///

  template <typename GR,

            typename WM = typename GR::template EdgeMap<int> >

  class MaxWeightedMatching {

  public:

    typedef GR Graph;

    typedef WM WeightMap;

    typedef typename WeightMap::Value Value;

    /// \brief Scaling factor for dual solution

    ///

    /// according to the value type.

    static const int dualScale =

      std::numeric_limits<Value>::is_integer ? 4 : 1;

  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;

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

    }

    /// The simplest way to execute the algorithm is to use the

    ///@{

    /// \brief Initialize the algorithm

    ///

    void init() {

      createStructures();

      for (ArcIt e(_graph); e != INVALID; ++e) {

        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;

      }

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

        (*_delta1_index)[n] = _delta1->PRE_HEAP;

      }

      for (EdgeIt e(_graph); e != INVALID; ++e) {

        (*_delta3_index)[e] = _delta3->PRE_HEAP;

      }

      for (int i = 0; i < _blossom_num; ++i) {

        (*_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_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);

        }

      }

    }

    ///

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

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

    }

    ///

    ///

    /// \note mwm.run() is just a shortcut of the following code.

    /// \code

    ///   mwm.init();

    ///   mwm.start();

    /// \endcode

    void run() {

      init();

      start();

    }

    /// @}

    /// @{

    ///

    Value matchingValue() const {

      Value sum = 0;

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

        if ((*_matching)[n] != INVALID) {

          sum += _weight[(*_matching)[n]];

        }

      }

      return sum /= 2;

    }

    ///

    int matchingSize() const {

      int num = 0;

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

        if ((*_matching)[n] != INVALID) {

          ++num;

        }

      }

      return num /= 2;

    }

    ///

    bool matching(const Edge& edge) const {

      return edge == (*_matching)[_graph.u(edge)];

    }

    ///

    Arc matching(const Node& node) const {

      return (*_matching)[node];

    }

    ///

    Node mate(const Node& node) const {

      return (*_matching)[node] != INVALID ?

        _graph.target((*_matching)[node]) : INVALID;

    }

    /// @}

    /// @{

    ///

    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;

    }

    ///

    Value nodeValue(const Node& n) const {

      return (*_node_potential)[n];

    }

    ///

    /// \see BlossomIt

    int blossomNum() const {

      return _blossom_potential.size();

    }

    ///

    int blossomSize(int k) const {

      return _blossom_potential[k].end - _blossom_potential[k].begin;

    }

    ///

    Value blossomValue(int k) const {

      return _blossom_potential[k].value;

    }

    ///

    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;

      }

      ///

      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.

  ///

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

  ///

  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;