RhodeCode

/* -*- mode: C++; indent-tabs-mode: nil; -*-

 *

 * This file is a part of LEMON, a generic C++ optimization library.

 *

 * Copyright (C) 2003-2009

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

 * (Egervary Research Group on Combinatorial Optimization, EGRES).

 *

 * Permission to use, modify and distribute this software is granted

 * provided that this copyright notice appears in all copies. For

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

#include <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 Maximum cardinality matching in general graphs

  ///

  /// This class implements Edmonds' alternating forest matching algorithm

  /// for finding a maximum cardinality matching in a general undirected graph.

  /// It can be started from an arbitrary initial matching 

      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 {

    };

    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;

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

      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;

           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 {

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

                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 {

            typename std::map<int, Arc>::iterator it =

              (*_node_data)[vi].heap_index.find(tree);

            if (it != (*_node_data)[vi].heap_index.end()) {

              (*_node_data)[vi].heap.erase(it->second);

              (*_node_data)[vi].heap_index.erase(it);

              if ((*_node_data)[vi].heap.empty()) {

                _blossom_set->increase(v, std::numeric_limits<Value>::max());

              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {

                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());

              }

              if ((*_blossom_data)[vb].status == MATCHED) {

                if (_blossom_set->classPrio(vb) ==

                    std::numeric_limits<Value>::max()) {

                  _delta2->erase(vb);

                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -

                           (*_blossom_data)[vb].offset) {

                  _delta2->increase(vb, _blossom_set->classPrio(vb) -

                                   (*_blossom_data)[vb].offset);

                }

              }

           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 +=

          2 * _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 {

            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 unmatchNode(const Node& node) {

      int blossom = _blossom_set->find(node);

      int tree = _tree_set->find(blossom);

      alternatePath(blossom, tree);

      destroyTree(tree);

      (*_blossom_data)[blossom].base = node;

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

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

        } else {

          Node base = (*_blossom_data)[blossoms[i]].base;

          extractBlossom(blossoms[i], base, INVALID);

        }

      }

    }

  public:

    /// \brief Constructor

    ///

    /// Constructor.

    MaxWeightedMatching(const Graph& graph, const WeightMap& weight)

    /// \brief Start the algorithm

    ///

    /// This function starts the algorithm.

    ///

    /// \pre \ref init() must be called before using this function.

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

        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }

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

          }

          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 {

              if (left_tree == right_tree) {

                shrinkOnEdge(e, left_tree);

              } else {

                augmentOnEdge(e);

                unmatched -= 2;

              }

            }

          } break;

        case D4:

          splitBlossom(_delta4->top());

          break;

        }

      }

      extractMatching();

    }

    /// \brief Run the algorithm.

    ///

    /// This method runs the \c %MaxWeightedMatching algorithm.

    ///

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

    /// \code

    ///   mwm.init();

    }

    /// @}

    /// \name Primal Solution

    /// Functions to get the primal solution, i.e. the maximum weighted 

    /// matching.\n

    /// Either \ref run() or \ref start() function should be called before

    /// using them.

    /// @{

    /// \brief Return the weight of the matching.

    ///

    /// This function returns the weight of the found matching.

    ///

    /// \pre Either run() or start() must be called before using this function.

    Value matchingWeight() const {

      Value sum = 0;

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

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

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

        }

      }

    }

    /// \brief Return the size (cardinality) of the matching.

    ///

    /// This function returns the size (cardinality) of the found matching.

    ///

    /// \pre Either run() or start() must be called before using this function.

    int matchingSize() const {

      int num = 0;

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

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

          ++num;

        }

      }

      return num /= 2;

    }

    /// \brief Return \c true if the given edge is in the matching.

    ///

    /// This function returns \c true if the given edge is in the found 

    /// matching.

    ///

    /// \pre Either run() or start() must be called before using this function.

    bool matching(const Edge& edge) const {

        _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 (!_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() {

      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 (_delta2) {

        delete _delta2_index;

      }

    }

    /// \brief Start the algorithm

    ///

    /// This function starts the algorithm.

    ///

    /// \pre \ref init() must be called before using this function.

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

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

    }

    /// @}

    /// \name Primal Solution

    /// Functions to get the primal solution, i.e. the maximum weighted 

    /// perfect matching.\n

    /// Either \ref run() or \ref start() function should be called before

    /// using them.

    /// @{

    /// \brief Return the weight of the matching.

    ///

    /// This function returns the weight of the found matching.

    ///

    /// \pre Either run() or start() must be called before using this function.

    Value matchingWeight() const {

      Value sum = 0;

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

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

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

        }

      }

    }

    /// \brief Return \c true if the given edge is in the matching.

    ///

    /// This function returns \c true if the given edge is in the found 

    /// matching.

    ///

    /// \pre Either run() or start() must be called before using this function.

    bool matching(const Edge& edge) const {

      return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;

    }

    /// \brief Return the matching arc (or edge) incident to the given node.

    ///

    /// This function returns the matching arc (or edge) incident to the

    /// given node in the found matching or \c INVALID if the node is 

    /// not covered by the matching.

    ///

    /// \pre Either run() or start() must be called before using this function.

    Arc matching(const Node& node) const {

      return (*_matching)[node];

    }

    /// \brief Return a const reference to the matching map.

      }

      /// \brief Validity checking

      ///

      /// This function checks whether the iterator is invalid.

      bool operator==(Invalid) const { return _index == _last; }

      /// \brief Validity checking

      ///

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