[Lemon-commits] deba: r3426 - in lemon/trunk: doc lemon
Lemon SVN
svn at lemon.cs.elte.hu
Fri Dec 28 12:00:56 CET 2007
Author: deba
Date: Fri Dec 28 12:00:51 2007
New Revision: 3426
Modified:
lemon/trunk/doc/groups.dox
lemon/trunk/lemon/bin_heap.h
lemon/trunk/lemon/max_matching.h
lemon/trunk/lemon/unionfind.h
Log:
Edmond's Blossom shrinking algroithm:
MaxWeightedMatching
MaxWeightedPerfectMatching
Modified: lemon/trunk/doc/groups.dox
==============================================================================
--- lemon/trunk/doc/groups.dox (original)
+++ lemon/trunk/doc/groups.dox Fri Dec 28 12:00:51 2007
@@ -318,8 +318,37 @@
\brief This group describes the algorithms
for find matchings in graphs and bipartite graphs.
-This group provides some algorithm objects and function
-to calculate matchings in graphs and bipartite graphs.
+This group provides some algorithm objects and function to calculate
+matchings in graphs and bipartite graphs. The general matching problem is
+finding a subset of the edges which does not shares common endpoints.
+
+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 the finding maximum cardinality, maximum weight or minimum cost
+matching. The search can be constrained to find perfect or
+maximum cardinality matching.
+
+Lemon contains the next algorithms:
+- \ref lemon::MaxBipartiteMatching "MaxBipartiteMatching" Hopcroft-Karp
+ augmenting path algorithm for calculate maximum cardinality matching in
+ bipartite graphs
+- \ref lemon::PrBipartiteMatching "PrBipartiteMatching" Push-Relabel
+ algorithm for calculate maximum cardinality matching in bipartite graphs
+- \ref lemon::MaxWeightedBipartiteMatching "MaxWeightedBipartiteMatching"
+ Successive shortest path algorithm for calculate maximum weighted matching
+ and maximum weighted bipartite matching in bipartite graph
+- \ref lemon::MinCostMaxBipartiteMatching "MinCostMaxBipartiteMatching"
+ Successive shortest path algorithm for calculate minimum cost maximum
+ matching in bipartite graph
+- \ref lemon::MaxMatching "MaxMatching" Edmond's blossom shrinking algorithm
+ for calculate maximum cardinality matching in general graph
+- \ref lemon::MaxWeightedMatching "MaxWeightedMatching" Edmond's blossom
+ shrinking algorithm for calculate maximum weighted matching in general
+ graph
+- \ref lemon::MaxWeightedPerfectMatching "MaxWeightedPerfectMatching"
+ Edmond's blossom shrinking algorithm for calculate maximum weighted
+ perfect matching in general graph
\image html bipartite_matching.png
\image latex bipartite_matching.eps "Bipartite Matching" width=\textwidth
Modified: lemon/trunk/lemon/bin_heap.h
==============================================================================
--- lemon/trunk/lemon/bin_heap.h (original)
+++ lemon/trunk/lemon/bin_heap.h Fri Dec 28 12:00:51 2007
@@ -52,10 +52,15 @@
class BinHeap {
public:
+ ///\e
typedef _ItemIntMap ItemIntMap;
+ ///\e
typedef _Prio Prio;
+ ///\e
typedef typename ItemIntMap::Key Item;
+ ///\e
typedef std::pair<Item,Prio> Pair;
+ ///\e
typedef _Compare Compare;
/// \brief Type to represent the items states.
@@ -321,6 +326,19 @@
}
}
+ /// \brief Replaces an item in the heap.
+ ///
+ /// The \c i item is replaced with \c j item. The \c i item should
+ /// be in the heap, while the \c j should be out of the heap. The
+ /// \c i item will out of the heap and \c j will be in the heap
+ /// with the same prioriority as prevoiusly the \c i item.
+ void replace(const Item& i, const Item& j) {
+ int idx = iim[i];
+ iim.set(i, iim[j]);
+ iim.set(j, idx);
+ data[idx].first = j;
+ }
+
}; // class BinHeap
} // namespace lemon
Modified: lemon/trunk/lemon/max_matching.h
==============================================================================
--- lemon/trunk/lemon/max_matching.h (original)
+++ lemon/trunk/lemon/max_matching.h Fri Dec 28 12:00:51 2007
@@ -19,10 +19,14 @@
#ifndef LEMON_MAX_MATCHING_H
#define LEMON_MAX_MATCHING_H
+#include <vector>
#include <queue>
+#include <set>
+
#include <lemon/bits/invalid.h>
#include <lemon/unionfind.h>
#include <lemon/graph_utils.h>
+#include <lemon/bin_heap.h>
///\ingroup matching
///\file
@@ -31,7 +35,7 @@
namespace lemon {
///\ingroup matching
-
+ ///
///\brief Edmonds' alternating forest maximum matching algorithm.
///
///This class provides Edmonds' alternating forest matching
@@ -618,6 +622,2500 @@
}
};
+
+ /// \ingroup matching
+ ///
+ /// \brief Weighted matching in general undirected 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 endpoints. The problem can be formulated with
+ /// the next 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 endpoints 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 next:
+ /// \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".
+ ///
+ /// \author Balazs Dezso
+ template <typename _UGraph,
+ typename _WeightMap = typename _UGraph::template UEdgeMap<int> >
+ class MaxWeightedMatching {
+ public:
+
+ typedef _UGraph UGraph;
+ typedef _WeightMap 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 UGraph::template NodeMap<typename UGraph::Edge>
+ MatchingMap;
+
+ private:
+
+ UGRAPH_TYPEDEFS(typename UGraph);
+
+ typedef typename UGraph::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 UGraph& _ugraph;
+ const WeightMap& _weight;
+
+ MatchingMap* _matching;
+
+ NodePotential* _node_potential;
+
+ BlossomPotential _blossom_potential;
+ BlossomNodeList _blossom_node_list;
+
+ int _node_num;
+ int _blossom_num;
+
+ typedef typename UGraph::template NodeMap<int> NodeIntMap;
+ typedef typename UGraph::template EdgeMap<int> EdgeIntMap;
+ typedef typename UGraph::template UEdgeMap<int> UEdgeIntMap;
+ typedef IntegerMap<int> IntIntMap;
+
+ enum Status {
+ EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
+ };
+
+ typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
+ struct BlossomData {
+ int tree;
+ Status status;
+ Edge pred, next;
+ Value pot, offset;
+ Node base;
+ };
+
+ NodeIntMap *_blossom_index;
+ BlossomSet *_blossom_set;
+ IntegerMap<BlossomData>* _blossom_data;
+
+ NodeIntMap *_node_index;
+ EdgeIntMap *_node_heap_index;
+
+ struct NodeData {
+
+ NodeData(EdgeIntMap& node_heap_index)
+ : heap(node_heap_index) {}
+
+ int blossom;
+ Value pot;
+ BinHeap<Value, EdgeIntMap> heap;
+ std::map<int, Edge> heap_index;
+
+ int tree;
+ };
+
+ IntegerMap<NodeData>* _node_data;
+
+ typedef ExtendFindEnum<IntIntMap> TreeSet;
+
+ IntIntMap *_tree_set_index;
+ TreeSet *_tree_set;
+
+ NodeIntMap *_delta1_index;
+ BinHeap<Value, NodeIntMap> *_delta1;
+
+ IntIntMap *_delta2_index;
+ BinHeap<Value, IntIntMap> *_delta2;
+
+ UEdgeIntMap *_delta3_index;
+ BinHeap<Value, UEdgeIntMap> *_delta3;
+
+ IntIntMap *_delta4_index;
+ BinHeap<Value, IntIntMap> *_delta4;
+
+ Value _delta_sum;
+
+ void createStructures() {
+ _node_num = countNodes(_ugraph);
+ _blossom_num = _node_num * 3 / 2;
+
+ if (!_matching) {
+ _matching = new MatchingMap(_ugraph);
+ }
+ if (!_node_potential) {
+ _node_potential = new NodePotential(_ugraph);
+ }
+ if (!_blossom_set) {
+ _blossom_index = new NodeIntMap(_ugraph);
+ _blossom_set = new BlossomSet(*_blossom_index);
+ _blossom_data = new IntegerMap<BlossomData>(_blossom_num);
+ }
+
+ if (!_node_index) {
+ _node_index = new NodeIntMap(_ugraph);
+ _node_heap_index = new EdgeIntMap(_ugraph);
+ _node_data = new IntegerMap<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 NodeIntMap(_ugraph);
+ _delta1 = new BinHeap<Value, NodeIntMap>(*_delta1_index);
+ }
+ if (!_delta2) {
+ _delta2_index = new IntIntMap(_blossom_num);
+ _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
+ }
+ if (!_delta3) {
+ _delta3_index = new UEdgeIntMap(_ugraph);
+ _delta3 = new BinHeap<Value, UEdgeIntMap>(*_delta3_index);
+ }
+ if (!_delta4) {
+ _delta4_index = new IntIntMap(_blossom_num);
+ _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
+ }
+ }
+
+ void destroyStructures() {
+ _node_num = countNodes(_ugraph);
+ _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 -=
+ 2 * (_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 (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+ Node v = _ugraph.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, Edge>::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 matchedToOdd(int blossom) {
+ if (_delta2->state(blossom) == _delta2->IN_HEAP) {
+ _delta2->erase(blossom);
+ }
+ (*_blossom_data)[blossom].offset += _delta_sum;
+ if (!_blossom_set->trivial(blossom)) {
+ _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
+ (*_blossom_data)[blossom].offset);
+ }
+ }
+
+ void evenToMatched(int blossom, int tree) {
+ if (!_blossom_set->trivial(blossom)) {
+ (*_blossom_data)[blossom].pot += 2 * _delta_sum;
+ }
+
+ for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
+ n != INVALID; ++n) {
+ int ni = (*_node_index)[n];
+ (*_node_data)[ni].pot -= _delta_sum;
+
+ _delta1->erase(n);
+
+ for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+ Node v = _ugraph.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 (vb == blossom) {
+ if (_delta3->state(e) == _delta3->IN_HEAP) {
+ _delta3->erase(e);
+ }
+ } else if ((*_blossom_data)[vb].status == EVEN) {
+
+ if (_delta3->state(e) == _delta3->IN_HEAP) {
+ _delta3->erase(e);
+ }
+
+ int vt = _tree_set->find(vb);
+
+ if (vt != tree) {
+
+ Edge r = _ugraph.oppositeEdge(e);
+
+ typename std::map<int, Edge>::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);
+ }
+ } else {
+
+ typename std::map<int, Edge>::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);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ void oddToMatched(int blossom) {
+ (*_blossom_data)[blossom].offset -= _delta_sum;
+
+ if (_blossom_set->classPrio(blossom) !=
+ std::numeric_limits<Value>::max()) {
+ _delta2->push(blossom, _blossom_set->classPrio(blossom) -
+ (*_blossom_data)[blossom].offset);
+ }
+
+ if (!_blossom_set->trivial(blossom)) {
+ _delta4->erase(blossom);
+ }
+ }
+
+ void oddToEven(int blossom, int tree) {
+ if (!_blossom_set->trivial(blossom)) {
+ _delta4->erase(blossom);
+ (*_blossom_data)[blossom].pot -=
+ 2 * (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 +=
+ 2 * _delta_sum - (*_blossom_data)[blossom].offset;
+
+ _delta1->push(n, (*_node_data)[ni].pot);
+
+ for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+ Node v = _ugraph.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, Edge>::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 matchedToUnmatched(int blossom) {
+ if (_delta2->state(blossom) == _delta2->IN_HEAP) {
+ _delta2->erase(blossom);
+ }
+
+ 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();
+
+ for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+ Node v = _ugraph.target(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) {
+ _delta3->push(e, rw);
+ }
+ }
+ }
+ }
+ }
+
+ void unmatchedToMatched(int blossom) {
+ for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
+ n != INVALID; ++n) {
+ int ni = (*_node_index)[n];
+
+ for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+ Node v = _ugraph.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 (vb == blossom) {
+ if (_delta3->state(e) == _delta3->IN_HEAP) {
+ _delta3->erase(e);
+ }
+ } else if ((*_blossom_data)[vb].status == EVEN) {
+
+ if (_delta3->state(e) == _delta3->IN_HEAP) {
+ _delta3->erase(e);
+ }
+
+ int vt = _tree_set->find(vb);
+
+ Edge r = _ugraph.oppositeEdge(e);
+
+ typename std::map<int, Edge>::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(_ugraph.target((*_blossom_data)[even].pred));
+ (*_blossom_data)[odd].status = MATCHED;
+ oddToMatched(odd);
+ (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
+
+ even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].pred));
+ (*_blossom_data)[even].status = MATCHED;
+ evenToMatched(even, tree);
+ (*_blossom_data)[even].next =
+ _ugraph.oppositeEdge((*_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 UEdge& edge) {
+
+ int left = _blossom_set->find(_ugraph.source(edge));
+ int right = _blossom_set->find(_ugraph.target(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 = _ugraph.direct(edge, true);
+ (*_blossom_data)[right].next = _ugraph.direct(edge, false);
+ }
+
+ void extendOnEdge(const Edge& edge) {
+ int base = _blossom_set->find(_ugraph.target(edge));
+ int tree = _tree_set->find(base);
+
+ int odd = _blossom_set->find(_ugraph.source(edge));
+ _tree_set->insert(odd, tree);
+ (*_blossom_data)[odd].status = ODD;
+ matchedToOdd(odd);
+ (*_blossom_data)[odd].pred = edge;
+
+ int even = _blossom_set->find(_ugraph.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 UEdge& uedge, int tree) {
+ int nca = -1;
+ std::vector<int> left_path, right_path;
+
+ {
+ std::set<int> left_set, right_set;
+ int left = _blossom_set->find(_ugraph.source(uedge));
+ left_path.push_back(left);
+ left_set.insert(left);
+
+ int right = _blossom_set->find(_ugraph.target(uedge));
+ right_path.push_back(right);
+ right_set.insert(right);
+
+ while (true) {
+
+ if ((*_blossom_data)[left].pred == INVALID) break;
+
+ left =
+ _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
+ left_path.push_back(left);
+ left =
+ _blossom_set->find(_ugraph.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(_ugraph.target((*_blossom_data)[right].pred));
+ right_path.push_back(right);
+ right =
+ _blossom_set->find(_ugraph.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(_ugraph.target((*_blossom_data)[nca].pred));
+ right_path.push_back(nca);
+ nca =
+ _blossom_set->find(_ugraph.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(_ugraph.target((*_blossom_data)[nca].pred));
+ left_path.push_back(nca);
+ nca =
+ _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
+ left_path.push_back(nca);
+ }
+ }
+ }
+ }
+
+ std::vector<int> subblossoms;
+ Edge prev;
+
+ prev = _ugraph.direct(uedge, 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 = _ugraph.oppositeEdge((*_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) {
+ Edge next = (*_blossom_data)[blossom].next;
+ Edge 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(_ugraph.source(pred));
+ int d = _blossom_set->find(_ugraph.source(next));
+
+ int ib, id;
+ 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 =
+ _ugraph.oppositeEdge((*_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 =
+ _ugraph.oppositeEdge((*_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 =
+ _ugraph.oppositeEdge((*_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 =
+ _ugraph.oppositeEdge((*_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 =
+ _ugraph.oppositeEdge((*_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 Edge& matching) {
+ if (_blossom_set->trivial(blossom)) {
+ int bi = (*_node_index)[base];
+ Value pot = (*_node_data)[bi].pot;
+
+ _matching->set(base, matching);
+ _blossom_node_list.push_back(base);
+ _node_potential->set(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()];
+
+ Edge m = (*_blossom_data)[tb].next;
+ extractBlossom(sb, _ugraph.target(m), _ugraph.oppositeEdge(m));
+ extractBlossom(tb, _ugraph.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;
+ }
+
+ Edge matching = (*_blossom_data)[blossoms[i]].next;
+ Node base = _ugraph.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 UGraph& ugraph, const WeightMap& weight)
+ : _ugraph(ugraph), _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 member
+ /// \c run() member function.
+
+ ///@{
+
+ /// \brief Initialize the algorithm
+ ///
+ /// Initialize the algorithm
+ void init() {
+ createStructures();
+
+ for (EdgeIt e(_ugraph); e != INVALID; ++e) {
+ _node_heap_index->set(e, BinHeap<Value, EdgeIntMap>::PRE_HEAP);
+ }
+ for (NodeIt n(_ugraph); n != INVALID; ++n) {
+ _delta1_index->set(n, _delta1->PRE_HEAP);
+ }
+ for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
+ _delta3_index->set(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);
+ }
+
+ int index = 0;
+ for (NodeIt n(_ugraph); n != INVALID; ++n) {
+ Value max = 0;
+ for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+ if (_ugraph.target(e) == n) continue;
+ if ((dualScale * _weight[e]) / 2 > max) {
+ max = (dualScale * _weight[e]) / 2;
+ }
+ }
+ _node_index->set(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 (UEdgeIt e(_ugraph); e != INVALID; ++e) {
+ int si = (*_node_index)[_ugraph.source(e)];
+ int ti = (*_node_index)[_ugraph.target(e)];
+ if (_ugraph.source(e) != _ugraph.target(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);
+ Edge e = (*_node_data)[(*_node_index)[n]].heap.top();
+ extendOnEdge(e);
+ }
+ break;
+ case D3:
+ {
+ UEdge e = _delta3->top();
+
+ int left_blossom = _blossom_set->find(_ugraph.source(e));
+ int right_blossom = _blossom_set->find(_ugraph.target(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 for 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(_ugraph); 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 UEdge& edge) const {
+ return (*_matching)[_ugraph.source(edge)] == _ugraph.direct(edge, true);
+ }
+
+ /// \brief Returns the incident matching edge.
+ ///
+ /// Returns the incident matching edge from given node. If the
+ /// node is not matched then it gives back \c INVALID.
+ Edge matching(const Node& node) const {
+ return (*_matching)[node];
+ }
+
+ /// \brief Returns the mate of the node.
+ ///
+ /// Returns the adjancent node in a mathcing edge. If the node is
+ /// not matched then it gives back \c INVALID.
+ Node mate(const Node& node) const {
+ return (*_matching)[node] != INVALID ?
+ _ugraph.target((*_matching)[node]) : INVALID;
+ }
+
+ /// @}
+
+ /// \name Dual solution
+ /// Functions for 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(_ugraph); 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 Lemon iterator for get the items of the blossom.
+ ///
+ /// Lemon iterator for get the nodes of the blossom. This class
+ /// provides a common style lemon iterator which gives back a
+ /// subset of the nodes.
+ class BlossomIt {
+ public:
+
+ /// \brief Constructor.
+ ///
+ /// Constructor for 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 Invalid constructor.
+ ///
+ /// Invalid constructor.
+ BlossomIt(Invalid) : _index(-1) {}
+
+ /// \brief Conversion to node.
+ ///
+ /// Conversion to node.
+ operator Node() const {
+ return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
+ }
+
+ /// \brief Increment operator.
+ ///
+ /// Increment operator.
+ BlossomIt& operator++() {
+ ++_index;
+ if (_index == _last) {
+ _index = -1;
+ }
+ return *this;
+ }
+
+ bool operator==(const BlossomIt& it) const {
+ return _index == it._index;
+ }
+ bool operator!=(const BlossomIt& it) const {
+ return _index != it._index;
+ }
+
+ private:
+ const MaxWeightedMatching* _algorithm;
+ int _last;
+ int _index;
+ };
+
+ /// @}
+
+ };
+
+ /// \ingroup matching
+ ///
+ /// \brief Weighted perfect matching in general undirected graphs
+ ///
+ /// This class provides an efficient implementation of Edmond's
+ /// maximum weighted perfecr 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 endpoints and covers all nodes. The problem
+ /// can be formulated with the next 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 endpoints 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 next:
+ /// \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".
+ ///
+ /// \author Balazs Dezso
+ template <typename _UGraph,
+ typename _WeightMap = typename _UGraph::template UEdgeMap<int> >
+ class MaxWeightedPerfectMatching {
+ public:
+
+ typedef _UGraph UGraph;
+ typedef _WeightMap 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 UGraph::template NodeMap<typename UGraph::Edge>
+ MatchingMap;
+
+ private:
+
+ UGRAPH_TYPEDEFS(typename UGraph);
+
+ typedef typename UGraph::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 UGraph& _ugraph;
+ const WeightMap& _weight;
+
+ MatchingMap* _matching;
+
+ NodePotential* _node_potential;
+
+ BlossomPotential _blossom_potential;
+ BlossomNodeList _blossom_node_list;
+
+ int _node_num;
+ int _blossom_num;
+
+ typedef typename UGraph::template NodeMap<int> NodeIntMap;
+ typedef typename UGraph::template EdgeMap<int> EdgeIntMap;
+ typedef typename UGraph::template UEdgeMap<int> UEdgeIntMap;
+ typedef IntegerMap<int> IntIntMap;
+
+ enum Status {
+ EVEN = -1, MATCHED = 0, ODD = 1
+ };
+
+ typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
+ struct BlossomData {
+ int tree;
+ Status status;
+ Edge pred, next;
+ Value pot, offset;
+ };
+
+ NodeIntMap *_blossom_index;
+ BlossomSet *_blossom_set;
+ IntegerMap<BlossomData>* _blossom_data;
+
+ NodeIntMap *_node_index;
+ EdgeIntMap *_node_heap_index;
+
+ struct NodeData {
+
+ NodeData(EdgeIntMap& node_heap_index)
+ : heap(node_heap_index) {}
+
+ int blossom;
+ Value pot;
+ BinHeap<Value, EdgeIntMap> heap;
+ std::map<int, Edge> heap_index;
+
+ int tree;
+ };
+
+ IntegerMap<NodeData>* _node_data;
+
+ typedef ExtendFindEnum<IntIntMap> TreeSet;
+
+ IntIntMap *_tree_set_index;
+ TreeSet *_tree_set;
+
+ IntIntMap *_delta2_index;
+ BinHeap<Value, IntIntMap> *_delta2;
+
+ UEdgeIntMap *_delta3_index;
+ BinHeap<Value, UEdgeIntMap> *_delta3;
+
+ IntIntMap *_delta4_index;
+ BinHeap<Value, IntIntMap> *_delta4;
+
+ Value _delta_sum;
+
+ void createStructures() {
+ _node_num = countNodes(_ugraph);
+ _blossom_num = _node_num * 3 / 2;
+
+ if (!_matching) {
+ _matching = new MatchingMap(_ugraph);
+ }
+ if (!_node_potential) {
+ _node_potential = new NodePotential(_ugraph);
+ }
+ if (!_blossom_set) {
+ _blossom_index = new NodeIntMap(_ugraph);
+ _blossom_set = new BlossomSet(*_blossom_index);
+ _blossom_data = new IntegerMap<BlossomData>(_blossom_num);
+ }
+
+ if (!_node_index) {
+ _node_index = new NodeIntMap(_ugraph);
+ _node_heap_index = new EdgeIntMap(_ugraph);
+ _node_data = new IntegerMap<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 UEdgeIntMap(_ugraph);
+ _delta3 = new BinHeap<Value, UEdgeIntMap>(*_delta3_index);
+ }
+ if (!_delta4) {
+ _delta4_index = new IntIntMap(_blossom_num);
+ _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
+ }
+ }
+
+ void destroyStructures() {
+ _node_num = countNodes(_ugraph);
+ _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 (_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 -=
+ 2 * (_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;
+
+ for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+ Node v = _ugraph.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, Edge>::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 matchedToOdd(int blossom) {
+ if (_delta2->state(blossom) == _delta2->IN_HEAP) {
+ _delta2->erase(blossom);
+ }
+ (*_blossom_data)[blossom].offset += _delta_sum;
+ if (!_blossom_set->trivial(blossom)) {
+ _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
+ (*_blossom_data)[blossom].offset);
+ }
+ }
+
+ void evenToMatched(int blossom, int tree) {
+ if (!_blossom_set->trivial(blossom)) {
+ (*_blossom_data)[blossom].pot += 2 * _delta_sum;
+ }
+
+ for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
+ n != INVALID; ++n) {
+ int ni = (*_node_index)[n];
+ (*_node_data)[ni].pot -= _delta_sum;
+
+ for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+ Node v = _ugraph.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 (vb == blossom) {
+ if (_delta3->state(e) == _delta3->IN_HEAP) {
+ _delta3->erase(e);
+ }
+ } else if ((*_blossom_data)[vb].status == EVEN) {
+
+ if (_delta3->state(e) == _delta3->IN_HEAP) {
+ _delta3->erase(e);
+ }
+
+ int vt = _tree_set->find(vb);
+
+ if (vt != tree) {
+
+ Edge r = _ugraph.oppositeEdge(e);
+
+ typename std::map<int, Edge>::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 {
+
+ typename std::map<int, Edge>::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);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ void oddToMatched(int blossom) {
+ (*_blossom_data)[blossom].offset -= _delta_sum;
+
+ if (_blossom_set->classPrio(blossom) !=
+ std::numeric_limits<Value>::max()) {
+ _delta2->push(blossom, _blossom_set->classPrio(blossom) -
+ (*_blossom_data)[blossom].offset);
+ }
+
+ if (!_blossom_set->trivial(blossom)) {
+ _delta4->erase(blossom);
+ }
+ }
+
+ void oddToEven(int blossom, int tree) {
+ if (!_blossom_set->trivial(blossom)) {
+ _delta4->erase(blossom);
+ (*_blossom_data)[blossom].pot -=
+ 2 * (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 +=
+ 2 * _delta_sum - (*_blossom_data)[blossom].offset;
+
+ for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+ Node v = _ugraph.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, Edge>::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(_ugraph.target((*_blossom_data)[even].pred));
+ (*_blossom_data)[odd].status = MATCHED;
+ oddToMatched(odd);
+ (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
+
+ even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].pred));
+ (*_blossom_data)[even].status = MATCHED;
+ evenToMatched(even, tree);
+ (*_blossom_data)[even].next =
+ _ugraph.oppositeEdge((*_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 UEdge& edge) {
+
+ int left = _blossom_set->find(_ugraph.source(edge));
+ int right = _blossom_set->find(_ugraph.target(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 = _ugraph.direct(edge, true);
+ (*_blossom_data)[right].next = _ugraph.direct(edge, false);
+ }
+
+ void extendOnEdge(const Edge& edge) {
+ int base = _blossom_set->find(_ugraph.target(edge));
+ int tree = _tree_set->find(base);
+
+ int odd = _blossom_set->find(_ugraph.source(edge));
+ _tree_set->insert(odd, tree);
+ (*_blossom_data)[odd].status = ODD;
+ matchedToOdd(odd);
+ (*_blossom_data)[odd].pred = edge;
+
+ int even = _blossom_set->find(_ugraph.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 UEdge& uedge, int tree) {
+ int nca = -1;
+ std::vector<int> left_path, right_path;
+
+ {
+ std::set<int> left_set, right_set;
+ int left = _blossom_set->find(_ugraph.source(uedge));
+ left_path.push_back(left);
+ left_set.insert(left);
+
+ int right = _blossom_set->find(_ugraph.target(uedge));
+ right_path.push_back(right);
+ right_set.insert(right);
+
+ while (true) {
+
+ if ((*_blossom_data)[left].pred == INVALID) break;
+
+ left =
+ _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
+ left_path.push_back(left);
+ left =
+ _blossom_set->find(_ugraph.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(_ugraph.target((*_blossom_data)[right].pred));
+ right_path.push_back(right);
+ right =
+ _blossom_set->find(_ugraph.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(_ugraph.target((*_blossom_data)[nca].pred));
+ right_path.push_back(nca);
+ nca =
+ _blossom_set->find(_ugraph.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(_ugraph.target((*_blossom_data)[nca].pred));
+ left_path.push_back(nca);
+ nca =
+ _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
+ left_path.push_back(nca);
+ }
+ }
+ }
+ }
+
+ std::vector<int> subblossoms;
+ Edge prev;
+
+ prev = _ugraph.direct(uedge, 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 = _ugraph.oppositeEdge((*_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) {
+ Edge next = (*_blossom_data)[blossom].next;
+ Edge 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(_ugraph.source(pred));
+ int d = _blossom_set->find(_ugraph.source(next));
+
+ int ib, id;
+ 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 =
+ _ugraph.oppositeEdge((*_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 =
+ _ugraph.oppositeEdge((*_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 =
+ _ugraph.oppositeEdge((*_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 =
+ _ugraph.oppositeEdge((*_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 =
+ _ugraph.oppositeEdge((*_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 Edge& matching) {
+ if (_blossom_set->trivial(blossom)) {
+ int bi = (*_node_index)[base];
+ Value pot = (*_node_data)[bi].pot;
+
+ _matching->set(base, matching);
+ _blossom_node_list.push_back(base);
+ _node_potential->set(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()];
+
+ Edge m = (*_blossom_data)[tb].next;
+ extractBlossom(sb, _ugraph.target(m), _ugraph.oppositeEdge(m));
+ extractBlossom(tb, _ugraph.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;
+ }
+
+ Edge matching = (*_blossom_data)[blossoms[i]].next;
+ Node base = _ugraph.source(matching);
+ extractBlossom(blossoms[i], base, matching);
+ }
+ }
+
+ public:
+
+ /// \brief Constructor
+ ///
+ /// Constructor.
+ MaxWeightedPerfectMatching(const UGraph& ugraph, const WeightMap& weight)
+ : _ugraph(ugraph), _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 member
+ /// \c run() member function.
+
+ ///@{
+
+ /// \brief Initialize the algorithm
+ ///
+ /// Initialize the algorithm
+ void init() {
+ createStructures();
+
+ for (EdgeIt e(_ugraph); e != INVALID; ++e) {
+ _node_heap_index->set(e, BinHeap<Value, EdgeIntMap>::PRE_HEAP);
+ }
+ for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
+ _delta3_index->set(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);
+ }
+
+ int index = 0;
+ for (NodeIt n(_ugraph); n != INVALID; ++n) {
+ Value max = std::numeric_limits<Value>::min();
+ for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+ if (_ugraph.target(e) == n) continue;
+ if ((dualScale * _weight[e]) / 2 > max) {
+ max = (dualScale * _weight[e]) / 2;
+ }
+ }
+ _node_index->set(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 (UEdgeIt e(_ugraph); e != INVALID; ++e) {
+ int si = (*_node_index)[_ugraph.source(e)];
+ int ti = (*_node_index)[_ugraph.target(e)];
+ if (_ugraph.source(e) != _ugraph.target(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);
+ Edge e = (*_node_data)[(*_node_index)[n]].heap.top();
+ extendOnEdge(e);
+ }
+ break;
+ case D3:
+ {
+ UEdge e = _delta3->top();
+
+ int left_blossom = _blossom_set->find(_ugraph.source(e));
+ int right_blossom = _blossom_set->find(_ugraph.target(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 for 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(_ugraph); 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 UEdge& edge) const {
+ return (*_matching)[_ugraph.source(edge)] == _ugraph.direct(edge, true);
+ }
+
+ /// \brief Returns the incident matching edge.
+ ///
+ /// Returns the incident matching edge from given node.
+ Edge matching(const Node& node) const {
+ return (*_matching)[node];
+ }
+
+ /// \brief Returns the mate of the node.
+ ///
+ /// Returns the adjancent node in a mathcing edge.
+ Node mate(const Node& node) const {
+ return _ugraph.target((*_matching)[node]);
+ }
+
+ /// @}
+
+ /// \name Dual solution
+ /// Functions for 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(_ugraph); 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 Lemon iterator for get the items of the blossom.
+ ///
+ /// Lemon iterator for get the nodes of the blossom. This class
+ /// provides a common style lemon iterator which gives back a
+ /// subset of the nodes.
+ class BlossomIt {
+ public:
+
+ /// \brief Constructor.
+ ///
+ /// Constructor for 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 Invalid constructor.
+ ///
+ /// Invalid constructor.
+ BlossomIt(Invalid) : _index(-1) {}
+
+ /// \brief Conversion to node.
+ ///
+ /// Conversion to node.
+ operator Node() const {
+ return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
+ }
+
+ /// \brief Increment operator.
+ ///
+ /// Increment operator.
+ BlossomIt& operator++() {
+ ++_index;
+ if (_index == _last) {
+ _index = -1;
+ }
+ return *this;
+ }
+
+ bool operator==(const BlossomIt& it) const {
+ return _index == it._index;
+ }
+ bool operator!=(const BlossomIt& it) const {
+ return _index != it._index;
+ }
+
+ private:
+ const MaxWeightedPerfectMatching* _algorithm;
+ int _last;
+ int _index;
+ };
+
+ /// @}
+
+ };
+
} //END OF NAMESPACE LEMON
Modified: lemon/trunk/lemon/unionfind.h
==============================================================================
--- lemon/trunk/lemon/unionfind.h (original)
+++ lemon/trunk/lemon/unionfind.h Fri Dec 28 12:00:51 2007
@@ -924,6 +924,869 @@
};
+ /// \ingroup auxdat
+ ///
+ /// \brief A \e Union-Find data structure implementation which
+ /// is able to store a priority for each item and retrieve the minimum of
+ /// each class.
+ ///
+ /// A \e Union-Find data structure implementation which is able to
+ /// store a priority for each item and retrieve the minimum of each
+ /// class. In addition, it supports the joining and splitting the
+ /// components. If you don't need this feature then you makes
+ /// better to use the \ref UnionFind class which is more efficient.
+ ///
+ /// The union-find data strcuture based on a (2, 16)-tree with a
+ /// tournament minimum selection on the internal nodes. The insert
+ /// operation takes O(1), the find, set, decrease and increase takes
+ /// O(log(n)), where n is the number of nodes in the current
+ /// component. The complexity of join and split is O(log(n)*k),
+ /// where n is the sum of the number of the nodes and k is the
+ /// number of joined components or the number of the components
+ /// after the split.
+ ///
+ /// \pre You need to add all the elements by the \ref insert()
+ /// method.
+ ///
+ template <typename _Value, typename _ItemIntMap,
+ typename _Comp = std::less<_Value> >
+ class HeapUnionFind {
+ public:
+
+ typedef _Value Value;
+ typedef typename _ItemIntMap::Key Item;
+
+ typedef _ItemIntMap ItemIntMap;
+
+ typedef _Comp Comp;
+
+ private:
+
+ static const int cmax = 3;
+
+ ItemIntMap& index;
+
+ struct ClassNode {
+ int parent;
+ int depth;
+
+ int left, right;
+ int next, prev;
+ };
+
+ int first_class;
+ int first_free_class;
+ std::vector<ClassNode> classes;
+
+ int newClass() {
+ if (first_free_class < 0) {
+ int id = classes.size();
+ classes.push_back(ClassNode());
+ return id;
+ } else {
+ int id = first_free_class;
+ first_free_class = classes[id].next;
+ return id;
+ }
+ }
+
+ void deleteClass(int id) {
+ classes[id].next = first_free_class;
+ first_free_class = id;
+ }
+
+ struct ItemNode {
+ int parent;
+ Item item;
+ Value prio;
+ int next, prev;
+ int left, right;
+ int size;
+ };
+
+ int first_free_node;
+ std::vector<ItemNode> nodes;
+
+ int newNode() {
+ if (first_free_node < 0) {
+ int id = nodes.size();
+ nodes.push_back(ItemNode());
+ return id;
+ } else {
+ int id = first_free_node;
+ first_free_node = nodes[id].next;
+ return id;
+ }
+ }
+
+ void deleteNode(int id) {
+ nodes[id].next = first_free_node;
+ first_free_node = id;
+ }
+
+ Comp comp;
+
+ int findClass(int id) const {
+ int kd = id;
+ while (kd >= 0) {
+ kd = nodes[kd].parent;
+ }
+ return ~kd;
+ }
+
+ int leftNode(int id) const {
+ int kd = ~(classes[id].parent);
+ for (int i = 0; i < classes[id].depth; ++i) {
+ kd = nodes[kd].left;
+ }
+ return kd;
+ }
+
+ int nextNode(int id) const {
+ int depth = 0;
+ while (id >= 0 && nodes[id].next == -1) {
+ id = nodes[id].parent;
+ ++depth;
+ }
+ if (id < 0) {
+ return -1;
+ }
+ id = nodes[id].next;
+ while (depth--) {
+ id = nodes[id].left;
+ }
+ return id;
+ }
+
+
+ void setPrio(int id) {
+ int jd = nodes[id].left;
+ nodes[id].prio = nodes[jd].prio;
+ nodes[id].item = nodes[jd].item;
+ jd = nodes[jd].next;
+ while (jd != -1) {
+ if (comp(nodes[jd].prio, nodes[id].prio)) {
+ nodes[id].prio = nodes[jd].prio;
+ nodes[id].item = nodes[jd].item;
+ }
+ jd = nodes[jd].next;
+ }
+ }
+
+ void push(int id, int jd) {
+ nodes[id].size = 1;
+ nodes[id].left = nodes[id].right = jd;
+ nodes[jd].next = nodes[jd].prev = -1;
+ nodes[jd].parent = id;
+ }
+
+ void pushAfter(int id, int jd) {
+ int kd = nodes[id].parent;
+ if (nodes[id].next != -1) {
+ nodes[nodes[id].next].prev = jd;
+ if (kd >= 0) {
+ nodes[kd].size += 1;
+ }
+ } else {
+ if (kd >= 0) {
+ nodes[kd].right = jd;
+ nodes[kd].size += 1;
+ }
+ }
+ nodes[jd].next = nodes[id].next;
+ nodes[jd].prev = id;
+ nodes[id].next = jd;
+ nodes[jd].parent = kd;
+ }
+
+ void pushRight(int id, int jd) {
+ nodes[id].size += 1;
+ nodes[jd].prev = nodes[id].right;
+ nodes[jd].next = -1;
+ nodes[nodes[id].right].next = jd;
+ nodes[id].right = jd;
+ nodes[jd].parent = id;
+ }
+
+ void popRight(int id) {
+ nodes[id].size -= 1;
+ int jd = nodes[id].right;
+ nodes[nodes[jd].prev].next = -1;
+ nodes[id].right = nodes[jd].prev;
+ }
+
+ void splice(int id, int jd) {
+ nodes[id].size += nodes[jd].size;
+ nodes[nodes[id].right].next = nodes[jd].left;
+ nodes[nodes[jd].left].prev = nodes[id].right;
+ int kd = nodes[jd].left;
+ while (kd != -1) {
+ nodes[kd].parent = id;
+ kd = nodes[kd].next;
+ }
+ }
+
+ void split(int id, int jd) {
+ int kd = nodes[id].parent;
+ nodes[kd].right = nodes[id].prev;
+ nodes[nodes[id].prev].next = -1;
+
+ nodes[jd].left = id;
+ nodes[id].prev = -1;
+ int num = 0;
+ while (id != -1) {
+ nodes[id].parent = jd;
+ nodes[jd].right = id;
+ id = nodes[id].next;
+ ++num;
+ }
+ nodes[kd].size -= num;
+ nodes[jd].size = num;
+ }
+
+ void pushLeft(int id, int jd) {
+ nodes[id].size += 1;
+ nodes[jd].next = nodes[id].left;
+ nodes[jd].prev = -1;
+ nodes[nodes[id].left].prev = jd;
+ nodes[id].left = jd;
+ nodes[jd].parent = id;
+ }
+
+ void popLeft(int id) {
+ nodes[id].size -= 1;
+ int jd = nodes[id].left;
+ nodes[nodes[jd].next].prev = -1;
+ nodes[id].left = nodes[jd].next;
+ }
+
+ void repairLeft(int id) {
+ int jd = ~(classes[id].parent);
+ while (nodes[jd].left != -1) {
+ int kd = nodes[jd].left;
+ if (nodes[jd].size == 1) {
+ if (nodes[jd].parent < 0) {
+ classes[id].parent = ~kd;
+ classes[id].depth -= 1;
+ nodes[kd].parent = ~id;
+ deleteNode(jd);
+ jd = kd;
+ } else {
+ int pd = nodes[jd].parent;
+ if (nodes[nodes[jd].next].size < cmax) {
+ pushLeft(nodes[jd].next, nodes[jd].left);
+ if (less(nodes[jd].left, nodes[jd].next)) {
+ nodes[nodes[jd].next].prio = nodes[nodes[jd].left].prio;
+ nodes[nodes[jd].next].item = nodes[nodes[jd].left].item;
+ }
+ popLeft(pd);
+ deleteNode(jd);
+ jd = pd;
+ } else {
+ int ld = nodes[nodes[jd].next].left;
+ popLeft(nodes[jd].next);
+ pushRight(jd, ld);
+ if (less(ld, nodes[jd].left)) {
+ nodes[jd].item = nodes[ld].item;
+ nodes[jd].prio = nodes[jd].prio;
+ }
+ if (nodes[nodes[jd].next].item == nodes[ld].item) {
+ setPrio(nodes[jd].next);
+ }
+ jd = nodes[jd].left;
+ }
+ }
+ } else {
+ jd = nodes[jd].left;
+ }
+ }
+ }
+
+ void repairRight(int id) {
+ int jd = ~(classes[id].parent);
+ while (nodes[jd].right != -1) {
+ int kd = nodes[jd].right;
+ if (nodes[jd].size == 1) {
+ if (nodes[jd].parent < 0) {
+ classes[id].parent = ~kd;
+ classes[id].depth -= 1;
+ nodes[kd].parent = ~id;
+ deleteNode(jd);
+ jd = kd;
+ } else {
+ int pd = nodes[jd].parent;
+ if (nodes[nodes[jd].prev].size < cmax) {
+ pushRight(nodes[jd].prev, nodes[jd].right);
+ if (less(nodes[jd].right, nodes[jd].prev)) {
+ nodes[nodes[jd].prev].prio = nodes[nodes[jd].right].prio;
+ nodes[nodes[jd].prev].item = nodes[nodes[jd].right].item;
+ }
+ popRight(pd);
+ deleteNode(jd);
+ jd = pd;
+ } else {
+ int ld = nodes[nodes[jd].prev].right;
+ popRight(nodes[jd].prev);
+ pushLeft(jd, ld);
+ if (less(ld, nodes[jd].right)) {
+ nodes[jd].item = nodes[ld].item;
+ nodes[jd].prio = nodes[jd].prio;
+ }
+ if (nodes[nodes[jd].prev].item == nodes[ld].item) {
+ setPrio(nodes[jd].prev);
+ }
+ jd = nodes[jd].right;
+ }
+ }
+ } else {
+ jd = nodes[jd].right;
+ }
+ }
+ }
+
+
+ bool less(int id, int jd) const {
+ return comp(nodes[id].prio, nodes[jd].prio);
+ }
+
+ bool equal(int id, int jd) const {
+ return !less(id, jd) && !less(jd, id);
+ }
+
+
+ public:
+
+ /// \brief Returns true when the given class is alive.
+ ///
+ /// Returns true when the given class is alive, ie. the class is
+ /// not nested into other class.
+ bool alive(int cls) const {
+ return classes[cls].parent < 0;
+ }
+
+ /// \brief Returns true when the given class is trivial.
+ ///
+ /// Returns true when the given class is trivial, ie. the class
+ /// contains just one item directly.
+ bool trivial(int cls) const {
+ return classes[cls].left == -1;
+ }
+
+ /// \brief Constructs the union-find.
+ ///
+ /// Constructs the union-find.
+ /// \brief _index The index map of the union-find. The data
+ /// structure uses internally for store references.
+ HeapUnionFind(ItemIntMap& _index)
+ : index(_index), first_class(-1),
+ first_free_class(-1), first_free_node(-1) {}
+
+ /// \brief Insert a new node into a new component.
+ ///
+ /// Insert a new node into a new component.
+ /// \param item The item of the new node.
+ /// \param prio The priority of the new node.
+ /// \return The class id of the one-item-heap.
+ int insert(const Item& item, const Value& prio) {
+ int id = newNode();
+ nodes[id].item = item;
+ nodes[id].prio = prio;
+ nodes[id].size = 0;
+
+ nodes[id].prev = -1;
+ nodes[id].next = -1;
+
+ nodes[id].left = -1;
+ nodes[id].right = -1;
+
+ nodes[id].item = item;
+ index[item] = id;
+
+ int class_id = newClass();
+ classes[class_id].parent = ~id;
+ classes[class_id].depth = 0;
+
+ classes[class_id].left = -1;
+ classes[class_id].right = -1;
+
+ if (first_class != -1) {
+ classes[first_class].prev = class_id;
+ }
+ classes[class_id].next = first_class;
+ classes[class_id].prev = -1;
+ first_class = class_id;
+
+ nodes[id].parent = ~class_id;
+
+ return class_id;
+ }
+
+ /// \brief The class of the item.
+ ///
+ /// \return The alive class id of the item, which is not nested into
+ /// other classes.
+ ///
+ /// The time complexity is O(log(n)).
+ int find(const Item& item) const {
+ return findClass(index[item]);
+ }
+
+ /// \brief Joins the classes.
+ ///
+ /// The current function joins the given classes. The parameter is
+ /// an STL range which should be contains valid class ids. The
+ /// time complexity is O(log(n)*k) where n is the overall number
+ /// of the joined nodes and k is the number of classes.
+ /// \return The class of the joined classes.
+ /// \pre The range should contain at least two class ids.
+ template <typename Iterator>
+ int join(Iterator begin, Iterator end) {
+ std::vector<int> cs;
+ for (Iterator it = begin; it != end; ++it) {
+ cs.push_back(*it);
+ }
+
+ int class_id = newClass();
+ { // creation union-find
+
+ if (first_class != -1) {
+ classes[first_class].prev = class_id;
+ }
+ classes[class_id].next = first_class;
+ classes[class_id].prev = -1;
+ first_class = class_id;
+
+ classes[class_id].depth = classes[cs[0]].depth;
+ classes[class_id].parent = classes[cs[0]].parent;
+ nodes[~(classes[class_id].parent)].parent = ~class_id;
+
+ int l = cs[0];
+
+ classes[class_id].left = l;
+ classes[class_id].right = l;
+
+ if (classes[l].next != -1) {
+ classes[classes[l].next].prev = classes[l].prev;
+ }
+ classes[classes[l].prev].next = classes[l].next;
+
+ classes[l].prev = -1;
+ classes[l].next = -1;
+
+ classes[l].depth = leftNode(l);
+ classes[l].parent = class_id;
+
+ }
+
+ { // merging of heap
+ int l = class_id;
+ for (int ci = 1; ci < int(cs.size()); ++ci) {
+ int r = cs[ci];
+ int rln = leftNode(r);
+ if (classes[l].depth > classes[r].depth) {
+ int id = ~(classes[l].parent);
+ for (int i = classes[r].depth + 1; i < classes[l].depth; ++i) {
+ id = nodes[id].right;
+ }
+ while (id >= 0 && nodes[id].size == cmax) {
+ int new_id = newNode();
+ int right_id = nodes[id].right;
+
+ popRight(id);
+ if (nodes[id].item == nodes[right_id].item) {
+ setPrio(id);
+ }
+ push(new_id, right_id);
+ pushRight(new_id, ~(classes[r].parent));
+ setPrio(new_id);
+
+ id = nodes[id].parent;
+ classes[r].parent = ~new_id;
+ }
+ if (id < 0) {
+ int new_parent = newNode();
+ nodes[new_parent].next = -1;
+ nodes[new_parent].prev = -1;
+ nodes[new_parent].parent = ~l;
+
+ push(new_parent, ~(classes[l].parent));
+ pushRight(new_parent, ~(classes[r].parent));
+ setPrio(new_parent);
+
+ classes[l].parent = ~new_parent;
+ classes[l].depth += 1;
+ } else {
+ pushRight(id, ~(classes[r].parent));
+ while (id >= 0 && less(~(classes[r].parent), id)) {
+ nodes[id].prio = nodes[~(classes[r].parent)].prio;
+ nodes[id].item = nodes[~(classes[r].parent)].item;
+ id = nodes[id].parent;
+ }
+ }
+ } else if (classes[r].depth > classes[l].depth) {
+ int id = ~(classes[r].parent);
+ for (int i = classes[l].depth + 1; i < classes[r].depth; ++i) {
+ id = nodes[id].left;
+ }
+ while (id >= 0 && nodes[id].size == cmax) {
+ int new_id = newNode();
+ int left_id = nodes[id].left;
+
+ popLeft(id);
+ if (nodes[id].prio == nodes[left_id].prio) {
+ setPrio(id);
+ }
+ push(new_id, left_id);
+ pushLeft(new_id, ~(classes[l].parent));
+ setPrio(new_id);
+
+ id = nodes[id].parent;
+ classes[l].parent = ~new_id;
+
+ }
+ if (id < 0) {
+ int new_parent = newNode();
+ nodes[new_parent].next = -1;
+ nodes[new_parent].prev = -1;
+ nodes[new_parent].parent = ~l;
+
+ push(new_parent, ~(classes[r].parent));
+ pushLeft(new_parent, ~(classes[l].parent));
+ setPrio(new_parent);
+
+ classes[r].parent = ~new_parent;
+ classes[r].depth += 1;
+ } else {
+ pushLeft(id, ~(classes[l].parent));
+ while (id >= 0 && less(~(classes[l].parent), id)) {
+ nodes[id].prio = nodes[~(classes[l].parent)].prio;
+ nodes[id].item = nodes[~(classes[l].parent)].item;
+ id = nodes[id].parent;
+ }
+ }
+ nodes[~(classes[r].parent)].parent = ~l;
+ classes[l].parent = classes[r].parent;
+ classes[l].depth = classes[r].depth;
+ } else {
+ if (classes[l].depth != 0 &&
+ nodes[~(classes[l].parent)].size +
+ nodes[~(classes[r].parent)].size <= cmax) {
+ splice(~(classes[l].parent), ~(classes[r].parent));
+ deleteNode(~(classes[r].parent));
+ if (less(~(classes[r].parent), ~(classes[l].parent))) {
+ nodes[~(classes[l].parent)].prio =
+ nodes[~(classes[r].parent)].prio;
+ nodes[~(classes[l].parent)].item =
+ nodes[~(classes[r].parent)].item;
+ }
+ } else {
+ int new_parent = newNode();
+ nodes[new_parent].next = nodes[new_parent].prev = -1;
+ push(new_parent, ~(classes[l].parent));
+ pushRight(new_parent, ~(classes[r].parent));
+ setPrio(new_parent);
+
+ classes[l].parent = ~new_parent;
+ classes[l].depth += 1;
+ nodes[new_parent].parent = ~l;
+ }
+ }
+ if (classes[r].next != -1) {
+ classes[classes[r].next].prev = classes[r].prev;
+ }
+ classes[classes[r].prev].next = classes[r].next;
+
+ classes[r].prev = classes[l].right;
+ classes[classes[l].right].next = r;
+ classes[l].right = r;
+ classes[r].parent = l;
+
+ classes[r].next = -1;
+ classes[r].depth = rln;
+ }
+ }
+ return class_id;
+ }
+
+ /// \brief Split the class to subclasses.
+ ///
+ /// The current function splits the given class. The join, which
+ /// made the current class, stored a reference to the
+ /// subclasses. The \c splitClass() member restores the classes
+ /// and creates the heaps. The parameter is an STL output iterator
+ /// which will be filled with the subclass ids. The time
+ /// complexity is O(log(n)*k) where n is the overall number of
+ /// nodes in the splitted classes and k is the number of the
+ /// classes.
+ template <typename Iterator>
+ void split(int cls, Iterator out) {
+ std::vector<int> cs;
+ { // splitting union-find
+ int id = cls;
+ int l = classes[id].left;
+
+ classes[l].parent = classes[id].parent;
+ classes[l].depth = classes[id].depth;
+
+ nodes[~(classes[l].parent)].parent = ~l;
+
+ *out++ = l;
+
+ while (l != -1) {
+ cs.push_back(l);
+ l = classes[l].next;
+ }
+
+ classes[classes[id].right].next = first_class;
+ classes[first_class].prev = classes[id].right;
+ first_class = classes[id].left;
+
+ if (classes[id].next != -1) {
+ classes[classes[id].next].prev = classes[id].prev;
+ }
+ classes[classes[id].prev].next = classes[id].next;
+
+ deleteClass(id);
+ }
+
+ {
+ for (int i = 1; i < int(cs.size()); ++i) {
+ int l = classes[cs[i]].depth;
+ while (nodes[nodes[l].parent].left == l) {
+ l = nodes[l].parent;
+ }
+ int r = l;
+ while (nodes[l].parent >= 0) {
+ l = nodes[l].parent;
+ int new_node = newNode();
+
+ nodes[new_node].prev = -1;
+ nodes[new_node].next = -1;
+
+ split(r, new_node);
+ pushAfter(l, new_node);
+ setPrio(l);
+ setPrio(new_node);
+ r = new_node;
+ }
+ classes[cs[i]].parent = ~r;
+ classes[cs[i]].depth = classes[~(nodes[l].parent)].depth;
+ nodes[r].parent = ~cs[i];
+
+ nodes[l].next = -1;
+ nodes[r].prev = -1;
+
+ repairRight(~(nodes[l].parent));
+ repairLeft(cs[i]);
+
+ *out++ = cs[i];
+ }
+ }
+ }
+
+ /// \brief Gives back the priority of the current item.
+ ///
+ /// \return Gives back the priority of the current item.
+ const Value& operator[](const Item& item) const {
+ return nodes[index[item]].prio;
+ }
+
+ /// \brief Sets the priority of the current item.
+ ///
+ /// Sets the priority of the current item.
+ void set(const Item& item, const Value& prio) {
+ if (comp(prio, nodes[index[item]].prio)) {
+ decrease(item, prio);
+ } else if (!comp(prio, nodes[index[item]].prio)) {
+ increase(item, prio);
+ }
+ }
+
+ /// \brief Increase the priority of the current item.
+ ///
+ /// Increase the priority of the current item.
+ void increase(const Item& item, const Value& prio) {
+ int id = index[item];
+ int kd = nodes[id].parent;
+ nodes[id].prio = prio;
+ while (kd >= 0 && nodes[kd].item == item) {
+ setPrio(kd);
+ kd = nodes[kd].parent;
+ }
+ }
+
+ /// \brief Increase the priority of the current item.
+ ///
+ /// Increase the priority of the current item.
+ void decrease(const Item& item, const Value& prio) {
+ int id = index[item];
+ int kd = nodes[id].parent;
+ nodes[id].prio = prio;
+ while (kd >= 0 && less(id, kd)) {
+ nodes[kd].prio = prio;
+ nodes[kd].item = item;
+ kd = nodes[kd].parent;
+ }
+ }
+
+ /// \brief Gives back the minimum priority of the class.
+ ///
+ /// \return Gives back the minimum priority of the class.
+ const Value& classPrio(int cls) const {
+ return nodes[~(classes[cls].parent)].prio;
+ }
+
+ /// \brief Gives back the minimum priority item of the class.
+ ///
+ /// \return Gives back the minimum priority item of the class.
+ const Item& classTop(int cls) const {
+ return nodes[~(classes[cls].parent)].item;
+ }
+
+ /// \brief Gives back a representant item of the class.
+ ///
+ /// The representant is indpendent from the priorities of the
+ /// items.
+ /// \return Gives back a representant item of the class.
+ const Item& classRep(int id) const {
+ int parent = classes[id].parent;
+ return nodes[parent >= 0 ? classes[id].depth : leftNode(id)].item;
+ }
+
+ /// \brief Lemon style iterator for the items of a class.
+ ///
+ /// ClassIt is a lemon style iterator for the components. It iterates
+ /// on the items of a class. By example if you want to iterate on
+ /// each items of each classes then you may write the next code.
+ ///\code
+ /// for (ClassIt cit(huf); cit != INVALID; ++cit) {
+ /// std::cout << "Class: ";
+ /// for (ItemIt iit(huf, cit); iit != INVALID; ++iit) {
+ /// std::cout << toString(iit) << ' ' << std::endl;
+ /// }
+ /// std::cout << std::endl;
+ /// }
+ ///\endcode
+ class ItemIt {
+ private:
+
+ const HeapUnionFind* _huf;
+ int _id, _lid;
+
+ public:
+
+ /// \brief Default constructor
+ ///
+ /// Default constructor
+ ItemIt() {}
+
+ ItemIt(const HeapUnionFind& huf, int cls) : _huf(&huf) {
+ int id = cls;
+ int parent = _huf->classes[id].parent;
+ if (parent >= 0) {
+ _id = _huf->classes[id].depth;
+ if (_huf->classes[id].next != -1) {
+ _lid = _huf->classes[_huf->classes[id].next].depth;
+ } else {
+ _lid = -1;
+ }
+ } else {
+ _id = _huf->leftNode(id);
+ _lid = -1;
+ }
+ }
+
+ /// \brief Increment operator
+ ///
+ /// It steps to the next item in the class.
+ ItemIt& operator++() {
+ _id = _huf->nextNode(_id);
+ return *this;
+ }
+
+ /// \brief Conversion operator
+ ///
+ /// It converts the iterator to the current item.
+ operator const Item&() const {
+ return _huf->nodes[_id].item;
+ }
+
+ /// \brief Equality operator
+ ///
+ /// Equality operator
+ bool operator==(const ItemIt& i) {
+ return i._id == _id;
+ }
+
+ /// \brief Inequality operator
+ ///
+ /// Inequality operator
+ bool operator!=(const ItemIt& i) {
+ return i._id != _id;
+ }
+
+ /// \brief Equality operator
+ ///
+ /// Equality operator
+ bool operator==(Invalid) {
+ return _id == _lid;
+ }
+
+ /// \brief Inequality operator
+ ///
+ /// Inequality operator
+ bool operator!=(Invalid) {
+ return _id != _lid;
+ }
+
+ };
+
+ /// \brief Class iterator
+ ///
+ /// The iterator stores
+ class ClassIt {
+ private:
+
+ const HeapUnionFind* _huf;
+ int _id;
+
+ public:
+
+ ClassIt(const HeapUnionFind& huf)
+ : _huf(&huf), _id(huf.first_class) {}
+
+ ClassIt(const HeapUnionFind& huf, int cls)
+ : _huf(&huf), _id(huf.classes[cls].left) {}
+
+ ClassIt(Invalid) : _huf(0), _id(-1) {}
+
+ const ClassIt& operator++() {
+ _id = _huf->classes[_id].next;
+ return *this;
+ }
+
+ /// \brief Equality operator
+ ///
+ /// Equality operator
+ bool operator==(const ClassIt& i) {
+ return i._id == _id;
+ }
+
+ /// \brief Inequality operator
+ ///
+ /// Inequality operator
+ bool operator!=(const ClassIt& i) {
+ return i._id != _id;
+ }
+
+ operator int() const {
+ return _id;
+ }
+
+ };
+
+ };
+
//! @}
} //namespace lemon
More information about the Lemon-commits
mailing list