1.1 --- a/doc/groups.dox Thu Dec 27 13:40:16 2007 +0000
1.2 +++ b/doc/groups.dox Fri Dec 28 11:00:51 2007 +0000
1.3 @@ -318,8 +318,37 @@
1.4 \brief This group describes the algorithms
1.5 for find matchings in graphs and bipartite graphs.
1.6
1.7 -This group provides some algorithm objects and function
1.8 -to calculate matchings in graphs and bipartite graphs.
1.9 +This group provides some algorithm objects and function to calculate
1.10 +matchings in graphs and bipartite graphs. The general matching problem is
1.11 +finding a subset of the edges which does not shares common endpoints.
1.12 +
1.13 +There are several different algorithms for calculate matchings in
1.14 +graphs. The matching problems in bipartite graphs are generally
1.15 +easier than in general graphs. The goal of the matching optimization
1.16 +can be the finding maximum cardinality, maximum weight or minimum cost
1.17 +matching. The search can be constrained to find perfect or
1.18 +maximum cardinality matching.
1.19 +
1.20 +Lemon contains the next algorithms:
1.21 +- \ref lemon::MaxBipartiteMatching "MaxBipartiteMatching" Hopcroft-Karp
1.22 + augmenting path algorithm for calculate maximum cardinality matching in
1.23 + bipartite graphs
1.24 +- \ref lemon::PrBipartiteMatching "PrBipartiteMatching" Push-Relabel
1.25 + algorithm for calculate maximum cardinality matching in bipartite graphs
1.26 +- \ref lemon::MaxWeightedBipartiteMatching "MaxWeightedBipartiteMatching"
1.27 + Successive shortest path algorithm for calculate maximum weighted matching
1.28 + and maximum weighted bipartite matching in bipartite graph
1.29 +- \ref lemon::MinCostMaxBipartiteMatching "MinCostMaxBipartiteMatching"
1.30 + Successive shortest path algorithm for calculate minimum cost maximum
1.31 + matching in bipartite graph
1.32 +- \ref lemon::MaxMatching "MaxMatching" Edmond's blossom shrinking algorithm
1.33 + for calculate maximum cardinality matching in general graph
1.34 +- \ref lemon::MaxWeightedMatching "MaxWeightedMatching" Edmond's blossom
1.35 + shrinking algorithm for calculate maximum weighted matching in general
1.36 + graph
1.37 +- \ref lemon::MaxWeightedPerfectMatching "MaxWeightedPerfectMatching"
1.38 + Edmond's blossom shrinking algorithm for calculate maximum weighted
1.39 + perfect matching in general graph
1.40
1.41 \image html bipartite_matching.png
1.42 \image latex bipartite_matching.eps "Bipartite Matching" width=\textwidth
2.1 --- a/lemon/bin_heap.h Thu Dec 27 13:40:16 2007 +0000
2.2 +++ b/lemon/bin_heap.h Fri Dec 28 11:00:51 2007 +0000
2.3 @@ -52,10 +52,15 @@
2.4 class BinHeap {
2.5
2.6 public:
2.7 + ///\e
2.8 typedef _ItemIntMap ItemIntMap;
2.9 + ///\e
2.10 typedef _Prio Prio;
2.11 + ///\e
2.12 typedef typename ItemIntMap::Key Item;
2.13 + ///\e
2.14 typedef std::pair<Item,Prio> Pair;
2.15 + ///\e
2.16 typedef _Compare Compare;
2.17
2.18 /// \brief Type to represent the items states.
2.19 @@ -321,6 +326,19 @@
2.20 }
2.21 }
2.22
2.23 + /// \brief Replaces an item in the heap.
2.24 + ///
2.25 + /// The \c i item is replaced with \c j item. The \c i item should
2.26 + /// be in the heap, while the \c j should be out of the heap. The
2.27 + /// \c i item will out of the heap and \c j will be in the heap
2.28 + /// with the same prioriority as prevoiusly the \c i item.
2.29 + void replace(const Item& i, const Item& j) {
2.30 + int idx = iim[i];
2.31 + iim.set(i, iim[j]);
2.32 + iim.set(j, idx);
2.33 + data[idx].first = j;
2.34 + }
2.35 +
2.36 }; // class BinHeap
2.37
2.38 } // namespace lemon
3.1 --- a/lemon/max_matching.h Thu Dec 27 13:40:16 2007 +0000
3.2 +++ b/lemon/max_matching.h Fri Dec 28 11:00:51 2007 +0000
3.3 @@ -19,10 +19,14 @@
3.4 #ifndef LEMON_MAX_MATCHING_H
3.5 #define LEMON_MAX_MATCHING_H
3.6
3.7 +#include <vector>
3.8 #include <queue>
3.9 +#include <set>
3.10 +
3.11 #include <lemon/bits/invalid.h>
3.12 #include <lemon/unionfind.h>
3.13 #include <lemon/graph_utils.h>
3.14 +#include <lemon/bin_heap.h>
3.15
3.16 ///\ingroup matching
3.17 ///\file
3.18 @@ -31,7 +35,7 @@
3.19 namespace lemon {
3.20
3.21 ///\ingroup matching
3.22 -
3.23 + ///
3.24 ///\brief Edmonds' alternating forest maximum matching algorithm.
3.25 ///
3.26 ///This class provides Edmonds' alternating forest matching
3.27 @@ -618,6 +622,2500 @@
3.28 }
3.29
3.30 };
3.31 +
3.32 + /// \ingroup matching
3.33 + ///
3.34 + /// \brief Weighted matching in general undirected graphs
3.35 + ///
3.36 + /// This class provides an efficient implementation of Edmond's
3.37 + /// maximum weighted matching algorithm. The implementation is based
3.38 + /// on extensive use of priority queues and provides
3.39 + /// \f$O(nm\log(n))\f$ time complexity.
3.40 + ///
3.41 + /// The maximum weighted matching problem is to find undirected
3.42 + /// edges in the graph with maximum overall weight and no two of
3.43 + /// them shares their endpoints. The problem can be formulated with
3.44 + /// the next linear program:
3.45 + /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
3.46 + ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
3.47 + /// \f[x_e \ge 0\quad \forall e\in E\f]
3.48 + /// \f[\max \sum_{e\in E}x_ew_e\f]
3.49 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
3.50 + /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both endpoints in
3.51 + /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
3.52 + /// the nodes.
3.53 + ///
3.54 + /// The algorithm calculates an optimal matching and a proof of the
3.55 + /// optimality. The solution of the dual problem can be used to check
3.56 + /// the result of the algorithm. The dual linear problem is the next:
3.57 + /// \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]
3.58 + /// \f[y_u \ge 0 \quad \forall u \in V\f]
3.59 + /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
3.60 + /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
3.61 + ///
3.62 + /// The algorithm can be executed with \c run() or the \c init() and
3.63 + /// then the \c start() member functions. After it the matching can
3.64 + /// be asked with \c matching() or mate() functions. The dual
3.65 + /// solution can be get with \c nodeValue(), \c blossomNum() and \c
3.66 + /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
3.67 + /// "BlossomIt" nested class which is able to iterate on the nodes
3.68 + /// of a blossom. If the value type is integral then the dual
3.69 + /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
3.70 + ///
3.71 + /// \author Balazs Dezso
3.72 + template <typename _UGraph,
3.73 + typename _WeightMap = typename _UGraph::template UEdgeMap<int> >
3.74 + class MaxWeightedMatching {
3.75 + public:
3.76 +
3.77 + typedef _UGraph UGraph;
3.78 + typedef _WeightMap WeightMap;
3.79 + typedef typename WeightMap::Value Value;
3.80 +
3.81 + /// \brief Scaling factor for dual solution
3.82 + ///
3.83 + /// Scaling factor for dual solution, it is equal to 4 or 1
3.84 + /// according to the value type.
3.85 + static const int dualScale =
3.86 + std::numeric_limits<Value>::is_integer ? 4 : 1;
3.87 +
3.88 + typedef typename UGraph::template NodeMap<typename UGraph::Edge>
3.89 + MatchingMap;
3.90 +
3.91 + private:
3.92 +
3.93 + UGRAPH_TYPEDEFS(typename UGraph);
3.94 +
3.95 + typedef typename UGraph::template NodeMap<Value> NodePotential;
3.96 + typedef std::vector<Node> BlossomNodeList;
3.97 +
3.98 + struct BlossomVariable {
3.99 + int begin, end;
3.100 + Value value;
3.101 +
3.102 + BlossomVariable(int _begin, int _end, Value _value)
3.103 + : begin(_begin), end(_end), value(_value) {}
3.104 +
3.105 + };
3.106 +
3.107 + typedef std::vector<BlossomVariable> BlossomPotential;
3.108 +
3.109 + const UGraph& _ugraph;
3.110 + const WeightMap& _weight;
3.111 +
3.112 + MatchingMap* _matching;
3.113 +
3.114 + NodePotential* _node_potential;
3.115 +
3.116 + BlossomPotential _blossom_potential;
3.117 + BlossomNodeList _blossom_node_list;
3.118 +
3.119 + int _node_num;
3.120 + int _blossom_num;
3.121 +
3.122 + typedef typename UGraph::template NodeMap<int> NodeIntMap;
3.123 + typedef typename UGraph::template EdgeMap<int> EdgeIntMap;
3.124 + typedef typename UGraph::template UEdgeMap<int> UEdgeIntMap;
3.125 + typedef IntegerMap<int> IntIntMap;
3.126 +
3.127 + enum Status {
3.128 + EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
3.129 + };
3.130 +
3.131 + typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
3.132 + struct BlossomData {
3.133 + int tree;
3.134 + Status status;
3.135 + Edge pred, next;
3.136 + Value pot, offset;
3.137 + Node base;
3.138 + };
3.139 +
3.140 + NodeIntMap *_blossom_index;
3.141 + BlossomSet *_blossom_set;
3.142 + IntegerMap<BlossomData>* _blossom_data;
3.143 +
3.144 + NodeIntMap *_node_index;
3.145 + EdgeIntMap *_node_heap_index;
3.146 +
3.147 + struct NodeData {
3.148 +
3.149 + NodeData(EdgeIntMap& node_heap_index)
3.150 + : heap(node_heap_index) {}
3.151 +
3.152 + int blossom;
3.153 + Value pot;
3.154 + BinHeap<Value, EdgeIntMap> heap;
3.155 + std::map<int, Edge> heap_index;
3.156 +
3.157 + int tree;
3.158 + };
3.159 +
3.160 + IntegerMap<NodeData>* _node_data;
3.161 +
3.162 + typedef ExtendFindEnum<IntIntMap> TreeSet;
3.163 +
3.164 + IntIntMap *_tree_set_index;
3.165 + TreeSet *_tree_set;
3.166 +
3.167 + NodeIntMap *_delta1_index;
3.168 + BinHeap<Value, NodeIntMap> *_delta1;
3.169 +
3.170 + IntIntMap *_delta2_index;
3.171 + BinHeap<Value, IntIntMap> *_delta2;
3.172 +
3.173 + UEdgeIntMap *_delta3_index;
3.174 + BinHeap<Value, UEdgeIntMap> *_delta3;
3.175 +
3.176 + IntIntMap *_delta4_index;
3.177 + BinHeap<Value, IntIntMap> *_delta4;
3.178 +
3.179 + Value _delta_sum;
3.180 +
3.181 + void createStructures() {
3.182 + _node_num = countNodes(_ugraph);
3.183 + _blossom_num = _node_num * 3 / 2;
3.184 +
3.185 + if (!_matching) {
3.186 + _matching = new MatchingMap(_ugraph);
3.187 + }
3.188 + if (!_node_potential) {
3.189 + _node_potential = new NodePotential(_ugraph);
3.190 + }
3.191 + if (!_blossom_set) {
3.192 + _blossom_index = new NodeIntMap(_ugraph);
3.193 + _blossom_set = new BlossomSet(*_blossom_index);
3.194 + _blossom_data = new IntegerMap<BlossomData>(_blossom_num);
3.195 + }
3.196 +
3.197 + if (!_node_index) {
3.198 + _node_index = new NodeIntMap(_ugraph);
3.199 + _node_heap_index = new EdgeIntMap(_ugraph);
3.200 + _node_data = new IntegerMap<NodeData>(_node_num,
3.201 + NodeData(*_node_heap_index));
3.202 + }
3.203 +
3.204 + if (!_tree_set) {
3.205 + _tree_set_index = new IntIntMap(_blossom_num);
3.206 + _tree_set = new TreeSet(*_tree_set_index);
3.207 + }
3.208 + if (!_delta1) {
3.209 + _delta1_index = new NodeIntMap(_ugraph);
3.210 + _delta1 = new BinHeap<Value, NodeIntMap>(*_delta1_index);
3.211 + }
3.212 + if (!_delta2) {
3.213 + _delta2_index = new IntIntMap(_blossom_num);
3.214 + _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
3.215 + }
3.216 + if (!_delta3) {
3.217 + _delta3_index = new UEdgeIntMap(_ugraph);
3.218 + _delta3 = new BinHeap<Value, UEdgeIntMap>(*_delta3_index);
3.219 + }
3.220 + if (!_delta4) {
3.221 + _delta4_index = new IntIntMap(_blossom_num);
3.222 + _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
3.223 + }
3.224 + }
3.225 +
3.226 + void destroyStructures() {
3.227 + _node_num = countNodes(_ugraph);
3.228 + _blossom_num = _node_num * 3 / 2;
3.229 +
3.230 + if (_matching) {
3.231 + delete _matching;
3.232 + }
3.233 + if (_node_potential) {
3.234 + delete _node_potential;
3.235 + }
3.236 + if (_blossom_set) {
3.237 + delete _blossom_index;
3.238 + delete _blossom_set;
3.239 + delete _blossom_data;
3.240 + }
3.241 +
3.242 + if (_node_index) {
3.243 + delete _node_index;
3.244 + delete _node_heap_index;
3.245 + delete _node_data;
3.246 + }
3.247 +
3.248 + if (_tree_set) {
3.249 + delete _tree_set_index;
3.250 + delete _tree_set;
3.251 + }
3.252 + if (_delta1) {
3.253 + delete _delta1_index;
3.254 + delete _delta1;
3.255 + }
3.256 + if (_delta2) {
3.257 + delete _delta2_index;
3.258 + delete _delta2;
3.259 + }
3.260 + if (_delta3) {
3.261 + delete _delta3_index;
3.262 + delete _delta3;
3.263 + }
3.264 + if (_delta4) {
3.265 + delete _delta4_index;
3.266 + delete _delta4;
3.267 + }
3.268 + }
3.269 +
3.270 + void matchedToEven(int blossom, int tree) {
3.271 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
3.272 + _delta2->erase(blossom);
3.273 + }
3.274 +
3.275 + if (!_blossom_set->trivial(blossom)) {
3.276 + (*_blossom_data)[blossom].pot -=
3.277 + 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
3.278 + }
3.279 +
3.280 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
3.281 + n != INVALID; ++n) {
3.282 +
3.283 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
3.284 + int ni = (*_node_index)[n];
3.285 +
3.286 + (*_node_data)[ni].heap.clear();
3.287 + (*_node_data)[ni].heap_index.clear();
3.288 +
3.289 + (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
3.290 +
3.291 + _delta1->push(n, (*_node_data)[ni].pot);
3.292 +
3.293 + for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
3.294 + Node v = _ugraph.source(e);
3.295 + int vb = _blossom_set->find(v);
3.296 + int vi = (*_node_index)[v];
3.297 +
3.298 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
3.299 + dualScale * _weight[e];
3.300 +
3.301 + if ((*_blossom_data)[vb].status == EVEN) {
3.302 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
3.303 + _delta3->push(e, rw / 2);
3.304 + }
3.305 + } else if ((*_blossom_data)[vb].status == UNMATCHED) {
3.306 + if (_delta3->state(e) != _delta3->IN_HEAP) {
3.307 + _delta3->push(e, rw);
3.308 + }
3.309 + } else {
3.310 + typename std::map<int, Edge>::iterator it =
3.311 + (*_node_data)[vi].heap_index.find(tree);
3.312 +
3.313 + if (it != (*_node_data)[vi].heap_index.end()) {
3.314 + if ((*_node_data)[vi].heap[it->second] > rw) {
3.315 + (*_node_data)[vi].heap.replace(it->second, e);
3.316 + (*_node_data)[vi].heap.decrease(e, rw);
3.317 + it->second = e;
3.318 + }
3.319 + } else {
3.320 + (*_node_data)[vi].heap.push(e, rw);
3.321 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
3.322 + }
3.323 +
3.324 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
3.325 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
3.326 +
3.327 + if ((*_blossom_data)[vb].status == MATCHED) {
3.328 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
3.329 + _delta2->push(vb, _blossom_set->classPrio(vb) -
3.330 + (*_blossom_data)[vb].offset);
3.331 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
3.332 + (*_blossom_data)[vb].offset){
3.333 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
3.334 + (*_blossom_data)[vb].offset);
3.335 + }
3.336 + }
3.337 + }
3.338 + }
3.339 + }
3.340 + }
3.341 + (*_blossom_data)[blossom].offset = 0;
3.342 + }
3.343 +
3.344 + void matchedToOdd(int blossom) {
3.345 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
3.346 + _delta2->erase(blossom);
3.347 + }
3.348 + (*_blossom_data)[blossom].offset += _delta_sum;
3.349 + if (!_blossom_set->trivial(blossom)) {
3.350 + _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
3.351 + (*_blossom_data)[blossom].offset);
3.352 + }
3.353 + }
3.354 +
3.355 + void evenToMatched(int blossom, int tree) {
3.356 + if (!_blossom_set->trivial(blossom)) {
3.357 + (*_blossom_data)[blossom].pot += 2 * _delta_sum;
3.358 + }
3.359 +
3.360 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
3.361 + n != INVALID; ++n) {
3.362 + int ni = (*_node_index)[n];
3.363 + (*_node_data)[ni].pot -= _delta_sum;
3.364 +
3.365 + _delta1->erase(n);
3.366 +
3.367 + for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
3.368 + Node v = _ugraph.source(e);
3.369 + int vb = _blossom_set->find(v);
3.370 + int vi = (*_node_index)[v];
3.371 +
3.372 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
3.373 + dualScale * _weight[e];
3.374 +
3.375 + if (vb == blossom) {
3.376 + if (_delta3->state(e) == _delta3->IN_HEAP) {
3.377 + _delta3->erase(e);
3.378 + }
3.379 + } else if ((*_blossom_data)[vb].status == EVEN) {
3.380 +
3.381 + if (_delta3->state(e) == _delta3->IN_HEAP) {
3.382 + _delta3->erase(e);
3.383 + }
3.384 +
3.385 + int vt = _tree_set->find(vb);
3.386 +
3.387 + if (vt != tree) {
3.388 +
3.389 + Edge r = _ugraph.oppositeEdge(e);
3.390 +
3.391 + typename std::map<int, Edge>::iterator it =
3.392 + (*_node_data)[ni].heap_index.find(vt);
3.393 +
3.394 + if (it != (*_node_data)[ni].heap_index.end()) {
3.395 + if ((*_node_data)[ni].heap[it->second] > rw) {
3.396 + (*_node_data)[ni].heap.replace(it->second, r);
3.397 + (*_node_data)[ni].heap.decrease(r, rw);
3.398 + it->second = r;
3.399 + }
3.400 + } else {
3.401 + (*_node_data)[ni].heap.push(r, rw);
3.402 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
3.403 + }
3.404 +
3.405 + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
3.406 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
3.407 +
3.408 + if (_delta2->state(blossom) != _delta2->IN_HEAP) {
3.409 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
3.410 + (*_blossom_data)[blossom].offset);
3.411 + } else if ((*_delta2)[blossom] >
3.412 + _blossom_set->classPrio(blossom) -
3.413 + (*_blossom_data)[blossom].offset){
3.414 + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
3.415 + (*_blossom_data)[blossom].offset);
3.416 + }
3.417 + }
3.418 + }
3.419 +
3.420 + } else if ((*_blossom_data)[vb].status == UNMATCHED) {
3.421 + if (_delta3->state(e) == _delta3->IN_HEAP) {
3.422 + _delta3->erase(e);
3.423 + }
3.424 + } else {
3.425 +
3.426 + typename std::map<int, Edge>::iterator it =
3.427 + (*_node_data)[vi].heap_index.find(tree);
3.428 +
3.429 + if (it != (*_node_data)[vi].heap_index.end()) {
3.430 + (*_node_data)[vi].heap.erase(it->second);
3.431 + (*_node_data)[vi].heap_index.erase(it);
3.432 + if ((*_node_data)[vi].heap.empty()) {
3.433 + _blossom_set->increase(v, std::numeric_limits<Value>::max());
3.434 + } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
3.435 + _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
3.436 + }
3.437 +
3.438 + if ((*_blossom_data)[vb].status == MATCHED) {
3.439 + if (_blossom_set->classPrio(vb) ==
3.440 + std::numeric_limits<Value>::max()) {
3.441 + _delta2->erase(vb);
3.442 + } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
3.443 + (*_blossom_data)[vb].offset) {
3.444 + _delta2->increase(vb, _blossom_set->classPrio(vb) -
3.445 + (*_blossom_data)[vb].offset);
3.446 + }
3.447 + }
3.448 + }
3.449 + }
3.450 + }
3.451 + }
3.452 + }
3.453 +
3.454 + void oddToMatched(int blossom) {
3.455 + (*_blossom_data)[blossom].offset -= _delta_sum;
3.456 +
3.457 + if (_blossom_set->classPrio(blossom) !=
3.458 + std::numeric_limits<Value>::max()) {
3.459 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
3.460 + (*_blossom_data)[blossom].offset);
3.461 + }
3.462 +
3.463 + if (!_blossom_set->trivial(blossom)) {
3.464 + _delta4->erase(blossom);
3.465 + }
3.466 + }
3.467 +
3.468 + void oddToEven(int blossom, int tree) {
3.469 + if (!_blossom_set->trivial(blossom)) {
3.470 + _delta4->erase(blossom);
3.471 + (*_blossom_data)[blossom].pot -=
3.472 + 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
3.473 + }
3.474 +
3.475 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
3.476 + n != INVALID; ++n) {
3.477 + int ni = (*_node_index)[n];
3.478 +
3.479 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
3.480 +
3.481 + (*_node_data)[ni].heap.clear();
3.482 + (*_node_data)[ni].heap_index.clear();
3.483 + (*_node_data)[ni].pot +=
3.484 + 2 * _delta_sum - (*_blossom_data)[blossom].offset;
3.485 +
3.486 + _delta1->push(n, (*_node_data)[ni].pot);
3.487 +
3.488 + for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
3.489 + Node v = _ugraph.source(e);
3.490 + int vb = _blossom_set->find(v);
3.491 + int vi = (*_node_index)[v];
3.492 +
3.493 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
3.494 + dualScale * _weight[e];
3.495 +
3.496 + if ((*_blossom_data)[vb].status == EVEN) {
3.497 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
3.498 + _delta3->push(e, rw / 2);
3.499 + }
3.500 + } else if ((*_blossom_data)[vb].status == UNMATCHED) {
3.501 + if (_delta3->state(e) != _delta3->IN_HEAP) {
3.502 + _delta3->push(e, rw);
3.503 + }
3.504 + } else {
3.505 +
3.506 + typename std::map<int, Edge>::iterator it =
3.507 + (*_node_data)[vi].heap_index.find(tree);
3.508 +
3.509 + if (it != (*_node_data)[vi].heap_index.end()) {
3.510 + if ((*_node_data)[vi].heap[it->second] > rw) {
3.511 + (*_node_data)[vi].heap.replace(it->second, e);
3.512 + (*_node_data)[vi].heap.decrease(e, rw);
3.513 + it->second = e;
3.514 + }
3.515 + } else {
3.516 + (*_node_data)[vi].heap.push(e, rw);
3.517 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
3.518 + }
3.519 +
3.520 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
3.521 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
3.522 +
3.523 + if ((*_blossom_data)[vb].status == MATCHED) {
3.524 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
3.525 + _delta2->push(vb, _blossom_set->classPrio(vb) -
3.526 + (*_blossom_data)[vb].offset);
3.527 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
3.528 + (*_blossom_data)[vb].offset) {
3.529 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
3.530 + (*_blossom_data)[vb].offset);
3.531 + }
3.532 + }
3.533 + }
3.534 + }
3.535 + }
3.536 + }
3.537 + (*_blossom_data)[blossom].offset = 0;
3.538 + }
3.539 +
3.540 +
3.541 + void matchedToUnmatched(int blossom) {
3.542 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
3.543 + _delta2->erase(blossom);
3.544 + }
3.545 +
3.546 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
3.547 + n != INVALID; ++n) {
3.548 + int ni = (*_node_index)[n];
3.549 +
3.550 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
3.551 +
3.552 + (*_node_data)[ni].heap.clear();
3.553 + (*_node_data)[ni].heap_index.clear();
3.554 +
3.555 + for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
3.556 + Node v = _ugraph.target(e);
3.557 + int vb = _blossom_set->find(v);
3.558 + int vi = (*_node_index)[v];
3.559 +
3.560 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
3.561 + dualScale * _weight[e];
3.562 +
3.563 + if ((*_blossom_data)[vb].status == EVEN) {
3.564 + if (_delta3->state(e) != _delta3->IN_HEAP) {
3.565 + _delta3->push(e, rw);
3.566 + }
3.567 + }
3.568 + }
3.569 + }
3.570 + }
3.571 +
3.572 + void unmatchedToMatched(int blossom) {
3.573 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
3.574 + n != INVALID; ++n) {
3.575 + int ni = (*_node_index)[n];
3.576 +
3.577 + for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
3.578 + Node v = _ugraph.source(e);
3.579 + int vb = _blossom_set->find(v);
3.580 + int vi = (*_node_index)[v];
3.581 +
3.582 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
3.583 + dualScale * _weight[e];
3.584 +
3.585 + if (vb == blossom) {
3.586 + if (_delta3->state(e) == _delta3->IN_HEAP) {
3.587 + _delta3->erase(e);
3.588 + }
3.589 + } else if ((*_blossom_data)[vb].status == EVEN) {
3.590 +
3.591 + if (_delta3->state(e) == _delta3->IN_HEAP) {
3.592 + _delta3->erase(e);
3.593 + }
3.594 +
3.595 + int vt = _tree_set->find(vb);
3.596 +
3.597 + Edge r = _ugraph.oppositeEdge(e);
3.598 +
3.599 + typename std::map<int, Edge>::iterator it =
3.600 + (*_node_data)[ni].heap_index.find(vt);
3.601 +
3.602 + if (it != (*_node_data)[ni].heap_index.end()) {
3.603 + if ((*_node_data)[ni].heap[it->second] > rw) {
3.604 + (*_node_data)[ni].heap.replace(it->second, r);
3.605 + (*_node_data)[ni].heap.decrease(r, rw);
3.606 + it->second = r;
3.607 + }
3.608 + } else {
3.609 + (*_node_data)[ni].heap.push(r, rw);
3.610 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
3.611 + }
3.612 +
3.613 + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
3.614 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
3.615 +
3.616 + if (_delta2->state(blossom) != _delta2->IN_HEAP) {
3.617 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
3.618 + (*_blossom_data)[blossom].offset);
3.619 + } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
3.620 + (*_blossom_data)[blossom].offset){
3.621 + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
3.622 + (*_blossom_data)[blossom].offset);
3.623 + }
3.624 + }
3.625 +
3.626 + } else if ((*_blossom_data)[vb].status == UNMATCHED) {
3.627 + if (_delta3->state(e) == _delta3->IN_HEAP) {
3.628 + _delta3->erase(e);
3.629 + }
3.630 + }
3.631 + }
3.632 + }
3.633 + }
3.634 +
3.635 + void alternatePath(int even, int tree) {
3.636 + int odd;
3.637 +
3.638 + evenToMatched(even, tree);
3.639 + (*_blossom_data)[even].status = MATCHED;
3.640 +
3.641 + while ((*_blossom_data)[even].pred != INVALID) {
3.642 + odd = _blossom_set->find(_ugraph.target((*_blossom_data)[even].pred));
3.643 + (*_blossom_data)[odd].status = MATCHED;
3.644 + oddToMatched(odd);
3.645 + (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
3.646 +
3.647 + even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].pred));
3.648 + (*_blossom_data)[even].status = MATCHED;
3.649 + evenToMatched(even, tree);
3.650 + (*_blossom_data)[even].next =
3.651 + _ugraph.oppositeEdge((*_blossom_data)[odd].pred);
3.652 + }
3.653 +
3.654 + }
3.655 +
3.656 + void destroyTree(int tree) {
3.657 + for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
3.658 + if ((*_blossom_data)[b].status == EVEN) {
3.659 + (*_blossom_data)[b].status = MATCHED;
3.660 + evenToMatched(b, tree);
3.661 + } else if ((*_blossom_data)[b].status == ODD) {
3.662 + (*_blossom_data)[b].status = MATCHED;
3.663 + oddToMatched(b);
3.664 + }
3.665 + }
3.666 + _tree_set->eraseClass(tree);
3.667 + }
3.668 +
3.669 +
3.670 + void unmatchNode(const Node& node) {
3.671 + int blossom = _blossom_set->find(node);
3.672 + int tree = _tree_set->find(blossom);
3.673 +
3.674 + alternatePath(blossom, tree);
3.675 + destroyTree(tree);
3.676 +
3.677 + (*_blossom_data)[blossom].status = UNMATCHED;
3.678 + (*_blossom_data)[blossom].base = node;
3.679 + matchedToUnmatched(blossom);
3.680 + }
3.681 +
3.682 +
3.683 + void augmentOnEdge(const UEdge& edge) {
3.684 +
3.685 + int left = _blossom_set->find(_ugraph.source(edge));
3.686 + int right = _blossom_set->find(_ugraph.target(edge));
3.687 +
3.688 + if ((*_blossom_data)[left].status == EVEN) {
3.689 + int left_tree = _tree_set->find(left);
3.690 + alternatePath(left, left_tree);
3.691 + destroyTree(left_tree);
3.692 + } else {
3.693 + (*_blossom_data)[left].status = MATCHED;
3.694 + unmatchedToMatched(left);
3.695 + }
3.696 +
3.697 + if ((*_blossom_data)[right].status == EVEN) {
3.698 + int right_tree = _tree_set->find(right);
3.699 + alternatePath(right, right_tree);
3.700 + destroyTree(right_tree);
3.701 + } else {
3.702 + (*_blossom_data)[right].status = MATCHED;
3.703 + unmatchedToMatched(right);
3.704 + }
3.705 +
3.706 + (*_blossom_data)[left].next = _ugraph.direct(edge, true);
3.707 + (*_blossom_data)[right].next = _ugraph.direct(edge, false);
3.708 + }
3.709 +
3.710 + void extendOnEdge(const Edge& edge) {
3.711 + int base = _blossom_set->find(_ugraph.target(edge));
3.712 + int tree = _tree_set->find(base);
3.713 +
3.714 + int odd = _blossom_set->find(_ugraph.source(edge));
3.715 + _tree_set->insert(odd, tree);
3.716 + (*_blossom_data)[odd].status = ODD;
3.717 + matchedToOdd(odd);
3.718 + (*_blossom_data)[odd].pred = edge;
3.719 +
3.720 + int even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].next));
3.721 + (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
3.722 + _tree_set->insert(even, tree);
3.723 + (*_blossom_data)[even].status = EVEN;
3.724 + matchedToEven(even, tree);
3.725 + }
3.726 +
3.727 + void shrinkOnEdge(const UEdge& uedge, int tree) {
3.728 + int nca = -1;
3.729 + std::vector<int> left_path, right_path;
3.730 +
3.731 + {
3.732 + std::set<int> left_set, right_set;
3.733 + int left = _blossom_set->find(_ugraph.source(uedge));
3.734 + left_path.push_back(left);
3.735 + left_set.insert(left);
3.736 +
3.737 + int right = _blossom_set->find(_ugraph.target(uedge));
3.738 + right_path.push_back(right);
3.739 + right_set.insert(right);
3.740 +
3.741 + while (true) {
3.742 +
3.743 + if ((*_blossom_data)[left].pred == INVALID) break;
3.744 +
3.745 + left =
3.746 + _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
3.747 + left_path.push_back(left);
3.748 + left =
3.749 + _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
3.750 + left_path.push_back(left);
3.751 +
3.752 + left_set.insert(left);
3.753 +
3.754 + if (right_set.find(left) != right_set.end()) {
3.755 + nca = left;
3.756 + break;
3.757 + }
3.758 +
3.759 + if ((*_blossom_data)[right].pred == INVALID) break;
3.760 +
3.761 + right =
3.762 + _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
3.763 + right_path.push_back(right);
3.764 + right =
3.765 + _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
3.766 + right_path.push_back(right);
3.767 +
3.768 + right_set.insert(right);
3.769 +
3.770 + if (left_set.find(right) != left_set.end()) {
3.771 + nca = right;
3.772 + break;
3.773 + }
3.774 +
3.775 + }
3.776 +
3.777 + if (nca == -1) {
3.778 + if ((*_blossom_data)[left].pred == INVALID) {
3.779 + nca = right;
3.780 + while (left_set.find(nca) == left_set.end()) {
3.781 + nca =
3.782 + _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
3.783 + right_path.push_back(nca);
3.784 + nca =
3.785 + _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
3.786 + right_path.push_back(nca);
3.787 + }
3.788 + } else {
3.789 + nca = left;
3.790 + while (right_set.find(nca) == right_set.end()) {
3.791 + nca =
3.792 + _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
3.793 + left_path.push_back(nca);
3.794 + nca =
3.795 + _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
3.796 + left_path.push_back(nca);
3.797 + }
3.798 + }
3.799 + }
3.800 + }
3.801 +
3.802 + std::vector<int> subblossoms;
3.803 + Edge prev;
3.804 +
3.805 + prev = _ugraph.direct(uedge, true);
3.806 + for (int i = 0; left_path[i] != nca; i += 2) {
3.807 + subblossoms.push_back(left_path[i]);
3.808 + (*_blossom_data)[left_path[i]].next = prev;
3.809 + _tree_set->erase(left_path[i]);
3.810 +
3.811 + subblossoms.push_back(left_path[i + 1]);
3.812 + (*_blossom_data)[left_path[i + 1]].status = EVEN;
3.813 + oddToEven(left_path[i + 1], tree);
3.814 + _tree_set->erase(left_path[i + 1]);
3.815 + prev = _ugraph.oppositeEdge((*_blossom_data)[left_path[i + 1]].pred);
3.816 + }
3.817 +
3.818 + int k = 0;
3.819 + while (right_path[k] != nca) ++k;
3.820 +
3.821 + subblossoms.push_back(nca);
3.822 + (*_blossom_data)[nca].next = prev;
3.823 +
3.824 + for (int i = k - 2; i >= 0; i -= 2) {
3.825 + subblossoms.push_back(right_path[i + 1]);
3.826 + (*_blossom_data)[right_path[i + 1]].status = EVEN;
3.827 + oddToEven(right_path[i + 1], tree);
3.828 + _tree_set->erase(right_path[i + 1]);
3.829 +
3.830 + (*_blossom_data)[right_path[i + 1]].next =
3.831 + (*_blossom_data)[right_path[i + 1]].pred;
3.832 +
3.833 + subblossoms.push_back(right_path[i]);
3.834 + _tree_set->erase(right_path[i]);
3.835 + }
3.836 +
3.837 + int surface =
3.838 + _blossom_set->join(subblossoms.begin(), subblossoms.end());
3.839 +
3.840 + for (int i = 0; i < int(subblossoms.size()); ++i) {
3.841 + if (!_blossom_set->trivial(subblossoms[i])) {
3.842 + (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
3.843 + }
3.844 + (*_blossom_data)[subblossoms[i]].status = MATCHED;
3.845 + }
3.846 +
3.847 + (*_blossom_data)[surface].pot = -2 * _delta_sum;
3.848 + (*_blossom_data)[surface].offset = 0;
3.849 + (*_blossom_data)[surface].status = EVEN;
3.850 + (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
3.851 + (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
3.852 +
3.853 + _tree_set->insert(surface, tree);
3.854 + _tree_set->erase(nca);
3.855 + }
3.856 +
3.857 + void splitBlossom(int blossom) {
3.858 + Edge next = (*_blossom_data)[blossom].next;
3.859 + Edge pred = (*_blossom_data)[blossom].pred;
3.860 +
3.861 + int tree = _tree_set->find(blossom);
3.862 +
3.863 + (*_blossom_data)[blossom].status = MATCHED;
3.864 + oddToMatched(blossom);
3.865 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
3.866 + _delta2->erase(blossom);
3.867 + }
3.868 +
3.869 + std::vector<int> subblossoms;
3.870 + _blossom_set->split(blossom, std::back_inserter(subblossoms));
3.871 +
3.872 + Value offset = (*_blossom_data)[blossom].offset;
3.873 + int b = _blossom_set->find(_ugraph.source(pred));
3.874 + int d = _blossom_set->find(_ugraph.source(next));
3.875 +
3.876 + int ib, id;
3.877 + for (int i = 0; i < int(subblossoms.size()); ++i) {
3.878 + if (subblossoms[i] == b) ib = i;
3.879 + if (subblossoms[i] == d) id = i;
3.880 +
3.881 + (*_blossom_data)[subblossoms[i]].offset = offset;
3.882 + if (!_blossom_set->trivial(subblossoms[i])) {
3.883 + (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
3.884 + }
3.885 + if (_blossom_set->classPrio(subblossoms[i]) !=
3.886 + std::numeric_limits<Value>::max()) {
3.887 + _delta2->push(subblossoms[i],
3.888 + _blossom_set->classPrio(subblossoms[i]) -
3.889 + (*_blossom_data)[subblossoms[i]].offset);
3.890 + }
3.891 + }
3.892 +
3.893 + if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
3.894 + for (int i = (id + 1) % subblossoms.size();
3.895 + i != ib; i = (i + 2) % subblossoms.size()) {
3.896 + int sb = subblossoms[i];
3.897 + int tb = subblossoms[(i + 1) % subblossoms.size()];
3.898 + (*_blossom_data)[sb].next =
3.899 + _ugraph.oppositeEdge((*_blossom_data)[tb].next);
3.900 + }
3.901 +
3.902 + for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
3.903 + int sb = subblossoms[i];
3.904 + int tb = subblossoms[(i + 1) % subblossoms.size()];
3.905 + int ub = subblossoms[(i + 2) % subblossoms.size()];
3.906 +
3.907 + (*_blossom_data)[sb].status = ODD;
3.908 + matchedToOdd(sb);
3.909 + _tree_set->insert(sb, tree);
3.910 + (*_blossom_data)[sb].pred = pred;
3.911 + (*_blossom_data)[sb].next =
3.912 + _ugraph.oppositeEdge((*_blossom_data)[tb].next);
3.913 +
3.914 + pred = (*_blossom_data)[ub].next;
3.915 +
3.916 + (*_blossom_data)[tb].status = EVEN;
3.917 + matchedToEven(tb, tree);
3.918 + _tree_set->insert(tb, tree);
3.919 + (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
3.920 + }
3.921 +
3.922 + (*_blossom_data)[subblossoms[id]].status = ODD;
3.923 + matchedToOdd(subblossoms[id]);
3.924 + _tree_set->insert(subblossoms[id], tree);
3.925 + (*_blossom_data)[subblossoms[id]].next = next;
3.926 + (*_blossom_data)[subblossoms[id]].pred = pred;
3.927 +
3.928 + } else {
3.929 +
3.930 + for (int i = (ib + 1) % subblossoms.size();
3.931 + i != id; i = (i + 2) % subblossoms.size()) {
3.932 + int sb = subblossoms[i];
3.933 + int tb = subblossoms[(i + 1) % subblossoms.size()];
3.934 + (*_blossom_data)[sb].next =
3.935 + _ugraph.oppositeEdge((*_blossom_data)[tb].next);
3.936 + }
3.937 +
3.938 + for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
3.939 + int sb = subblossoms[i];
3.940 + int tb = subblossoms[(i + 1) % subblossoms.size()];
3.941 + int ub = subblossoms[(i + 2) % subblossoms.size()];
3.942 +
3.943 + (*_blossom_data)[sb].status = ODD;
3.944 + matchedToOdd(sb);
3.945 + _tree_set->insert(sb, tree);
3.946 + (*_blossom_data)[sb].next = next;
3.947 + (*_blossom_data)[sb].pred =
3.948 + _ugraph.oppositeEdge((*_blossom_data)[tb].next);
3.949 +
3.950 + (*_blossom_data)[tb].status = EVEN;
3.951 + matchedToEven(tb, tree);
3.952 + _tree_set->insert(tb, tree);
3.953 + (*_blossom_data)[tb].pred =
3.954 + (*_blossom_data)[tb].next =
3.955 + _ugraph.oppositeEdge((*_blossom_data)[ub].next);
3.956 + next = (*_blossom_data)[ub].next;
3.957 + }
3.958 +
3.959 + (*_blossom_data)[subblossoms[ib]].status = ODD;
3.960 + matchedToOdd(subblossoms[ib]);
3.961 + _tree_set->insert(subblossoms[ib], tree);
3.962 + (*_blossom_data)[subblossoms[ib]].next = next;
3.963 + (*_blossom_data)[subblossoms[ib]].pred = pred;
3.964 + }
3.965 + _tree_set->erase(blossom);
3.966 + }
3.967 +
3.968 + void extractBlossom(int blossom, const Node& base, const Edge& matching) {
3.969 + if (_blossom_set->trivial(blossom)) {
3.970 + int bi = (*_node_index)[base];
3.971 + Value pot = (*_node_data)[bi].pot;
3.972 +
3.973 + _matching->set(base, matching);
3.974 + _blossom_node_list.push_back(base);
3.975 + _node_potential->set(base, pot);
3.976 + } else {
3.977 +
3.978 + Value pot = (*_blossom_data)[blossom].pot;
3.979 + int bn = _blossom_node_list.size();
3.980 +
3.981 + std::vector<int> subblossoms;
3.982 + _blossom_set->split(blossom, std::back_inserter(subblossoms));
3.983 + int b = _blossom_set->find(base);
3.984 + int ib = -1;
3.985 + for (int i = 0; i < int(subblossoms.size()); ++i) {
3.986 + if (subblossoms[i] == b) { ib = i; break; }
3.987 + }
3.988 +
3.989 + for (int i = 1; i < int(subblossoms.size()); i += 2) {
3.990 + int sb = subblossoms[(ib + i) % subblossoms.size()];
3.991 + int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
3.992 +
3.993 + Edge m = (*_blossom_data)[tb].next;
3.994 + extractBlossom(sb, _ugraph.target(m), _ugraph.oppositeEdge(m));
3.995 + extractBlossom(tb, _ugraph.source(m), m);
3.996 + }
3.997 + extractBlossom(subblossoms[ib], base, matching);
3.998 +
3.999 + int en = _blossom_node_list.size();
3.1000 +
3.1001 + _blossom_potential.push_back(BlossomVariable(bn, en, pot));
3.1002 + }
3.1003 + }
3.1004 +
3.1005 + void extractMatching() {
3.1006 + std::vector<int> blossoms;
3.1007 + for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
3.1008 + blossoms.push_back(c);
3.1009 + }
3.1010 +
3.1011 + for (int i = 0; i < int(blossoms.size()); ++i) {
3.1012 + if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
3.1013 +
3.1014 + Value offset = (*_blossom_data)[blossoms[i]].offset;
3.1015 + (*_blossom_data)[blossoms[i]].pot += 2 * offset;
3.1016 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
3.1017 + n != INVALID; ++n) {
3.1018 + (*_node_data)[(*_node_index)[n]].pot -= offset;
3.1019 + }
3.1020 +
3.1021 + Edge matching = (*_blossom_data)[blossoms[i]].next;
3.1022 + Node base = _ugraph.source(matching);
3.1023 + extractBlossom(blossoms[i], base, matching);
3.1024 + } else {
3.1025 + Node base = (*_blossom_data)[blossoms[i]].base;
3.1026 + extractBlossom(blossoms[i], base, INVALID);
3.1027 + }
3.1028 + }
3.1029 + }
3.1030 +
3.1031 + public:
3.1032 +
3.1033 + /// \brief Constructor
3.1034 + ///
3.1035 + /// Constructor.
3.1036 + MaxWeightedMatching(const UGraph& ugraph, const WeightMap& weight)
3.1037 + : _ugraph(ugraph), _weight(weight), _matching(0),
3.1038 + _node_potential(0), _blossom_potential(), _blossom_node_list(),
3.1039 + _node_num(0), _blossom_num(0),
3.1040 +
3.1041 + _blossom_index(0), _blossom_set(0), _blossom_data(0),
3.1042 + _node_index(0), _node_heap_index(0), _node_data(0),
3.1043 + _tree_set_index(0), _tree_set(0),
3.1044 +
3.1045 + _delta1_index(0), _delta1(0),
3.1046 + _delta2_index(0), _delta2(0),
3.1047 + _delta3_index(0), _delta3(0),
3.1048 + _delta4_index(0), _delta4(0),
3.1049 +
3.1050 + _delta_sum() {}
3.1051 +
3.1052 + ~MaxWeightedMatching() {
3.1053 + destroyStructures();
3.1054 + }
3.1055 +
3.1056 + /// \name Execution control
3.1057 + /// The simplest way to execute the algorithm is to use the member
3.1058 + /// \c run() member function.
3.1059 +
3.1060 + ///@{
3.1061 +
3.1062 + /// \brief Initialize the algorithm
3.1063 + ///
3.1064 + /// Initialize the algorithm
3.1065 + void init() {
3.1066 + createStructures();
3.1067 +
3.1068 + for (EdgeIt e(_ugraph); e != INVALID; ++e) {
3.1069 + _node_heap_index->set(e, BinHeap<Value, EdgeIntMap>::PRE_HEAP);
3.1070 + }
3.1071 + for (NodeIt n(_ugraph); n != INVALID; ++n) {
3.1072 + _delta1_index->set(n, _delta1->PRE_HEAP);
3.1073 + }
3.1074 + for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
3.1075 + _delta3_index->set(e, _delta3->PRE_HEAP);
3.1076 + }
3.1077 + for (int i = 0; i < _blossom_num; ++i) {
3.1078 + _delta2_index->set(i, _delta2->PRE_HEAP);
3.1079 + _delta4_index->set(i, _delta4->PRE_HEAP);
3.1080 + }
3.1081 +
3.1082 + int index = 0;
3.1083 + for (NodeIt n(_ugraph); n != INVALID; ++n) {
3.1084 + Value max = 0;
3.1085 + for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
3.1086 + if (_ugraph.target(e) == n) continue;
3.1087 + if ((dualScale * _weight[e]) / 2 > max) {
3.1088 + max = (dualScale * _weight[e]) / 2;
3.1089 + }
3.1090 + }
3.1091 + _node_index->set(n, index);
3.1092 + (*_node_data)[index].pot = max;
3.1093 + _delta1->push(n, max);
3.1094 + int blossom =
3.1095 + _blossom_set->insert(n, std::numeric_limits<Value>::max());
3.1096 +
3.1097 + _tree_set->insert(blossom);
3.1098 +
3.1099 + (*_blossom_data)[blossom].status = EVEN;
3.1100 + (*_blossom_data)[blossom].pred = INVALID;
3.1101 + (*_blossom_data)[blossom].next = INVALID;
3.1102 + (*_blossom_data)[blossom].pot = 0;
3.1103 + (*_blossom_data)[blossom].offset = 0;
3.1104 + ++index;
3.1105 + }
3.1106 + for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
3.1107 + int si = (*_node_index)[_ugraph.source(e)];
3.1108 + int ti = (*_node_index)[_ugraph.target(e)];
3.1109 + if (_ugraph.source(e) != _ugraph.target(e)) {
3.1110 + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
3.1111 + dualScale * _weight[e]) / 2);
3.1112 + }
3.1113 + }
3.1114 + }
3.1115 +
3.1116 + /// \brief Starts the algorithm
3.1117 + ///
3.1118 + /// Starts the algorithm
3.1119 + void start() {
3.1120 + enum OpType {
3.1121 + D1, D2, D3, D4
3.1122 + };
3.1123 +
3.1124 + int unmatched = _node_num;
3.1125 + while (unmatched > 0) {
3.1126 + Value d1 = !_delta1->empty() ?
3.1127 + _delta1->prio() : std::numeric_limits<Value>::max();
3.1128 +
3.1129 + Value d2 = !_delta2->empty() ?
3.1130 + _delta2->prio() : std::numeric_limits<Value>::max();
3.1131 +
3.1132 + Value d3 = !_delta3->empty() ?
3.1133 + _delta3->prio() : std::numeric_limits<Value>::max();
3.1134 +
3.1135 + Value d4 = !_delta4->empty() ?
3.1136 + _delta4->prio() : std::numeric_limits<Value>::max();
3.1137 +
3.1138 + _delta_sum = d1; OpType ot = D1;
3.1139 + if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
3.1140 + if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
3.1141 + if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
3.1142 +
3.1143 +
3.1144 + switch (ot) {
3.1145 + case D1:
3.1146 + {
3.1147 + Node n = _delta1->top();
3.1148 + unmatchNode(n);
3.1149 + --unmatched;
3.1150 + }
3.1151 + break;
3.1152 + case D2:
3.1153 + {
3.1154 + int blossom = _delta2->top();
3.1155 + Node n = _blossom_set->classTop(blossom);
3.1156 + Edge e = (*_node_data)[(*_node_index)[n]].heap.top();
3.1157 + extendOnEdge(e);
3.1158 + }
3.1159 + break;
3.1160 + case D3:
3.1161 + {
3.1162 + UEdge e = _delta3->top();
3.1163 +
3.1164 + int left_blossom = _blossom_set->find(_ugraph.source(e));
3.1165 + int right_blossom = _blossom_set->find(_ugraph.target(e));
3.1166 +
3.1167 + if (left_blossom == right_blossom) {
3.1168 + _delta3->pop();
3.1169 + } else {
3.1170 + int left_tree;
3.1171 + if ((*_blossom_data)[left_blossom].status == EVEN) {
3.1172 + left_tree = _tree_set->find(left_blossom);
3.1173 + } else {
3.1174 + left_tree = -1;
3.1175 + ++unmatched;
3.1176 + }
3.1177 + int right_tree;
3.1178 + if ((*_blossom_data)[right_blossom].status == EVEN) {
3.1179 + right_tree = _tree_set->find(right_blossom);
3.1180 + } else {
3.1181 + right_tree = -1;
3.1182 + ++unmatched;
3.1183 + }
3.1184 +
3.1185 + if (left_tree == right_tree) {
3.1186 + shrinkOnEdge(e, left_tree);
3.1187 + } else {
3.1188 + augmentOnEdge(e);
3.1189 + unmatched -= 2;
3.1190 + }
3.1191 + }
3.1192 + } break;
3.1193 + case D4:
3.1194 + splitBlossom(_delta4->top());
3.1195 + break;
3.1196 + }
3.1197 + }
3.1198 + extractMatching();
3.1199 + }
3.1200 +
3.1201 + /// \brief Runs %MaxWeightedMatching algorithm.
3.1202 + ///
3.1203 + /// This method runs the %MaxWeightedMatching algorithm.
3.1204 + ///
3.1205 + /// \note mwm.run() is just a shortcut of the following code.
3.1206 + /// \code
3.1207 + /// mwm.init();
3.1208 + /// mwm.start();
3.1209 + /// \endcode
3.1210 + void run() {
3.1211 + init();
3.1212 + start();
3.1213 + }
3.1214 +
3.1215 + /// @}
3.1216 +
3.1217 + /// \name Primal solution
3.1218 + /// Functions for get the primal solution, ie. the matching.
3.1219 +
3.1220 + /// @{
3.1221 +
3.1222 + /// \brief Returns the matching value.
3.1223 + ///
3.1224 + /// Returns the matching value.
3.1225 + Value matchingValue() const {
3.1226 + Value sum = 0;
3.1227 + for (NodeIt n(_ugraph); n != INVALID; ++n) {
3.1228 + if ((*_matching)[n] != INVALID) {
3.1229 + sum += _weight[(*_matching)[n]];
3.1230 + }
3.1231 + }
3.1232 + return sum /= 2;
3.1233 + }
3.1234 +
3.1235 + /// \brief Returns true when the edge is in the matching.
3.1236 + ///
3.1237 + /// Returns true when the edge is in the matching.
3.1238 + bool matching(const UEdge& edge) const {
3.1239 + return (*_matching)[_ugraph.source(edge)] == _ugraph.direct(edge, true);
3.1240 + }
3.1241 +
3.1242 + /// \brief Returns the incident matching edge.
3.1243 + ///
3.1244 + /// Returns the incident matching edge from given node. If the
3.1245 + /// node is not matched then it gives back \c INVALID.
3.1246 + Edge matching(const Node& node) const {
3.1247 + return (*_matching)[node];
3.1248 + }
3.1249 +
3.1250 + /// \brief Returns the mate of the node.
3.1251 + ///
3.1252 + /// Returns the adjancent node in a mathcing edge. If the node is
3.1253 + /// not matched then it gives back \c INVALID.
3.1254 + Node mate(const Node& node) const {
3.1255 + return (*_matching)[node] != INVALID ?
3.1256 + _ugraph.target((*_matching)[node]) : INVALID;
3.1257 + }
3.1258 +
3.1259 + /// @}
3.1260 +
3.1261 + /// \name Dual solution
3.1262 + /// Functions for get the dual solution.
3.1263 +
3.1264 + /// @{
3.1265 +
3.1266 + /// \brief Returns the value of the dual solution.
3.1267 + ///
3.1268 + /// Returns the value of the dual solution. It should be equal to
3.1269 + /// the primal value scaled by \ref dualScale "dual scale".
3.1270 + Value dualValue() const {
3.1271 + Value sum = 0;
3.1272 + for (NodeIt n(_ugraph); n != INVALID; ++n) {
3.1273 + sum += nodeValue(n);
3.1274 + }
3.1275 + for (int i = 0; i < blossomNum(); ++i) {
3.1276 + sum += blossomValue(i) * (blossomSize(i) / 2);
3.1277 + }
3.1278 + return sum;
3.1279 + }
3.1280 +
3.1281 + /// \brief Returns the value of the node.
3.1282 + ///
3.1283 + /// Returns the the value of the node.
3.1284 + Value nodeValue(const Node& n) const {
3.1285 + return (*_node_potential)[n];
3.1286 + }
3.1287 +
3.1288 + /// \brief Returns the number of the blossoms in the basis.
3.1289 + ///
3.1290 + /// Returns the number of the blossoms in the basis.
3.1291 + /// \see BlossomIt
3.1292 + int blossomNum() const {
3.1293 + return _blossom_potential.size();
3.1294 + }
3.1295 +
3.1296 +
3.1297 + /// \brief Returns the number of the nodes in the blossom.
3.1298 + ///
3.1299 + /// Returns the number of the nodes in the blossom.
3.1300 + int blossomSize(int k) const {
3.1301 + return _blossom_potential[k].end - _blossom_potential[k].begin;
3.1302 + }
3.1303 +
3.1304 + /// \brief Returns the value of the blossom.
3.1305 + ///
3.1306 + /// Returns the the value of the blossom.
3.1307 + /// \see BlossomIt
3.1308 + Value blossomValue(int k) const {
3.1309 + return _blossom_potential[k].value;
3.1310 + }
3.1311 +
3.1312 + /// \brief Lemon iterator for get the items of the blossom.
3.1313 + ///
3.1314 + /// Lemon iterator for get the nodes of the blossom. This class
3.1315 + /// provides a common style lemon iterator which gives back a
3.1316 + /// subset of the nodes.
3.1317 + class BlossomIt {
3.1318 + public:
3.1319 +
3.1320 + /// \brief Constructor.
3.1321 + ///
3.1322 + /// Constructor for get the nodes of the variable.
3.1323 + BlossomIt(const MaxWeightedMatching& algorithm, int variable)
3.1324 + : _algorithm(&algorithm)
3.1325 + {
3.1326 + _index = _algorithm->_blossom_potential[variable].begin;
3.1327 + _last = _algorithm->_blossom_potential[variable].end;
3.1328 + }
3.1329 +
3.1330 + /// \brief Invalid constructor.
3.1331 + ///
3.1332 + /// Invalid constructor.
3.1333 + BlossomIt(Invalid) : _index(-1) {}
3.1334 +
3.1335 + /// \brief Conversion to node.
3.1336 + ///
3.1337 + /// Conversion to node.
3.1338 + operator Node() const {
3.1339 + return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
3.1340 + }
3.1341 +
3.1342 + /// \brief Increment operator.
3.1343 + ///
3.1344 + /// Increment operator.
3.1345 + BlossomIt& operator++() {
3.1346 + ++_index;
3.1347 + if (_index == _last) {
3.1348 + _index = -1;
3.1349 + }
3.1350 + return *this;
3.1351 + }
3.1352 +
3.1353 + bool operator==(const BlossomIt& it) const {
3.1354 + return _index == it._index;
3.1355 + }
3.1356 + bool operator!=(const BlossomIt& it) const {
3.1357 + return _index != it._index;
3.1358 + }
3.1359 +
3.1360 + private:
3.1361 + const MaxWeightedMatching* _algorithm;
3.1362 + int _last;
3.1363 + int _index;
3.1364 + };
3.1365 +
3.1366 + /// @}
3.1367 +
3.1368 + };
3.1369 +
3.1370 + /// \ingroup matching
3.1371 + ///
3.1372 + /// \brief Weighted perfect matching in general undirected graphs
3.1373 + ///
3.1374 + /// This class provides an efficient implementation of Edmond's
3.1375 + /// maximum weighted perfecr matching algorithm. The implementation
3.1376 + /// is based on extensive use of priority queues and provides
3.1377 + /// \f$O(nm\log(n))\f$ time complexity.
3.1378 + ///
3.1379 + /// The maximum weighted matching problem is to find undirected
3.1380 + /// edges in the graph with maximum overall weight and no two of
3.1381 + /// them shares their endpoints and covers all nodes. The problem
3.1382 + /// can be formulated with the next linear program:
3.1383 + /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
3.1384 + ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
3.1385 + /// \f[x_e \ge 0\quad \forall e\in E\f]
3.1386 + /// \f[\max \sum_{e\in E}x_ew_e\f]
3.1387 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
3.1388 + /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both endpoints in
3.1389 + /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
3.1390 + /// the nodes.
3.1391 + ///
3.1392 + /// The algorithm calculates an optimal matching and a proof of the
3.1393 + /// optimality. The solution of the dual problem can be used to check
3.1394 + /// the result of the algorithm. The dual linear problem is the next:
3.1395 + /// \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]
3.1396 + /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
3.1397 + /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
3.1398 + ///
3.1399 + /// The algorithm can be executed with \c run() or the \c init() and
3.1400 + /// then the \c start() member functions. After it the matching can
3.1401 + /// be asked with \c matching() or mate() functions. The dual
3.1402 + /// solution can be get with \c nodeValue(), \c blossomNum() and \c
3.1403 + /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
3.1404 + /// "BlossomIt" nested class which is able to iterate on the nodes
3.1405 + /// of a blossom. If the value type is integral then the dual
3.1406 + /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
3.1407 + ///
3.1408 + /// \author Balazs Dezso
3.1409 + template <typename _UGraph,
3.1410 + typename _WeightMap = typename _UGraph::template UEdgeMap<int> >
3.1411 + class MaxWeightedPerfectMatching {
3.1412 + public:
3.1413 +
3.1414 + typedef _UGraph UGraph;
3.1415 + typedef _WeightMap WeightMap;
3.1416 + typedef typename WeightMap::Value Value;
3.1417 +
3.1418 + /// \brief Scaling factor for dual solution
3.1419 + ///
3.1420 + /// Scaling factor for dual solution, it is equal to 4 or 1
3.1421 + /// according to the value type.
3.1422 + static const int dualScale =
3.1423 + std::numeric_limits<Value>::is_integer ? 4 : 1;
3.1424 +
3.1425 + typedef typename UGraph::template NodeMap<typename UGraph::Edge>
3.1426 + MatchingMap;
3.1427 +
3.1428 + private:
3.1429 +
3.1430 + UGRAPH_TYPEDEFS(typename UGraph);
3.1431 +
3.1432 + typedef typename UGraph::template NodeMap<Value> NodePotential;
3.1433 + typedef std::vector<Node> BlossomNodeList;
3.1434 +
3.1435 + struct BlossomVariable {
3.1436 + int begin, end;
3.1437 + Value value;
3.1438 +
3.1439 + BlossomVariable(int _begin, int _end, Value _value)
3.1440 + : begin(_begin), end(_end), value(_value) {}
3.1441 +
3.1442 + };
3.1443 +
3.1444 + typedef std::vector<BlossomVariable> BlossomPotential;
3.1445 +
3.1446 + const UGraph& _ugraph;
3.1447 + const WeightMap& _weight;
3.1448 +
3.1449 + MatchingMap* _matching;
3.1450 +
3.1451 + NodePotential* _node_potential;
3.1452 +
3.1453 + BlossomPotential _blossom_potential;
3.1454 + BlossomNodeList _blossom_node_list;
3.1455 +
3.1456 + int _node_num;
3.1457 + int _blossom_num;
3.1458 +
3.1459 + typedef typename UGraph::template NodeMap<int> NodeIntMap;
3.1460 + typedef typename UGraph::template EdgeMap<int> EdgeIntMap;
3.1461 + typedef typename UGraph::template UEdgeMap<int> UEdgeIntMap;
3.1462 + typedef IntegerMap<int> IntIntMap;
3.1463 +
3.1464 + enum Status {
3.1465 + EVEN = -1, MATCHED = 0, ODD = 1
3.1466 + };
3.1467 +
3.1468 + typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
3.1469 + struct BlossomData {
3.1470 + int tree;
3.1471 + Status status;
3.1472 + Edge pred, next;
3.1473 + Value pot, offset;
3.1474 + };
3.1475 +
3.1476 + NodeIntMap *_blossom_index;
3.1477 + BlossomSet *_blossom_set;
3.1478 + IntegerMap<BlossomData>* _blossom_data;
3.1479 +
3.1480 + NodeIntMap *_node_index;
3.1481 + EdgeIntMap *_node_heap_index;
3.1482 +
3.1483 + struct NodeData {
3.1484 +
3.1485 + NodeData(EdgeIntMap& node_heap_index)
3.1486 + : heap(node_heap_index) {}
3.1487 +
3.1488 + int blossom;
3.1489 + Value pot;
3.1490 + BinHeap<Value, EdgeIntMap> heap;
3.1491 + std::map<int, Edge> heap_index;
3.1492 +
3.1493 + int tree;
3.1494 + };
3.1495 +
3.1496 + IntegerMap<NodeData>* _node_data;
3.1497 +
3.1498 + typedef ExtendFindEnum<IntIntMap> TreeSet;
3.1499 +
3.1500 + IntIntMap *_tree_set_index;
3.1501 + TreeSet *_tree_set;
3.1502 +
3.1503 + IntIntMap *_delta2_index;
3.1504 + BinHeap<Value, IntIntMap> *_delta2;
3.1505 +
3.1506 + UEdgeIntMap *_delta3_index;
3.1507 + BinHeap<Value, UEdgeIntMap> *_delta3;
3.1508 +
3.1509 + IntIntMap *_delta4_index;
3.1510 + BinHeap<Value, IntIntMap> *_delta4;
3.1511 +
3.1512 + Value _delta_sum;
3.1513 +
3.1514 + void createStructures() {
3.1515 + _node_num = countNodes(_ugraph);
3.1516 + _blossom_num = _node_num * 3 / 2;
3.1517 +
3.1518 + if (!_matching) {
3.1519 + _matching = new MatchingMap(_ugraph);
3.1520 + }
3.1521 + if (!_node_potential) {
3.1522 + _node_potential = new NodePotential(_ugraph);
3.1523 + }
3.1524 + if (!_blossom_set) {
3.1525 + _blossom_index = new NodeIntMap(_ugraph);
3.1526 + _blossom_set = new BlossomSet(*_blossom_index);
3.1527 + _blossom_data = new IntegerMap<BlossomData>(_blossom_num);
3.1528 + }
3.1529 +
3.1530 + if (!_node_index) {
3.1531 + _node_index = new NodeIntMap(_ugraph);
3.1532 + _node_heap_index = new EdgeIntMap(_ugraph);
3.1533 + _node_data = new IntegerMap<NodeData>(_node_num,
3.1534 + NodeData(*_node_heap_index));
3.1535 + }
3.1536 +
3.1537 + if (!_tree_set) {
3.1538 + _tree_set_index = new IntIntMap(_blossom_num);
3.1539 + _tree_set = new TreeSet(*_tree_set_index);
3.1540 + }
3.1541 + if (!_delta2) {
3.1542 + _delta2_index = new IntIntMap(_blossom_num);
3.1543 + _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
3.1544 + }
3.1545 + if (!_delta3) {
3.1546 + _delta3_index = new UEdgeIntMap(_ugraph);
3.1547 + _delta3 = new BinHeap<Value, UEdgeIntMap>(*_delta3_index);
3.1548 + }
3.1549 + if (!_delta4) {
3.1550 + _delta4_index = new IntIntMap(_blossom_num);
3.1551 + _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
3.1552 + }
3.1553 + }
3.1554 +
3.1555 + void destroyStructures() {
3.1556 + _node_num = countNodes(_ugraph);
3.1557 + _blossom_num = _node_num * 3 / 2;
3.1558 +
3.1559 + if (_matching) {
3.1560 + delete _matching;
3.1561 + }
3.1562 + if (_node_potential) {
3.1563 + delete _node_potential;
3.1564 + }
3.1565 + if (_blossom_set) {
3.1566 + delete _blossom_index;
3.1567 + delete _blossom_set;
3.1568 + delete _blossom_data;
3.1569 + }
3.1570 +
3.1571 + if (_node_index) {
3.1572 + delete _node_index;
3.1573 + delete _node_heap_index;
3.1574 + delete _node_data;
3.1575 + }
3.1576 +
3.1577 + if (_tree_set) {
3.1578 + delete _tree_set_index;
3.1579 + delete _tree_set;
3.1580 + }
3.1581 + if (_delta2) {
3.1582 + delete _delta2_index;
3.1583 + delete _delta2;
3.1584 + }
3.1585 + if (_delta3) {
3.1586 + delete _delta3_index;
3.1587 + delete _delta3;
3.1588 + }
3.1589 + if (_delta4) {
3.1590 + delete _delta4_index;
3.1591 + delete _delta4;
3.1592 + }
3.1593 + }
3.1594 +
3.1595 + void matchedToEven(int blossom, int tree) {
3.1596 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
3.1597 + _delta2->erase(blossom);
3.1598 + }
3.1599 +
3.1600 + if (!_blossom_set->trivial(blossom)) {
3.1601 + (*_blossom_data)[blossom].pot -=
3.1602 + 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
3.1603 + }
3.1604 +
3.1605 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
3.1606 + n != INVALID; ++n) {
3.1607 +
3.1608 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
3.1609 + int ni = (*_node_index)[n];
3.1610 +
3.1611 + (*_node_data)[ni].heap.clear();
3.1612 + (*_node_data)[ni].heap_index.clear();
3.1613 +
3.1614 + (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
3.1615 +
3.1616 + for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
3.1617 + Node v = _ugraph.source(e);
3.1618 + int vb = _blossom_set->find(v);
3.1619 + int vi = (*_node_index)[v];
3.1620 +
3.1621 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
3.1622 + dualScale * _weight[e];
3.1623 +
3.1624 + if ((*_blossom_data)[vb].status == EVEN) {
3.1625 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
3.1626 + _delta3->push(e, rw / 2);
3.1627 + }
3.1628 + } else {
3.1629 + typename std::map<int, Edge>::iterator it =
3.1630 + (*_node_data)[vi].heap_index.find(tree);
3.1631 +
3.1632 + if (it != (*_node_data)[vi].heap_index.end()) {
3.1633 + if ((*_node_data)[vi].heap[it->second] > rw) {
3.1634 + (*_node_data)[vi].heap.replace(it->second, e);
3.1635 + (*_node_data)[vi].heap.decrease(e, rw);
3.1636 + it->second = e;
3.1637 + }
3.1638 + } else {
3.1639 + (*_node_data)[vi].heap.push(e, rw);
3.1640 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
3.1641 + }
3.1642 +
3.1643 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
3.1644 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
3.1645 +
3.1646 + if ((*_blossom_data)[vb].status == MATCHED) {
3.1647 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
3.1648 + _delta2->push(vb, _blossom_set->classPrio(vb) -
3.1649 + (*_blossom_data)[vb].offset);
3.1650 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
3.1651 + (*_blossom_data)[vb].offset){
3.1652 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
3.1653 + (*_blossom_data)[vb].offset);
3.1654 + }
3.1655 + }
3.1656 + }
3.1657 + }
3.1658 + }
3.1659 + }
3.1660 + (*_blossom_data)[blossom].offset = 0;
3.1661 + }
3.1662 +
3.1663 + void matchedToOdd(int blossom) {
3.1664 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
3.1665 + _delta2->erase(blossom);
3.1666 + }
3.1667 + (*_blossom_data)[blossom].offset += _delta_sum;
3.1668 + if (!_blossom_set->trivial(blossom)) {
3.1669 + _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
3.1670 + (*_blossom_data)[blossom].offset);
3.1671 + }
3.1672 + }
3.1673 +
3.1674 + void evenToMatched(int blossom, int tree) {
3.1675 + if (!_blossom_set->trivial(blossom)) {
3.1676 + (*_blossom_data)[blossom].pot += 2 * _delta_sum;
3.1677 + }
3.1678 +
3.1679 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
3.1680 + n != INVALID; ++n) {
3.1681 + int ni = (*_node_index)[n];
3.1682 + (*_node_data)[ni].pot -= _delta_sum;
3.1683 +
3.1684 + for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
3.1685 + Node v = _ugraph.source(e);
3.1686 + int vb = _blossom_set->find(v);
3.1687 + int vi = (*_node_index)[v];
3.1688 +
3.1689 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
3.1690 + dualScale * _weight[e];
3.1691 +
3.1692 + if (vb == blossom) {
3.1693 + if (_delta3->state(e) == _delta3->IN_HEAP) {
3.1694 + _delta3->erase(e);
3.1695 + }
3.1696 + } else if ((*_blossom_data)[vb].status == EVEN) {
3.1697 +
3.1698 + if (_delta3->state(e) == _delta3->IN_HEAP) {
3.1699 + _delta3->erase(e);
3.1700 + }
3.1701 +
3.1702 + int vt = _tree_set->find(vb);
3.1703 +
3.1704 + if (vt != tree) {
3.1705 +
3.1706 + Edge r = _ugraph.oppositeEdge(e);
3.1707 +
3.1708 + typename std::map<int, Edge>::iterator it =
3.1709 + (*_node_data)[ni].heap_index.find(vt);
3.1710 +
3.1711 + if (it != (*_node_data)[ni].heap_index.end()) {
3.1712 + if ((*_node_data)[ni].heap[it->second] > rw) {
3.1713 + (*_node_data)[ni].heap.replace(it->second, r);
3.1714 + (*_node_data)[ni].heap.decrease(r, rw);
3.1715 + it->second = r;
3.1716 + }
3.1717 + } else {
3.1718 + (*_node_data)[ni].heap.push(r, rw);
3.1719 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
3.1720 + }
3.1721 +
3.1722 + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
3.1723 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
3.1724 +
3.1725 + if (_delta2->state(blossom) != _delta2->IN_HEAP) {
3.1726 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
3.1727 + (*_blossom_data)[blossom].offset);
3.1728 + } else if ((*_delta2)[blossom] >
3.1729 + _blossom_set->classPrio(blossom) -
3.1730 + (*_blossom_data)[blossom].offset){
3.1731 + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
3.1732 + (*_blossom_data)[blossom].offset);
3.1733 + }
3.1734 + }
3.1735 + }
3.1736 + } else {
3.1737 +
3.1738 + typename std::map<int, Edge>::iterator it =
3.1739 + (*_node_data)[vi].heap_index.find(tree);
3.1740 +
3.1741 + if (it != (*_node_data)[vi].heap_index.end()) {
3.1742 + (*_node_data)[vi].heap.erase(it->second);
3.1743 + (*_node_data)[vi].heap_index.erase(it);
3.1744 + if ((*_node_data)[vi].heap.empty()) {
3.1745 + _blossom_set->increase(v, std::numeric_limits<Value>::max());
3.1746 + } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
3.1747 + _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
3.1748 + }
3.1749 +
3.1750 + if ((*_blossom_data)[vb].status == MATCHED) {
3.1751 + if (_blossom_set->classPrio(vb) ==
3.1752 + std::numeric_limits<Value>::max()) {
3.1753 + _delta2->erase(vb);
3.1754 + } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
3.1755 + (*_blossom_data)[vb].offset) {
3.1756 + _delta2->increase(vb, _blossom_set->classPrio(vb) -
3.1757 + (*_blossom_data)[vb].offset);
3.1758 + }
3.1759 + }
3.1760 + }
3.1761 + }
3.1762 + }
3.1763 + }
3.1764 + }
3.1765 +
3.1766 + void oddToMatched(int blossom) {
3.1767 + (*_blossom_data)[blossom].offset -= _delta_sum;
3.1768 +
3.1769 + if (_blossom_set->classPrio(blossom) !=
3.1770 + std::numeric_limits<Value>::max()) {
3.1771 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
3.1772 + (*_blossom_data)[blossom].offset);
3.1773 + }
3.1774 +
3.1775 + if (!_blossom_set->trivial(blossom)) {
3.1776 + _delta4->erase(blossom);
3.1777 + }
3.1778 + }
3.1779 +
3.1780 + void oddToEven(int blossom, int tree) {
3.1781 + if (!_blossom_set->trivial(blossom)) {
3.1782 + _delta4->erase(blossom);
3.1783 + (*_blossom_data)[blossom].pot -=
3.1784 + 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
3.1785 + }
3.1786 +
3.1787 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
3.1788 + n != INVALID; ++n) {
3.1789 + int ni = (*_node_index)[n];
3.1790 +
3.1791 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
3.1792 +
3.1793 + (*_node_data)[ni].heap.clear();
3.1794 + (*_node_data)[ni].heap_index.clear();
3.1795 + (*_node_data)[ni].pot +=
3.1796 + 2 * _delta_sum - (*_blossom_data)[blossom].offset;
3.1797 +
3.1798 + for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
3.1799 + Node v = _ugraph.source(e);
3.1800 + int vb = _blossom_set->find(v);
3.1801 + int vi = (*_node_index)[v];
3.1802 +
3.1803 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
3.1804 + dualScale * _weight[e];
3.1805 +
3.1806 + if ((*_blossom_data)[vb].status == EVEN) {
3.1807 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
3.1808 + _delta3->push(e, rw / 2);
3.1809 + }
3.1810 + } else {
3.1811 +
3.1812 + typename std::map<int, Edge>::iterator it =
3.1813 + (*_node_data)[vi].heap_index.find(tree);
3.1814 +
3.1815 + if (it != (*_node_data)[vi].heap_index.end()) {
3.1816 + if ((*_node_data)[vi].heap[it->second] > rw) {
3.1817 + (*_node_data)[vi].heap.replace(it->second, e);
3.1818 + (*_node_data)[vi].heap.decrease(e, rw);
3.1819 + it->second = e;
3.1820 + }
3.1821 + } else {
3.1822 + (*_node_data)[vi].heap.push(e, rw);
3.1823 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
3.1824 + }
3.1825 +
3.1826 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
3.1827 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
3.1828 +
3.1829 + if ((*_blossom_data)[vb].status == MATCHED) {
3.1830 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
3.1831 + _delta2->push(vb, _blossom_set->classPrio(vb) -
3.1832 + (*_blossom_data)[vb].offset);
3.1833 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
3.1834 + (*_blossom_data)[vb].offset) {
3.1835 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
3.1836 + (*_blossom_data)[vb].offset);
3.1837 + }
3.1838 + }
3.1839 + }
3.1840 + }
3.1841 + }
3.1842 + }
3.1843 + (*_blossom_data)[blossom].offset = 0;
3.1844 + }
3.1845 +
3.1846 + void alternatePath(int even, int tree) {
3.1847 + int odd;
3.1848 +
3.1849 + evenToMatched(even, tree);
3.1850 + (*_blossom_data)[even].status = MATCHED;
3.1851 +
3.1852 + while ((*_blossom_data)[even].pred != INVALID) {
3.1853 + odd = _blossom_set->find(_ugraph.target((*_blossom_data)[even].pred));
3.1854 + (*_blossom_data)[odd].status = MATCHED;
3.1855 + oddToMatched(odd);
3.1856 + (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
3.1857 +
3.1858 + even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].pred));
3.1859 + (*_blossom_data)[even].status = MATCHED;
3.1860 + evenToMatched(even, tree);
3.1861 + (*_blossom_data)[even].next =
3.1862 + _ugraph.oppositeEdge((*_blossom_data)[odd].pred);
3.1863 + }
3.1864 +
3.1865 + }
3.1866 +
3.1867 + void destroyTree(int tree) {
3.1868 + for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
3.1869 + if ((*_blossom_data)[b].status == EVEN) {
3.1870 + (*_blossom_data)[b].status = MATCHED;
3.1871 + evenToMatched(b, tree);
3.1872 + } else if ((*_blossom_data)[b].status == ODD) {
3.1873 + (*_blossom_data)[b].status = MATCHED;
3.1874 + oddToMatched(b);
3.1875 + }
3.1876 + }
3.1877 + _tree_set->eraseClass(tree);
3.1878 + }
3.1879 +
3.1880 + void augmentOnEdge(const UEdge& edge) {
3.1881 +
3.1882 + int left = _blossom_set->find(_ugraph.source(edge));
3.1883 + int right = _blossom_set->find(_ugraph.target(edge));
3.1884 +
3.1885 + int left_tree = _tree_set->find(left);
3.1886 + alternatePath(left, left_tree);
3.1887 + destroyTree(left_tree);
3.1888 +
3.1889 + int right_tree = _tree_set->find(right);
3.1890 + alternatePath(right, right_tree);
3.1891 + destroyTree(right_tree);
3.1892 +
3.1893 + (*_blossom_data)[left].next = _ugraph.direct(edge, true);
3.1894 + (*_blossom_data)[right].next = _ugraph.direct(edge, false);
3.1895 + }
3.1896 +
3.1897 + void extendOnEdge(const Edge& edge) {
3.1898 + int base = _blossom_set->find(_ugraph.target(edge));
3.1899 + int tree = _tree_set->find(base);
3.1900 +
3.1901 + int odd = _blossom_set->find(_ugraph.source(edge));
3.1902 + _tree_set->insert(odd, tree);
3.1903 + (*_blossom_data)[odd].status = ODD;
3.1904 + matchedToOdd(odd);
3.1905 + (*_blossom_data)[odd].pred = edge;
3.1906 +
3.1907 + int even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].next));
3.1908 + (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
3.1909 + _tree_set->insert(even, tree);
3.1910 + (*_blossom_data)[even].status = EVEN;
3.1911 + matchedToEven(even, tree);
3.1912 + }
3.1913 +
3.1914 + void shrinkOnEdge(const UEdge& uedge, int tree) {
3.1915 + int nca = -1;
3.1916 + std::vector<int> left_path, right_path;
3.1917 +
3.1918 + {
3.1919 + std::set<int> left_set, right_set;
3.1920 + int left = _blossom_set->find(_ugraph.source(uedge));
3.1921 + left_path.push_back(left);
3.1922 + left_set.insert(left);
3.1923 +
3.1924 + int right = _blossom_set->find(_ugraph.target(uedge));
3.1925 + right_path.push_back(right);
3.1926 + right_set.insert(right);
3.1927 +
3.1928 + while (true) {
3.1929 +
3.1930 + if ((*_blossom_data)[left].pred == INVALID) break;
3.1931 +
3.1932 + left =
3.1933 + _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
3.1934 + left_path.push_back(left);
3.1935 + left =
3.1936 + _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
3.1937 + left_path.push_back(left);
3.1938 +
3.1939 + left_set.insert(left);
3.1940 +
3.1941 + if (right_set.find(left) != right_set.end()) {
3.1942 + nca = left;
3.1943 + break;
3.1944 + }
3.1945 +
3.1946 + if ((*_blossom_data)[right].pred == INVALID) break;
3.1947 +
3.1948 + right =
3.1949 + _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
3.1950 + right_path.push_back(right);
3.1951 + right =
3.1952 + _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
3.1953 + right_path.push_back(right);
3.1954 +
3.1955 + right_set.insert(right);
3.1956 +
3.1957 + if (left_set.find(right) != left_set.end()) {
3.1958 + nca = right;
3.1959 + break;
3.1960 + }
3.1961 +
3.1962 + }
3.1963 +
3.1964 + if (nca == -1) {
3.1965 + if ((*_blossom_data)[left].pred == INVALID) {
3.1966 + nca = right;
3.1967 + while (left_set.find(nca) == left_set.end()) {
3.1968 + nca =
3.1969 + _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
3.1970 + right_path.push_back(nca);
3.1971 + nca =
3.1972 + _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
3.1973 + right_path.push_back(nca);
3.1974 + }
3.1975 + } else {
3.1976 + nca = left;
3.1977 + while (right_set.find(nca) == right_set.end()) {
3.1978 + nca =
3.1979 + _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
3.1980 + left_path.push_back(nca);
3.1981 + nca =
3.1982 + _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
3.1983 + left_path.push_back(nca);
3.1984 + }
3.1985 + }
3.1986 + }
3.1987 + }
3.1988 +
3.1989 + std::vector<int> subblossoms;
3.1990 + Edge prev;
3.1991 +
3.1992 + prev = _ugraph.direct(uedge, true);
3.1993 + for (int i = 0; left_path[i] != nca; i += 2) {
3.1994 + subblossoms.push_back(left_path[i]);
3.1995 + (*_blossom_data)[left_path[i]].next = prev;
3.1996 + _tree_set->erase(left_path[i]);
3.1997 +
3.1998 + subblossoms.push_back(left_path[i + 1]);
3.1999 + (*_blossom_data)[left_path[i + 1]].status = EVEN;
3.2000 + oddToEven(left_path[i + 1], tree);
3.2001 + _tree_set->erase(left_path[i + 1]);
3.2002 + prev = _ugraph.oppositeEdge((*_blossom_data)[left_path[i + 1]].pred);
3.2003 + }
3.2004 +
3.2005 + int k = 0;
3.2006 + while (right_path[k] != nca) ++k;
3.2007 +
3.2008 + subblossoms.push_back(nca);
3.2009 + (*_blossom_data)[nca].next = prev;
3.2010 +
3.2011 + for (int i = k - 2; i >= 0; i -= 2) {
3.2012 + subblossoms.push_back(right_path[i + 1]);
3.2013 + (*_blossom_data)[right_path[i + 1]].status = EVEN;
3.2014 + oddToEven(right_path[i + 1], tree);
3.2015 + _tree_set->erase(right_path[i + 1]);
3.2016 +
3.2017 + (*_blossom_data)[right_path[i + 1]].next =
3.2018 + (*_blossom_data)[right_path[i + 1]].pred;
3.2019 +
3.2020 + subblossoms.push_back(right_path[i]);
3.2021 + _tree_set->erase(right_path[i]);
3.2022 + }
3.2023 +
3.2024 + int surface =
3.2025 + _blossom_set->join(subblossoms.begin(), subblossoms.end());
3.2026 +
3.2027 + for (int i = 0; i < int(subblossoms.size()); ++i) {
3.2028 + if (!_blossom_set->trivial(subblossoms[i])) {
3.2029 + (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
3.2030 + }
3.2031 + (*_blossom_data)[subblossoms[i]].status = MATCHED;
3.2032 + }
3.2033 +
3.2034 + (*_blossom_data)[surface].pot = -2 * _delta_sum;
3.2035 + (*_blossom_data)[surface].offset = 0;
3.2036 + (*_blossom_data)[surface].status = EVEN;
3.2037 + (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
3.2038 + (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
3.2039 +
3.2040 + _tree_set->insert(surface, tree);
3.2041 + _tree_set->erase(nca);
3.2042 + }
3.2043 +
3.2044 + void splitBlossom(int blossom) {
3.2045 + Edge next = (*_blossom_data)[blossom].next;
3.2046 + Edge pred = (*_blossom_data)[blossom].pred;
3.2047 +
3.2048 + int tree = _tree_set->find(blossom);
3.2049 +
3.2050 + (*_blossom_data)[blossom].status = MATCHED;
3.2051 + oddToMatched(blossom);
3.2052 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
3.2053 + _delta2->erase(blossom);
3.2054 + }
3.2055 +
3.2056 + std::vector<int> subblossoms;
3.2057 + _blossom_set->split(blossom, std::back_inserter(subblossoms));
3.2058 +
3.2059 + Value offset = (*_blossom_data)[blossom].offset;
3.2060 + int b = _blossom_set->find(_ugraph.source(pred));
3.2061 + int d = _blossom_set->find(_ugraph.source(next));
3.2062 +
3.2063 + int ib, id;
3.2064 + for (int i = 0; i < int(subblossoms.size()); ++i) {
3.2065 + if (subblossoms[i] == b) ib = i;
3.2066 + if (subblossoms[i] == d) id = i;
3.2067 +
3.2068 + (*_blossom_data)[subblossoms[i]].offset = offset;
3.2069 + if (!_blossom_set->trivial(subblossoms[i])) {
3.2070 + (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
3.2071 + }
3.2072 + if (_blossom_set->classPrio(subblossoms[i]) !=
3.2073 + std::numeric_limits<Value>::max()) {
3.2074 + _delta2->push(subblossoms[i],
3.2075 + _blossom_set->classPrio(subblossoms[i]) -
3.2076 + (*_blossom_data)[subblossoms[i]].offset);
3.2077 + }
3.2078 + }
3.2079 +
3.2080 + if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
3.2081 + for (int i = (id + 1) % subblossoms.size();
3.2082 + i != ib; i = (i + 2) % subblossoms.size()) {
3.2083 + int sb = subblossoms[i];
3.2084 + int tb = subblossoms[(i + 1) % subblossoms.size()];
3.2085 + (*_blossom_data)[sb].next =
3.2086 + _ugraph.oppositeEdge((*_blossom_data)[tb].next);
3.2087 + }
3.2088 +
3.2089 + for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
3.2090 + int sb = subblossoms[i];
3.2091 + int tb = subblossoms[(i + 1) % subblossoms.size()];
3.2092 + int ub = subblossoms[(i + 2) % subblossoms.size()];
3.2093 +
3.2094 + (*_blossom_data)[sb].status = ODD;
3.2095 + matchedToOdd(sb);
3.2096 + _tree_set->insert(sb, tree);
3.2097 + (*_blossom_data)[sb].pred = pred;
3.2098 + (*_blossom_data)[sb].next =
3.2099 + _ugraph.oppositeEdge((*_blossom_data)[tb].next);
3.2100 +
3.2101 + pred = (*_blossom_data)[ub].next;
3.2102 +
3.2103 + (*_blossom_data)[tb].status = EVEN;
3.2104 + matchedToEven(tb, tree);
3.2105 + _tree_set->insert(tb, tree);
3.2106 + (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
3.2107 + }
3.2108 +
3.2109 + (*_blossom_data)[subblossoms[id]].status = ODD;
3.2110 + matchedToOdd(subblossoms[id]);
3.2111 + _tree_set->insert(subblossoms[id], tree);
3.2112 + (*_blossom_data)[subblossoms[id]].next = next;
3.2113 + (*_blossom_data)[subblossoms[id]].pred = pred;
3.2114 +
3.2115 + } else {
3.2116 +
3.2117 + for (int i = (ib + 1) % subblossoms.size();
3.2118 + i != id; i = (i + 2) % subblossoms.size()) {
3.2119 + int sb = subblossoms[i];
3.2120 + int tb = subblossoms[(i + 1) % subblossoms.size()];
3.2121 + (*_blossom_data)[sb].next =
3.2122 + _ugraph.oppositeEdge((*_blossom_data)[tb].next);
3.2123 + }
3.2124 +
3.2125 + for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
3.2126 + int sb = subblossoms[i];
3.2127 + int tb = subblossoms[(i + 1) % subblossoms.size()];
3.2128 + int ub = subblossoms[(i + 2) % subblossoms.size()];
3.2129 +
3.2130 + (*_blossom_data)[sb].status = ODD;
3.2131 + matchedToOdd(sb);
3.2132 + _tree_set->insert(sb, tree);
3.2133 + (*_blossom_data)[sb].next = next;
3.2134 + (*_blossom_data)[sb].pred =
3.2135 + _ugraph.oppositeEdge((*_blossom_data)[tb].next);
3.2136 +
3.2137 + (*_blossom_data)[tb].status = EVEN;
3.2138 + matchedToEven(tb, tree);
3.2139 + _tree_set->insert(tb, tree);
3.2140 + (*_blossom_data)[tb].pred =
3.2141 + (*_blossom_data)[tb].next =
3.2142 + _ugraph.oppositeEdge((*_blossom_data)[ub].next);
3.2143 + next = (*_blossom_data)[ub].next;
3.2144 + }
3.2145 +
3.2146 + (*_blossom_data)[subblossoms[ib]].status = ODD;
3.2147 + matchedToOdd(subblossoms[ib]);
3.2148 + _tree_set->insert(subblossoms[ib], tree);
3.2149 + (*_blossom_data)[subblossoms[ib]].next = next;
3.2150 + (*_blossom_data)[subblossoms[ib]].pred = pred;
3.2151 + }
3.2152 + _tree_set->erase(blossom);
3.2153 + }
3.2154 +
3.2155 + void extractBlossom(int blossom, const Node& base, const Edge& matching) {
3.2156 + if (_blossom_set->trivial(blossom)) {
3.2157 + int bi = (*_node_index)[base];
3.2158 + Value pot = (*_node_data)[bi].pot;
3.2159 +
3.2160 + _matching->set(base, matching);
3.2161 + _blossom_node_list.push_back(base);
3.2162 + _node_potential->set(base, pot);
3.2163 + } else {
3.2164 +
3.2165 + Value pot = (*_blossom_data)[blossom].pot;
3.2166 + int bn = _blossom_node_list.size();
3.2167 +
3.2168 + std::vector<int> subblossoms;
3.2169 + _blossom_set->split(blossom, std::back_inserter(subblossoms));
3.2170 + int b = _blossom_set->find(base);
3.2171 + int ib = -1;
3.2172 + for (int i = 0; i < int(subblossoms.size()); ++i) {
3.2173 + if (subblossoms[i] == b) { ib = i; break; }
3.2174 + }
3.2175 +
3.2176 + for (int i = 1; i < int(subblossoms.size()); i += 2) {
3.2177 + int sb = subblossoms[(ib + i) % subblossoms.size()];
3.2178 + int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
3.2179 +
3.2180 + Edge m = (*_blossom_data)[tb].next;
3.2181 + extractBlossom(sb, _ugraph.target(m), _ugraph.oppositeEdge(m));
3.2182 + extractBlossom(tb, _ugraph.source(m), m);
3.2183 + }
3.2184 + extractBlossom(subblossoms[ib], base, matching);
3.2185 +
3.2186 + int en = _blossom_node_list.size();
3.2187 +
3.2188 + _blossom_potential.push_back(BlossomVariable(bn, en, pot));
3.2189 + }
3.2190 + }
3.2191 +
3.2192 + void extractMatching() {
3.2193 + std::vector<int> blossoms;
3.2194 + for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
3.2195 + blossoms.push_back(c);
3.2196 + }
3.2197 +
3.2198 + for (int i = 0; i < int(blossoms.size()); ++i) {
3.2199 +
3.2200 + Value offset = (*_blossom_data)[blossoms[i]].offset;
3.2201 + (*_blossom_data)[blossoms[i]].pot += 2 * offset;
3.2202 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
3.2203 + n != INVALID; ++n) {
3.2204 + (*_node_data)[(*_node_index)[n]].pot -= offset;
3.2205 + }
3.2206 +
3.2207 + Edge matching = (*_blossom_data)[blossoms[i]].next;
3.2208 + Node base = _ugraph.source(matching);
3.2209 + extractBlossom(blossoms[i], base, matching);
3.2210 + }
3.2211 + }
3.2212 +
3.2213 + public:
3.2214 +
3.2215 + /// \brief Constructor
3.2216 + ///
3.2217 + /// Constructor.
3.2218 + MaxWeightedPerfectMatching(const UGraph& ugraph, const WeightMap& weight)
3.2219 + : _ugraph(ugraph), _weight(weight), _matching(0),
3.2220 + _node_potential(0), _blossom_potential(), _blossom_node_list(),
3.2221 + _node_num(0), _blossom_num(0),
3.2222 +
3.2223 + _blossom_index(0), _blossom_set(0), _blossom_data(0),
3.2224 + _node_index(0), _node_heap_index(0), _node_data(0),
3.2225 + _tree_set_index(0), _tree_set(0),
3.2226 +
3.2227 + _delta2_index(0), _delta2(0),
3.2228 + _delta3_index(0), _delta3(0),
3.2229 + _delta4_index(0), _delta4(0),
3.2230 +
3.2231 + _delta_sum() {}
3.2232 +
3.2233 + ~MaxWeightedPerfectMatching() {
3.2234 + destroyStructures();
3.2235 + }
3.2236 +
3.2237 + /// \name Execution control
3.2238 + /// The simplest way to execute the algorithm is to use the member
3.2239 + /// \c run() member function.
3.2240 +
3.2241 + ///@{
3.2242 +
3.2243 + /// \brief Initialize the algorithm
3.2244 + ///
3.2245 + /// Initialize the algorithm
3.2246 + void init() {
3.2247 + createStructures();
3.2248 +
3.2249 + for (EdgeIt e(_ugraph); e != INVALID; ++e) {
3.2250 + _node_heap_index->set(e, BinHeap<Value, EdgeIntMap>::PRE_HEAP);
3.2251 + }
3.2252 + for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
3.2253 + _delta3_index->set(e, _delta3->PRE_HEAP);
3.2254 + }
3.2255 + for (int i = 0; i < _blossom_num; ++i) {
3.2256 + _delta2_index->set(i, _delta2->PRE_HEAP);
3.2257 + _delta4_index->set(i, _delta4->PRE_HEAP);
3.2258 + }
3.2259 +
3.2260 + int index = 0;
3.2261 + for (NodeIt n(_ugraph); n != INVALID; ++n) {
3.2262 + Value max = std::numeric_limits<Value>::min();
3.2263 + for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
3.2264 + if (_ugraph.target(e) == n) continue;
3.2265 + if ((dualScale * _weight[e]) / 2 > max) {
3.2266 + max = (dualScale * _weight[e]) / 2;
3.2267 + }
3.2268 + }
3.2269 + _node_index->set(n, index);
3.2270 + (*_node_data)[index].pot = max;
3.2271 + int blossom =
3.2272 + _blossom_set->insert(n, std::numeric_limits<Value>::max());
3.2273 +
3.2274 + _tree_set->insert(blossom);
3.2275 +
3.2276 + (*_blossom_data)[blossom].status = EVEN;
3.2277 + (*_blossom_data)[blossom].pred = INVALID;
3.2278 + (*_blossom_data)[blossom].next = INVALID;
3.2279 + (*_blossom_data)[blossom].pot = 0;
3.2280 + (*_blossom_data)[blossom].offset = 0;
3.2281 + ++index;
3.2282 + }
3.2283 + for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
3.2284 + int si = (*_node_index)[_ugraph.source(e)];
3.2285 + int ti = (*_node_index)[_ugraph.target(e)];
3.2286 + if (_ugraph.source(e) != _ugraph.target(e)) {
3.2287 + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
3.2288 + dualScale * _weight[e]) / 2);
3.2289 + }
3.2290 + }
3.2291 + }
3.2292 +
3.2293 + /// \brief Starts the algorithm
3.2294 + ///
3.2295 + /// Starts the algorithm
3.2296 + bool start() {
3.2297 + enum OpType {
3.2298 + D2, D3, D4
3.2299 + };
3.2300 +
3.2301 + int unmatched = _node_num;
3.2302 + while (unmatched > 0) {
3.2303 + Value d2 = !_delta2->empty() ?
3.2304 + _delta2->prio() : std::numeric_limits<Value>::max();
3.2305 +
3.2306 + Value d3 = !_delta3->empty() ?
3.2307 + _delta3->prio() : std::numeric_limits<Value>::max();
3.2308 +
3.2309 + Value d4 = !_delta4->empty() ?
3.2310 + _delta4->prio() : std::numeric_limits<Value>::max();
3.2311 +
3.2312 + _delta_sum = d2; OpType ot = D2;
3.2313 + if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
3.2314 + if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
3.2315 +
3.2316 + if (_delta_sum == std::numeric_limits<Value>::max()) {
3.2317 + return false;
3.2318 + }
3.2319 +
3.2320 + switch (ot) {
3.2321 + case D2:
3.2322 + {
3.2323 + int blossom = _delta2->top();
3.2324 + Node n = _blossom_set->classTop(blossom);
3.2325 + Edge e = (*_node_data)[(*_node_index)[n]].heap.top();
3.2326 + extendOnEdge(e);
3.2327 + }
3.2328 + break;
3.2329 + case D3:
3.2330 + {
3.2331 + UEdge e = _delta3->top();
3.2332 +
3.2333 + int left_blossom = _blossom_set->find(_ugraph.source(e));
3.2334 + int right_blossom = _blossom_set->find(_ugraph.target(e));
3.2335 +
3.2336 + if (left_blossom == right_blossom) {
3.2337 + _delta3->pop();
3.2338 + } else {
3.2339 + int left_tree = _tree_set->find(left_blossom);
3.2340 + int right_tree = _tree_set->find(right_blossom);
3.2341 +
3.2342 + if (left_tree == right_tree) {
3.2343 + shrinkOnEdge(e, left_tree);
3.2344 + } else {
3.2345 + augmentOnEdge(e);
3.2346 + unmatched -= 2;
3.2347 + }
3.2348 + }
3.2349 + } break;
3.2350 + case D4:
3.2351 + splitBlossom(_delta4->top());
3.2352 + break;
3.2353 + }
3.2354 + }
3.2355 + extractMatching();
3.2356 + return true;
3.2357 + }
3.2358 +
3.2359 + /// \brief Runs %MaxWeightedPerfectMatching algorithm.
3.2360 + ///
3.2361 + /// This method runs the %MaxWeightedPerfectMatching algorithm.
3.2362 + ///
3.2363 + /// \note mwm.run() is just a shortcut of the following code.
3.2364 + /// \code
3.2365 + /// mwm.init();
3.2366 + /// mwm.start();
3.2367 + /// \endcode
3.2368 + bool run() {
3.2369 + init();
3.2370 + return start();
3.2371 + }
3.2372 +
3.2373 + /// @}
3.2374 +
3.2375 + /// \name Primal solution
3.2376 + /// Functions for get the primal solution, ie. the matching.
3.2377 +
3.2378 + /// @{
3.2379 +
3.2380 + /// \brief Returns the matching value.
3.2381 + ///
3.2382 + /// Returns the matching value.
3.2383 + Value matchingValue() const {
3.2384 + Value sum = 0;
3.2385 + for (NodeIt n(_ugraph); n != INVALID; ++n) {
3.2386 + if ((*_matching)[n] != INVALID) {
3.2387 + sum += _weight[(*_matching)[n]];
3.2388 + }
3.2389 + }
3.2390 + return sum /= 2;
3.2391 + }
3.2392 +
3.2393 + /// \brief Returns true when the edge is in the matching.
3.2394 + ///
3.2395 + /// Returns true when the edge is in the matching.
3.2396 + bool matching(const UEdge& edge) const {
3.2397 + return (*_matching)[_ugraph.source(edge)] == _ugraph.direct(edge, true);
3.2398 + }
3.2399 +
3.2400 + /// \brief Returns the incident matching edge.
3.2401 + ///
3.2402 + /// Returns the incident matching edge from given node.
3.2403 + Edge matching(const Node& node) const {
3.2404 + return (*_matching)[node];
3.2405 + }
3.2406 +
3.2407 + /// \brief Returns the mate of the node.
3.2408 + ///
3.2409 + /// Returns the adjancent node in a mathcing edge.
3.2410 + Node mate(const Node& node) const {
3.2411 + return _ugraph.target((*_matching)[node]);
3.2412 + }
3.2413 +
3.2414 + /// @}
3.2415 +
3.2416 + /// \name Dual solution
3.2417 + /// Functions for get the dual solution.
3.2418 +
3.2419 + /// @{
3.2420 +
3.2421 + /// \brief Returns the value of the dual solution.
3.2422 + ///
3.2423 + /// Returns the value of the dual solution. It should be equal to
3.2424 + /// the primal value scaled by \ref dualScale "dual scale".
3.2425 + Value dualValue() const {
3.2426 + Value sum = 0;
3.2427 + for (NodeIt n(_ugraph); n != INVALID; ++n) {
3.2428 + sum += nodeValue(n);
3.2429 + }
3.2430 + for (int i = 0; i < blossomNum(); ++i) {
3.2431 + sum += blossomValue(i) * (blossomSize(i) / 2);
3.2432 + }
3.2433 + return sum;
3.2434 + }
3.2435 +
3.2436 + /// \brief Returns the value of the node.
3.2437 + ///
3.2438 + /// Returns the the value of the node.
3.2439 + Value nodeValue(const Node& n) const {
3.2440 + return (*_node_potential)[n];
3.2441 + }
3.2442 +
3.2443 + /// \brief Returns the number of the blossoms in the basis.
3.2444 + ///
3.2445 + /// Returns the number of the blossoms in the basis.
3.2446 + /// \see BlossomIt
3.2447 + int blossomNum() const {
3.2448 + return _blossom_potential.size();
3.2449 + }
3.2450 +
3.2451 +
3.2452 + /// \brief Returns the number of the nodes in the blossom.
3.2453 + ///
3.2454 + /// Returns the number of the nodes in the blossom.
3.2455 + int blossomSize(int k) const {
3.2456 + return _blossom_potential[k].end - _blossom_potential[k].begin;
3.2457 + }
3.2458 +
3.2459 + /// \brief Returns the value of the blossom.
3.2460 + ///
3.2461 + /// Returns the the value of the blossom.
3.2462 + /// \see BlossomIt
3.2463 + Value blossomValue(int k) const {
3.2464 + return _blossom_potential[k].value;
3.2465 + }
3.2466 +
3.2467 + /// \brief Lemon iterator for get the items of the blossom.
3.2468 + ///
3.2469 + /// Lemon iterator for get the nodes of the blossom. This class
3.2470 + /// provides a common style lemon iterator which gives back a
3.2471 + /// subset of the nodes.
3.2472 + class BlossomIt {
3.2473 + public:
3.2474 +
3.2475 + /// \brief Constructor.
3.2476 + ///
3.2477 + /// Constructor for get the nodes of the variable.
3.2478 + BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
3.2479 + : _algorithm(&algorithm)
3.2480 + {
3.2481 + _index = _algorithm->_blossom_potential[variable].begin;
3.2482 + _last = _algorithm->_blossom_potential[variable].end;
3.2483 + }
3.2484 +
3.2485 + /// \brief Invalid constructor.
3.2486 + ///
3.2487 + /// Invalid constructor.
3.2488 + BlossomIt(Invalid) : _index(-1) {}
3.2489 +
3.2490 + /// \brief Conversion to node.
3.2491 + ///
3.2492 + /// Conversion to node.
3.2493 + operator Node() const {
3.2494 + return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
3.2495 + }
3.2496 +
3.2497 + /// \brief Increment operator.
3.2498 + ///
3.2499 + /// Increment operator.
3.2500 + BlossomIt& operator++() {
3.2501 + ++_index;
3.2502 + if (_index == _last) {
3.2503 + _index = -1;
3.2504 + }
3.2505 + return *this;
3.2506 + }
3.2507 +
3.2508 + bool operator==(const BlossomIt& it) const {
3.2509 + return _index == it._index;
3.2510 + }
3.2511 + bool operator!=(const BlossomIt& it) const {
3.2512 + return _index != it._index;
3.2513 + }
3.2514 +
3.2515 + private:
3.2516 + const MaxWeightedPerfectMatching* _algorithm;
3.2517 + int _last;
3.2518 + int _index;
3.2519 + };
3.2520 +
3.2521 + /// @}
3.2522 +
3.2523 + };
3.2524 +
3.2525
3.2526 } //END OF NAMESPACE LEMON
3.2527
4.1 --- a/lemon/unionfind.h Thu Dec 27 13:40:16 2007 +0000
4.2 +++ b/lemon/unionfind.h Fri Dec 28 11:00:51 2007 +0000
4.3 @@ -924,6 +924,869 @@
4.4
4.5 };
4.6
4.7 + /// \ingroup auxdat
4.8 + ///
4.9 + /// \brief A \e Union-Find data structure implementation which
4.10 + /// is able to store a priority for each item and retrieve the minimum of
4.11 + /// each class.
4.12 + ///
4.13 + /// A \e Union-Find data structure implementation which is able to
4.14 + /// store a priority for each item and retrieve the minimum of each
4.15 + /// class. In addition, it supports the joining and splitting the
4.16 + /// components. If you don't need this feature then you makes
4.17 + /// better to use the \ref UnionFind class which is more efficient.
4.18 + ///
4.19 + /// The union-find data strcuture based on a (2, 16)-tree with a
4.20 + /// tournament minimum selection on the internal nodes. The insert
4.21 + /// operation takes O(1), the find, set, decrease and increase takes
4.22 + /// O(log(n)), where n is the number of nodes in the current
4.23 + /// component. The complexity of join and split is O(log(n)*k),
4.24 + /// where n is the sum of the number of the nodes and k is the
4.25 + /// number of joined components or the number of the components
4.26 + /// after the split.
4.27 + ///
4.28 + /// \pre You need to add all the elements by the \ref insert()
4.29 + /// method.
4.30 + ///
4.31 + template <typename _Value, typename _ItemIntMap,
4.32 + typename _Comp = std::less<_Value> >
4.33 + class HeapUnionFind {
4.34 + public:
4.35 +
4.36 + typedef _Value Value;
4.37 + typedef typename _ItemIntMap::Key Item;
4.38 +
4.39 + typedef _ItemIntMap ItemIntMap;
4.40 +
4.41 + typedef _Comp Comp;
4.42 +
4.43 + private:
4.44 +
4.45 + static const int cmax = 3;
4.46 +
4.47 + ItemIntMap& index;
4.48 +
4.49 + struct ClassNode {
4.50 + int parent;
4.51 + int depth;
4.52 +
4.53 + int left, right;
4.54 + int next, prev;
4.55 + };
4.56 +
4.57 + int first_class;
4.58 + int first_free_class;
4.59 + std::vector<ClassNode> classes;
4.60 +
4.61 + int newClass() {
4.62 + if (first_free_class < 0) {
4.63 + int id = classes.size();
4.64 + classes.push_back(ClassNode());
4.65 + return id;
4.66 + } else {
4.67 + int id = first_free_class;
4.68 + first_free_class = classes[id].next;
4.69 + return id;
4.70 + }
4.71 + }
4.72 +
4.73 + void deleteClass(int id) {
4.74 + classes[id].next = first_free_class;
4.75 + first_free_class = id;
4.76 + }
4.77 +
4.78 + struct ItemNode {
4.79 + int parent;
4.80 + Item item;
4.81 + Value prio;
4.82 + int next, prev;
4.83 + int left, right;
4.84 + int size;
4.85 + };
4.86 +
4.87 + int first_free_node;
4.88 + std::vector<ItemNode> nodes;
4.89 +
4.90 + int newNode() {
4.91 + if (first_free_node < 0) {
4.92 + int id = nodes.size();
4.93 + nodes.push_back(ItemNode());
4.94 + return id;
4.95 + } else {
4.96 + int id = first_free_node;
4.97 + first_free_node = nodes[id].next;
4.98 + return id;
4.99 + }
4.100 + }
4.101 +
4.102 + void deleteNode(int id) {
4.103 + nodes[id].next = first_free_node;
4.104 + first_free_node = id;
4.105 + }
4.106 +
4.107 + Comp comp;
4.108 +
4.109 + int findClass(int id) const {
4.110 + int kd = id;
4.111 + while (kd >= 0) {
4.112 + kd = nodes[kd].parent;
4.113 + }
4.114 + return ~kd;
4.115 + }
4.116 +
4.117 + int leftNode(int id) const {
4.118 + int kd = ~(classes[id].parent);
4.119 + for (int i = 0; i < classes[id].depth; ++i) {
4.120 + kd = nodes[kd].left;
4.121 + }
4.122 + return kd;
4.123 + }
4.124 +
4.125 + int nextNode(int id) const {
4.126 + int depth = 0;
4.127 + while (id >= 0 && nodes[id].next == -1) {
4.128 + id = nodes[id].parent;
4.129 + ++depth;
4.130 + }
4.131 + if (id < 0) {
4.132 + return -1;
4.133 + }
4.134 + id = nodes[id].next;
4.135 + while (depth--) {
4.136 + id = nodes[id].left;
4.137 + }
4.138 + return id;
4.139 + }
4.140 +
4.141 +
4.142 + void setPrio(int id) {
4.143 + int jd = nodes[id].left;
4.144 + nodes[id].prio = nodes[jd].prio;
4.145 + nodes[id].item = nodes[jd].item;
4.146 + jd = nodes[jd].next;
4.147 + while (jd != -1) {
4.148 + if (comp(nodes[jd].prio, nodes[id].prio)) {
4.149 + nodes[id].prio = nodes[jd].prio;
4.150 + nodes[id].item = nodes[jd].item;
4.151 + }
4.152 + jd = nodes[jd].next;
4.153 + }
4.154 + }
4.155 +
4.156 + void push(int id, int jd) {
4.157 + nodes[id].size = 1;
4.158 + nodes[id].left = nodes[id].right = jd;
4.159 + nodes[jd].next = nodes[jd].prev = -1;
4.160 + nodes[jd].parent = id;
4.161 + }
4.162 +
4.163 + void pushAfter(int id, int jd) {
4.164 + int kd = nodes[id].parent;
4.165 + if (nodes[id].next != -1) {
4.166 + nodes[nodes[id].next].prev = jd;
4.167 + if (kd >= 0) {
4.168 + nodes[kd].size += 1;
4.169 + }
4.170 + } else {
4.171 + if (kd >= 0) {
4.172 + nodes[kd].right = jd;
4.173 + nodes[kd].size += 1;
4.174 + }
4.175 + }
4.176 + nodes[jd].next = nodes[id].next;
4.177 + nodes[jd].prev = id;
4.178 + nodes[id].next = jd;
4.179 + nodes[jd].parent = kd;
4.180 + }
4.181 +
4.182 + void pushRight(int id, int jd) {
4.183 + nodes[id].size += 1;
4.184 + nodes[jd].prev = nodes[id].right;
4.185 + nodes[jd].next = -1;
4.186 + nodes[nodes[id].right].next = jd;
4.187 + nodes[id].right = jd;
4.188 + nodes[jd].parent = id;
4.189 + }
4.190 +
4.191 + void popRight(int id) {
4.192 + nodes[id].size -= 1;
4.193 + int jd = nodes[id].right;
4.194 + nodes[nodes[jd].prev].next = -1;
4.195 + nodes[id].right = nodes[jd].prev;
4.196 + }
4.197 +
4.198 + void splice(int id, int jd) {
4.199 + nodes[id].size += nodes[jd].size;
4.200 + nodes[nodes[id].right].next = nodes[jd].left;
4.201 + nodes[nodes[jd].left].prev = nodes[id].right;
4.202 + int kd = nodes[jd].left;
4.203 + while (kd != -1) {
4.204 + nodes[kd].parent = id;
4.205 + kd = nodes[kd].next;
4.206 + }
4.207 + }
4.208 +
4.209 + void split(int id, int jd) {
4.210 + int kd = nodes[id].parent;
4.211 + nodes[kd].right = nodes[id].prev;
4.212 + nodes[nodes[id].prev].next = -1;
4.213 +
4.214 + nodes[jd].left = id;
4.215 + nodes[id].prev = -1;
4.216 + int num = 0;
4.217 + while (id != -1) {
4.218 + nodes[id].parent = jd;
4.219 + nodes[jd].right = id;
4.220 + id = nodes[id].next;
4.221 + ++num;
4.222 + }
4.223 + nodes[kd].size -= num;
4.224 + nodes[jd].size = num;
4.225 + }
4.226 +
4.227 + void pushLeft(int id, int jd) {
4.228 + nodes[id].size += 1;
4.229 + nodes[jd].next = nodes[id].left;
4.230 + nodes[jd].prev = -1;
4.231 + nodes[nodes[id].left].prev = jd;
4.232 + nodes[id].left = jd;
4.233 + nodes[jd].parent = id;
4.234 + }
4.235 +
4.236 + void popLeft(int id) {
4.237 + nodes[id].size -= 1;
4.238 + int jd = nodes[id].left;
4.239 + nodes[nodes[jd].next].prev = -1;
4.240 + nodes[id].left = nodes[jd].next;
4.241 + }
4.242 +
4.243 + void repairLeft(int id) {
4.244 + int jd = ~(classes[id].parent);
4.245 + while (nodes[jd].left != -1) {
4.246 + int kd = nodes[jd].left;
4.247 + if (nodes[jd].size == 1) {
4.248 + if (nodes[jd].parent < 0) {
4.249 + classes[id].parent = ~kd;
4.250 + classes[id].depth -= 1;
4.251 + nodes[kd].parent = ~id;
4.252 + deleteNode(jd);
4.253 + jd = kd;
4.254 + } else {
4.255 + int pd = nodes[jd].parent;
4.256 + if (nodes[nodes[jd].next].size < cmax) {
4.257 + pushLeft(nodes[jd].next, nodes[jd].left);
4.258 + if (less(nodes[jd].left, nodes[jd].next)) {
4.259 + nodes[nodes[jd].next].prio = nodes[nodes[jd].left].prio;
4.260 + nodes[nodes[jd].next].item = nodes[nodes[jd].left].item;
4.261 + }
4.262 + popLeft(pd);
4.263 + deleteNode(jd);
4.264 + jd = pd;
4.265 + } else {
4.266 + int ld = nodes[nodes[jd].next].left;
4.267 + popLeft(nodes[jd].next);
4.268 + pushRight(jd, ld);
4.269 + if (less(ld, nodes[jd].left)) {
4.270 + nodes[jd].item = nodes[ld].item;
4.271 + nodes[jd].prio = nodes[jd].prio;
4.272 + }
4.273 + if (nodes[nodes[jd].next].item == nodes[ld].item) {
4.274 + setPrio(nodes[jd].next);
4.275 + }
4.276 + jd = nodes[jd].left;
4.277 + }
4.278 + }
4.279 + } else {
4.280 + jd = nodes[jd].left;
4.281 + }
4.282 + }
4.283 + }
4.284 +
4.285 + void repairRight(int id) {
4.286 + int jd = ~(classes[id].parent);
4.287 + while (nodes[jd].right != -1) {
4.288 + int kd = nodes[jd].right;
4.289 + if (nodes[jd].size == 1) {
4.290 + if (nodes[jd].parent < 0) {
4.291 + classes[id].parent = ~kd;
4.292 + classes[id].depth -= 1;
4.293 + nodes[kd].parent = ~id;
4.294 + deleteNode(jd);
4.295 + jd = kd;
4.296 + } else {
4.297 + int pd = nodes[jd].parent;
4.298 + if (nodes[nodes[jd].prev].size < cmax) {
4.299 + pushRight(nodes[jd].prev, nodes[jd].right);
4.300 + if (less(nodes[jd].right, nodes[jd].prev)) {
4.301 + nodes[nodes[jd].prev].prio = nodes[nodes[jd].right].prio;
4.302 + nodes[nodes[jd].prev].item = nodes[nodes[jd].right].item;
4.303 + }
4.304 + popRight(pd);
4.305 + deleteNode(jd);
4.306 + jd = pd;
4.307 + } else {
4.308 + int ld = nodes[nodes[jd].prev].right;
4.309 + popRight(nodes[jd].prev);
4.310 + pushLeft(jd, ld);
4.311 + if (less(ld, nodes[jd].right)) {
4.312 + nodes[jd].item = nodes[ld].item;
4.313 + nodes[jd].prio = nodes[jd].prio;
4.314 + }
4.315 + if (nodes[nodes[jd].prev].item == nodes[ld].item) {
4.316 + setPrio(nodes[jd].prev);
4.317 + }
4.318 + jd = nodes[jd].right;
4.319 + }
4.320 + }
4.321 + } else {
4.322 + jd = nodes[jd].right;
4.323 + }
4.324 + }
4.325 + }
4.326 +
4.327 +
4.328 + bool less(int id, int jd) const {
4.329 + return comp(nodes[id].prio, nodes[jd].prio);
4.330 + }
4.331 +
4.332 + bool equal(int id, int jd) const {
4.333 + return !less(id, jd) && !less(jd, id);
4.334 + }
4.335 +
4.336 +
4.337 + public:
4.338 +
4.339 + /// \brief Returns true when the given class is alive.
4.340 + ///
4.341 + /// Returns true when the given class is alive, ie. the class is
4.342 + /// not nested into other class.
4.343 + bool alive(int cls) const {
4.344 + return classes[cls].parent < 0;
4.345 + }
4.346 +
4.347 + /// \brief Returns true when the given class is trivial.
4.348 + ///
4.349 + /// Returns true when the given class is trivial, ie. the class
4.350 + /// contains just one item directly.
4.351 + bool trivial(int cls) const {
4.352 + return classes[cls].left == -1;
4.353 + }
4.354 +
4.355 + /// \brief Constructs the union-find.
4.356 + ///
4.357 + /// Constructs the union-find.
4.358 + /// \brief _index The index map of the union-find. The data
4.359 + /// structure uses internally for store references.
4.360 + HeapUnionFind(ItemIntMap& _index)
4.361 + : index(_index), first_class(-1),
4.362 + first_free_class(-1), first_free_node(-1) {}
4.363 +
4.364 + /// \brief Insert a new node into a new component.
4.365 + ///
4.366 + /// Insert a new node into a new component.
4.367 + /// \param item The item of the new node.
4.368 + /// \param prio The priority of the new node.
4.369 + /// \return The class id of the one-item-heap.
4.370 + int insert(const Item& item, const Value& prio) {
4.371 + int id = newNode();
4.372 + nodes[id].item = item;
4.373 + nodes[id].prio = prio;
4.374 + nodes[id].size = 0;
4.375 +
4.376 + nodes[id].prev = -1;
4.377 + nodes[id].next = -1;
4.378 +
4.379 + nodes[id].left = -1;
4.380 + nodes[id].right = -1;
4.381 +
4.382 + nodes[id].item = item;
4.383 + index[item] = id;
4.384 +
4.385 + int class_id = newClass();
4.386 + classes[class_id].parent = ~id;
4.387 + classes[class_id].depth = 0;
4.388 +
4.389 + classes[class_id].left = -1;
4.390 + classes[class_id].right = -1;
4.391 +
4.392 + if (first_class != -1) {
4.393 + classes[first_class].prev = class_id;
4.394 + }
4.395 + classes[class_id].next = first_class;
4.396 + classes[class_id].prev = -1;
4.397 + first_class = class_id;
4.398 +
4.399 + nodes[id].parent = ~class_id;
4.400 +
4.401 + return class_id;
4.402 + }
4.403 +
4.404 + /// \brief The class of the item.
4.405 + ///
4.406 + /// \return The alive class id of the item, which is not nested into
4.407 + /// other classes.
4.408 + ///
4.409 + /// The time complexity is O(log(n)).
4.410 + int find(const Item& item) const {
4.411 + return findClass(index[item]);
4.412 + }
4.413 +
4.414 + /// \brief Joins the classes.
4.415 + ///
4.416 + /// The current function joins the given classes. The parameter is
4.417 + /// an STL range which should be contains valid class ids. The
4.418 + /// time complexity is O(log(n)*k) where n is the overall number
4.419 + /// of the joined nodes and k is the number of classes.
4.420 + /// \return The class of the joined classes.
4.421 + /// \pre The range should contain at least two class ids.
4.422 + template <typename Iterator>
4.423 + int join(Iterator begin, Iterator end) {
4.424 + std::vector<int> cs;
4.425 + for (Iterator it = begin; it != end; ++it) {
4.426 + cs.push_back(*it);
4.427 + }
4.428 +
4.429 + int class_id = newClass();
4.430 + { // creation union-find
4.431 +
4.432 + if (first_class != -1) {
4.433 + classes[first_class].prev = class_id;
4.434 + }
4.435 + classes[class_id].next = first_class;
4.436 + classes[class_id].prev = -1;
4.437 + first_class = class_id;
4.438 +
4.439 + classes[class_id].depth = classes[cs[0]].depth;
4.440 + classes[class_id].parent = classes[cs[0]].parent;
4.441 + nodes[~(classes[class_id].parent)].parent = ~class_id;
4.442 +
4.443 + int l = cs[0];
4.444 +
4.445 + classes[class_id].left = l;
4.446 + classes[class_id].right = l;
4.447 +
4.448 + if (classes[l].next != -1) {
4.449 + classes[classes[l].next].prev = classes[l].prev;
4.450 + }
4.451 + classes[classes[l].prev].next = classes[l].next;
4.452 +
4.453 + classes[l].prev = -1;
4.454 + classes[l].next = -1;
4.455 +
4.456 + classes[l].depth = leftNode(l);
4.457 + classes[l].parent = class_id;
4.458 +
4.459 + }
4.460 +
4.461 + { // merging of heap
4.462 + int l = class_id;
4.463 + for (int ci = 1; ci < int(cs.size()); ++ci) {
4.464 + int r = cs[ci];
4.465 + int rln = leftNode(r);
4.466 + if (classes[l].depth > classes[r].depth) {
4.467 + int id = ~(classes[l].parent);
4.468 + for (int i = classes[r].depth + 1; i < classes[l].depth; ++i) {
4.469 + id = nodes[id].right;
4.470 + }
4.471 + while (id >= 0 && nodes[id].size == cmax) {
4.472 + int new_id = newNode();
4.473 + int right_id = nodes[id].right;
4.474 +
4.475 + popRight(id);
4.476 + if (nodes[id].item == nodes[right_id].item) {
4.477 + setPrio(id);
4.478 + }
4.479 + push(new_id, right_id);
4.480 + pushRight(new_id, ~(classes[r].parent));
4.481 + setPrio(new_id);
4.482 +
4.483 + id = nodes[id].parent;
4.484 + classes[r].parent = ~new_id;
4.485 + }
4.486 + if (id < 0) {
4.487 + int new_parent = newNode();
4.488 + nodes[new_parent].next = -1;
4.489 + nodes[new_parent].prev = -1;
4.490 + nodes[new_parent].parent = ~l;
4.491 +
4.492 + push(new_parent, ~(classes[l].parent));
4.493 + pushRight(new_parent, ~(classes[r].parent));
4.494 + setPrio(new_parent);
4.495 +
4.496 + classes[l].parent = ~new_parent;
4.497 + classes[l].depth += 1;
4.498 + } else {
4.499 + pushRight(id, ~(classes[r].parent));
4.500 + while (id >= 0 && less(~(classes[r].parent), id)) {
4.501 + nodes[id].prio = nodes[~(classes[r].parent)].prio;
4.502 + nodes[id].item = nodes[~(classes[r].parent)].item;
4.503 + id = nodes[id].parent;
4.504 + }
4.505 + }
4.506 + } else if (classes[r].depth > classes[l].depth) {
4.507 + int id = ~(classes[r].parent);
4.508 + for (int i = classes[l].depth + 1; i < classes[r].depth; ++i) {
4.509 + id = nodes[id].left;
4.510 + }
4.511 + while (id >= 0 && nodes[id].size == cmax) {
4.512 + int new_id = newNode();
4.513 + int left_id = nodes[id].left;
4.514 +
4.515 + popLeft(id);
4.516 + if (nodes[id].prio == nodes[left_id].prio) {
4.517 + setPrio(id);
4.518 + }
4.519 + push(new_id, left_id);
4.520 + pushLeft(new_id, ~(classes[l].parent));
4.521 + setPrio(new_id);
4.522 +
4.523 + id = nodes[id].parent;
4.524 + classes[l].parent = ~new_id;
4.525 +
4.526 + }
4.527 + if (id < 0) {
4.528 + int new_parent = newNode();
4.529 + nodes[new_parent].next = -1;
4.530 + nodes[new_parent].prev = -1;
4.531 + nodes[new_parent].parent = ~l;
4.532 +
4.533 + push(new_parent, ~(classes[r].parent));
4.534 + pushLeft(new_parent, ~(classes[l].parent));
4.535 + setPrio(new_parent);
4.536 +
4.537 + classes[r].parent = ~new_parent;
4.538 + classes[r].depth += 1;
4.539 + } else {
4.540 + pushLeft(id, ~(classes[l].parent));
4.541 + while (id >= 0 && less(~(classes[l].parent), id)) {
4.542 + nodes[id].prio = nodes[~(classes[l].parent)].prio;
4.543 + nodes[id].item = nodes[~(classes[l].parent)].item;
4.544 + id = nodes[id].parent;
4.545 + }
4.546 + }
4.547 + nodes[~(classes[r].parent)].parent = ~l;
4.548 + classes[l].parent = classes[r].parent;
4.549 + classes[l].depth = classes[r].depth;
4.550 + } else {
4.551 + if (classes[l].depth != 0 &&
4.552 + nodes[~(classes[l].parent)].size +
4.553 + nodes[~(classes[r].parent)].size <= cmax) {
4.554 + splice(~(classes[l].parent), ~(classes[r].parent));
4.555 + deleteNode(~(classes[r].parent));
4.556 + if (less(~(classes[r].parent), ~(classes[l].parent))) {
4.557 + nodes[~(classes[l].parent)].prio =
4.558 + nodes[~(classes[r].parent)].prio;
4.559 + nodes[~(classes[l].parent)].item =
4.560 + nodes[~(classes[r].parent)].item;
4.561 + }
4.562 + } else {
4.563 + int new_parent = newNode();
4.564 + nodes[new_parent].next = nodes[new_parent].prev = -1;
4.565 + push(new_parent, ~(classes[l].parent));
4.566 + pushRight(new_parent, ~(classes[r].parent));
4.567 + setPrio(new_parent);
4.568 +
4.569 + classes[l].parent = ~new_parent;
4.570 + classes[l].depth += 1;
4.571 + nodes[new_parent].parent = ~l;
4.572 + }
4.573 + }
4.574 + if (classes[r].next != -1) {
4.575 + classes[classes[r].next].prev = classes[r].prev;
4.576 + }
4.577 + classes[classes[r].prev].next = classes[r].next;
4.578 +
4.579 + classes[r].prev = classes[l].right;
4.580 + classes[classes[l].right].next = r;
4.581 + classes[l].right = r;
4.582 + classes[r].parent = l;
4.583 +
4.584 + classes[r].next = -1;
4.585 + classes[r].depth = rln;
4.586 + }
4.587 + }
4.588 + return class_id;
4.589 + }
4.590 +
4.591 + /// \brief Split the class to subclasses.
4.592 + ///
4.593 + /// The current function splits the given class. The join, which
4.594 + /// made the current class, stored a reference to the
4.595 + /// subclasses. The \c splitClass() member restores the classes
4.596 + /// and creates the heaps. The parameter is an STL output iterator
4.597 + /// which will be filled with the subclass ids. The time
4.598 + /// complexity is O(log(n)*k) where n is the overall number of
4.599 + /// nodes in the splitted classes and k is the number of the
4.600 + /// classes.
4.601 + template <typename Iterator>
4.602 + void split(int cls, Iterator out) {
4.603 + std::vector<int> cs;
4.604 + { // splitting union-find
4.605 + int id = cls;
4.606 + int l = classes[id].left;
4.607 +
4.608 + classes[l].parent = classes[id].parent;
4.609 + classes[l].depth = classes[id].depth;
4.610 +
4.611 + nodes[~(classes[l].parent)].parent = ~l;
4.612 +
4.613 + *out++ = l;
4.614 +
4.615 + while (l != -1) {
4.616 + cs.push_back(l);
4.617 + l = classes[l].next;
4.618 + }
4.619 +
4.620 + classes[classes[id].right].next = first_class;
4.621 + classes[first_class].prev = classes[id].right;
4.622 + first_class = classes[id].left;
4.623 +
4.624 + if (classes[id].next != -1) {
4.625 + classes[classes[id].next].prev = classes[id].prev;
4.626 + }
4.627 + classes[classes[id].prev].next = classes[id].next;
4.628 +
4.629 + deleteClass(id);
4.630 + }
4.631 +
4.632 + {
4.633 + for (int i = 1; i < int(cs.size()); ++i) {
4.634 + int l = classes[cs[i]].depth;
4.635 + while (nodes[nodes[l].parent].left == l) {
4.636 + l = nodes[l].parent;
4.637 + }
4.638 + int r = l;
4.639 + while (nodes[l].parent >= 0) {
4.640 + l = nodes[l].parent;
4.641 + int new_node = newNode();
4.642 +
4.643 + nodes[new_node].prev = -1;
4.644 + nodes[new_node].next = -1;
4.645 +
4.646 + split(r, new_node);
4.647 + pushAfter(l, new_node);
4.648 + setPrio(l);
4.649 + setPrio(new_node);
4.650 + r = new_node;
4.651 + }
4.652 + classes[cs[i]].parent = ~r;
4.653 + classes[cs[i]].depth = classes[~(nodes[l].parent)].depth;
4.654 + nodes[r].parent = ~cs[i];
4.655 +
4.656 + nodes[l].next = -1;
4.657 + nodes[r].prev = -1;
4.658 +
4.659 + repairRight(~(nodes[l].parent));
4.660 + repairLeft(cs[i]);
4.661 +
4.662 + *out++ = cs[i];
4.663 + }
4.664 + }
4.665 + }
4.666 +
4.667 + /// \brief Gives back the priority of the current item.
4.668 + ///
4.669 + /// \return Gives back the priority of the current item.
4.670 + const Value& operator[](const Item& item) const {
4.671 + return nodes[index[item]].prio;
4.672 + }
4.673 +
4.674 + /// \brief Sets the priority of the current item.
4.675 + ///
4.676 + /// Sets the priority of the current item.
4.677 + void set(const Item& item, const Value& prio) {
4.678 + if (comp(prio, nodes[index[item]].prio)) {
4.679 + decrease(item, prio);
4.680 + } else if (!comp(prio, nodes[index[item]].prio)) {
4.681 + increase(item, prio);
4.682 + }
4.683 + }
4.684 +
4.685 + /// \brief Increase the priority of the current item.
4.686 + ///
4.687 + /// Increase the priority of the current item.
4.688 + void increase(const Item& item, const Value& prio) {
4.689 + int id = index[item];
4.690 + int kd = nodes[id].parent;
4.691 + nodes[id].prio = prio;
4.692 + while (kd >= 0 && nodes[kd].item == item) {
4.693 + setPrio(kd);
4.694 + kd = nodes[kd].parent;
4.695 + }
4.696 + }
4.697 +
4.698 + /// \brief Increase the priority of the current item.
4.699 + ///
4.700 + /// Increase the priority of the current item.
4.701 + void decrease(const Item& item, const Value& prio) {
4.702 + int id = index[item];
4.703 + int kd = nodes[id].parent;
4.704 + nodes[id].prio = prio;
4.705 + while (kd >= 0 && less(id, kd)) {
4.706 + nodes[kd].prio = prio;
4.707 + nodes[kd].item = item;
4.708 + kd = nodes[kd].parent;
4.709 + }
4.710 + }
4.711 +
4.712 + /// \brief Gives back the minimum priority of the class.
4.713 + ///
4.714 + /// \return Gives back the minimum priority of the class.
4.715 + const Value& classPrio(int cls) const {
4.716 + return nodes[~(classes[cls].parent)].prio;
4.717 + }
4.718 +
4.719 + /// \brief Gives back the minimum priority item of the class.
4.720 + ///
4.721 + /// \return Gives back the minimum priority item of the class.
4.722 + const Item& classTop(int cls) const {
4.723 + return nodes[~(classes[cls].parent)].item;
4.724 + }
4.725 +
4.726 + /// \brief Gives back a representant item of the class.
4.727 + ///
4.728 + /// The representant is indpendent from the priorities of the
4.729 + /// items.
4.730 + /// \return Gives back a representant item of the class.
4.731 + const Item& classRep(int id) const {
4.732 + int parent = classes[id].parent;
4.733 + return nodes[parent >= 0 ? classes[id].depth : leftNode(id)].item;
4.734 + }
4.735 +
4.736 + /// \brief Lemon style iterator for the items of a class.
4.737 + ///
4.738 + /// ClassIt is a lemon style iterator for the components. It iterates
4.739 + /// on the items of a class. By example if you want to iterate on
4.740 + /// each items of each classes then you may write the next code.
4.741 + ///\code
4.742 + /// for (ClassIt cit(huf); cit != INVALID; ++cit) {
4.743 + /// std::cout << "Class: ";
4.744 + /// for (ItemIt iit(huf, cit); iit != INVALID; ++iit) {
4.745 + /// std::cout << toString(iit) << ' ' << std::endl;
4.746 + /// }
4.747 + /// std::cout << std::endl;
4.748 + /// }
4.749 + ///\endcode
4.750 + class ItemIt {
4.751 + private:
4.752 +
4.753 + const HeapUnionFind* _huf;
4.754 + int _id, _lid;
4.755 +
4.756 + public:
4.757 +
4.758 + /// \brief Default constructor
4.759 + ///
4.760 + /// Default constructor
4.761 + ItemIt() {}
4.762 +
4.763 + ItemIt(const HeapUnionFind& huf, int cls) : _huf(&huf) {
4.764 + int id = cls;
4.765 + int parent = _huf->classes[id].parent;
4.766 + if (parent >= 0) {
4.767 + _id = _huf->classes[id].depth;
4.768 + if (_huf->classes[id].next != -1) {
4.769 + _lid = _huf->classes[_huf->classes[id].next].depth;
4.770 + } else {
4.771 + _lid = -1;
4.772 + }
4.773 + } else {
4.774 + _id = _huf->leftNode(id);
4.775 + _lid = -1;
4.776 + }
4.777 + }
4.778 +
4.779 + /// \brief Increment operator
4.780 + ///
4.781 + /// It steps to the next item in the class.
4.782 + ItemIt& operator++() {
4.783 + _id = _huf->nextNode(_id);
4.784 + return *this;
4.785 + }
4.786 +
4.787 + /// \brief Conversion operator
4.788 + ///
4.789 + /// It converts the iterator to the current item.
4.790 + operator const Item&() const {
4.791 + return _huf->nodes[_id].item;
4.792 + }
4.793 +
4.794 + /// \brief Equality operator
4.795 + ///
4.796 + /// Equality operator
4.797 + bool operator==(const ItemIt& i) {
4.798 + return i._id == _id;
4.799 + }
4.800 +
4.801 + /// \brief Inequality operator
4.802 + ///
4.803 + /// Inequality operator
4.804 + bool operator!=(const ItemIt& i) {
4.805 + return i._id != _id;
4.806 + }
4.807 +
4.808 + /// \brief Equality operator
4.809 + ///
4.810 + /// Equality operator
4.811 + bool operator==(Invalid) {
4.812 + return _id == _lid;
4.813 + }
4.814 +
4.815 + /// \brief Inequality operator
4.816 + ///
4.817 + /// Inequality operator
4.818 + bool operator!=(Invalid) {
4.819 + return _id != _lid;
4.820 + }
4.821 +
4.822 + };
4.823 +
4.824 + /// \brief Class iterator
4.825 + ///
4.826 + /// The iterator stores
4.827 + class ClassIt {
4.828 + private:
4.829 +
4.830 + const HeapUnionFind* _huf;
4.831 + int _id;
4.832 +
4.833 + public:
4.834 +
4.835 + ClassIt(const HeapUnionFind& huf)
4.836 + : _huf(&huf), _id(huf.first_class) {}
4.837 +
4.838 + ClassIt(const HeapUnionFind& huf, int cls)
4.839 + : _huf(&huf), _id(huf.classes[cls].left) {}
4.840 +
4.841 + ClassIt(Invalid) : _huf(0), _id(-1) {}
4.842 +
4.843 + const ClassIt& operator++() {
4.844 + _id = _huf->classes[_id].next;
4.845 + return *this;
4.846 + }
4.847 +
4.848 + /// \brief Equality operator
4.849 + ///
4.850 + /// Equality operator
4.851 + bool operator==(const ClassIt& i) {
4.852 + return i._id == _id;
4.853 + }
4.854 +
4.855 + /// \brief Inequality operator
4.856 + ///
4.857 + /// Inequality operator
4.858 + bool operator!=(const ClassIt& i) {
4.859 + return i._id != _id;
4.860 + }
4.861 +
4.862 + operator int() const {
4.863 + return _id;
4.864 + }
4.865 +
4.866 + };
4.867 +
4.868 + };
4.869 +
4.870 //! @}
4.871
4.872 } //namespace lemon