RhodeCode

\f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function. The minimum

cut is the \f$X\f$ solution of the next optimization problem:

\f[ \min_{X \subset V, X\not\in \{\emptyset, V\}}

    \sum_{uv\in A, u\in X, v\not\in X}cap(uv) \f]

LEMON contains several algorithms related to minimum cut problems:

- \ref HaoOrlin "Hao-Orlin algorithm" for calculating minimum cut

  in directed graphs.

- \ref NagamochiIbaraki "Nagamochi-Ibaraki algorithm" for

  calculating minimum cut in undirected graphs.

- \ref GomoryHu "Gomory-Hu tree computation" for calculating

  all-pairs minimum cut in undirected graphs.

If you want to find minimum cut just between two distinict nodes,

see the \ref max_flow "maximum flow problem".

*/

/**

@defgroup graph_properties Connectivity and Other Graph Properties

@ingroup algs

\brief Algorithms for discovering the graph properties

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

  Edmond's blossom shrinking algorithm for calculating maximum weighted

  perfect matching in general graphs.

\image html bipartite_matching.png

\image latex bipartite_matching.eps "Bipartite Matching" width=\textwidth

*/

/**

@defgroup spantree Minimum Spanning Tree Algorithms

@ingroup algs

\brief Algorithms for finding a minimum cost spanning tree in a graph.

This group contains the algorithms for finding a minimum cost spanning

tree in a graph.

*/

/**

@defgroup auxalg Auxiliary Algorithms

@ingroup algs

\brief Auxiliary algorithms implemented in LEMON.

This group contains some algorithms implemented in LEMON

in order to make it easier to implement complex algorithms.

*/

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

  ///

  ///

  ///

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

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

          }

          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;

      {

        Node node = _graph.v(e);

        Arc arc = _graph.direct(e, false);

        Node base = (*_blossom_rep)[_blossom_set->find(node)];

        while (base != nca) {

          (*_ear)[node] = arc;

          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;

      Arc arc = (*_matching)[even];

      (*_matching)[even] = a;

      while (arc != INVALID) {

        odd = _graph.target(arc);

        arc = (*_ear)[odd];

        even = _graph.target(arc);

        (*_matching)[odd] = arc;

        arc = (*_matching)[even];

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

      }

      for (typename TreeSet::ItemIt it(*_tree_set, tree);

           it != INVALID; ++it) {

        if ((*_status)[it] == ODD) {

          (*_status)[it] = MATCHED;

        } else {

          int blossom = _blossom_set->find(it);

          for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);

               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;

    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;

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

    }

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

        }

      }

    }

    ///

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

    }

    ///

    ///

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

    }

    /// @}

    /// @{

    ///