1.1 --- a/lemon/max_matching.h Thu Dec 27 13:40:16 2007 +0000
1.2 +++ b/lemon/max_matching.h Fri Dec 28 11:00:51 2007 +0000
1.3 @@ -19,10 +19,14 @@
1.4 #ifndef LEMON_MAX_MATCHING_H
1.5 #define LEMON_MAX_MATCHING_H
1.6
1.7 +#include <vector>
1.8 #include <queue>
1.9 +#include <set>
1.10 +
1.11 #include <lemon/bits/invalid.h>
1.12 #include <lemon/unionfind.h>
1.13 #include <lemon/graph_utils.h>
1.14 +#include <lemon/bin_heap.h>
1.15
1.16 ///\ingroup matching
1.17 ///\file
1.18 @@ -31,7 +35,7 @@
1.19 namespace lemon {
1.20
1.21 ///\ingroup matching
1.22 -
1.23 + ///
1.24 ///\brief Edmonds' alternating forest maximum matching algorithm.
1.25 ///
1.26 ///This class provides Edmonds' alternating forest matching
1.27 @@ -618,6 +622,2500 @@
1.28 }
1.29
1.30 };
1.31 +
1.32 + /// \ingroup matching
1.33 + ///
1.34 + /// \brief Weighted matching in general undirected graphs
1.35 + ///
1.36 + /// This class provides an efficient implementation of Edmond's
1.37 + /// maximum weighted matching algorithm. The implementation is based
1.38 + /// on extensive use of priority queues and provides
1.39 + /// \f$O(nm\log(n))\f$ time complexity.
1.40 + ///
1.41 + /// The maximum weighted matching problem is to find undirected
1.42 + /// edges in the graph with maximum overall weight and no two of
1.43 + /// them shares their endpoints. The problem can be formulated with
1.44 + /// the next linear program:
1.45 + /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
1.46 + ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
1.47 + /// \f[x_e \ge 0\quad \forall e\in E\f]
1.48 + /// \f[\max \sum_{e\in E}x_ew_e\f]
1.49 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1.50 + /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both endpoints in
1.51 + /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
1.52 + /// the nodes.
1.53 + ///
1.54 + /// The algorithm calculates an optimal matching and a proof of the
1.55 + /// optimality. The solution of the dual problem can be used to check
1.56 + /// the result of the algorithm. The dual linear problem is the next:
1.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]
1.58 + /// \f[y_u \ge 0 \quad \forall u \in V\f]
1.59 + /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
1.60 + /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
1.61 + ///
1.62 + /// The algorithm can be executed with \c run() or the \c init() and
1.63 + /// then the \c start() member functions. After it the matching can
1.64 + /// be asked with \c matching() or mate() functions. The dual
1.65 + /// solution can be get with \c nodeValue(), \c blossomNum() and \c
1.66 + /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
1.67 + /// "BlossomIt" nested class which is able to iterate on the nodes
1.68 + /// of a blossom. If the value type is integral then the dual
1.69 + /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
1.70 + ///
1.71 + /// \author Balazs Dezso
1.72 + template <typename _UGraph,
1.73 + typename _WeightMap = typename _UGraph::template UEdgeMap<int> >
1.74 + class MaxWeightedMatching {
1.75 + public:
1.76 +
1.77 + typedef _UGraph UGraph;
1.78 + typedef _WeightMap WeightMap;
1.79 + typedef typename WeightMap::Value Value;
1.80 +
1.81 + /// \brief Scaling factor for dual solution
1.82 + ///
1.83 + /// Scaling factor for dual solution, it is equal to 4 or 1
1.84 + /// according to the value type.
1.85 + static const int dualScale =
1.86 + std::numeric_limits<Value>::is_integer ? 4 : 1;
1.87 +
1.88 + typedef typename UGraph::template NodeMap<typename UGraph::Edge>
1.89 + MatchingMap;
1.90 +
1.91 + private:
1.92 +
1.93 + UGRAPH_TYPEDEFS(typename UGraph);
1.94 +
1.95 + typedef typename UGraph::template NodeMap<Value> NodePotential;
1.96 + typedef std::vector<Node> BlossomNodeList;
1.97 +
1.98 + struct BlossomVariable {
1.99 + int begin, end;
1.100 + Value value;
1.101 +
1.102 + BlossomVariable(int _begin, int _end, Value _value)
1.103 + : begin(_begin), end(_end), value(_value) {}
1.104 +
1.105 + };
1.106 +
1.107 + typedef std::vector<BlossomVariable> BlossomPotential;
1.108 +
1.109 + const UGraph& _ugraph;
1.110 + const WeightMap& _weight;
1.111 +
1.112 + MatchingMap* _matching;
1.113 +
1.114 + NodePotential* _node_potential;
1.115 +
1.116 + BlossomPotential _blossom_potential;
1.117 + BlossomNodeList _blossom_node_list;
1.118 +
1.119 + int _node_num;
1.120 + int _blossom_num;
1.121 +
1.122 + typedef typename UGraph::template NodeMap<int> NodeIntMap;
1.123 + typedef typename UGraph::template EdgeMap<int> EdgeIntMap;
1.124 + typedef typename UGraph::template UEdgeMap<int> UEdgeIntMap;
1.125 + typedef IntegerMap<int> IntIntMap;
1.126 +
1.127 + enum Status {
1.128 + EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
1.129 + };
1.130 +
1.131 + typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
1.132 + struct BlossomData {
1.133 + int tree;
1.134 + Status status;
1.135 + Edge pred, next;
1.136 + Value pot, offset;
1.137 + Node base;
1.138 + };
1.139 +
1.140 + NodeIntMap *_blossom_index;
1.141 + BlossomSet *_blossom_set;
1.142 + IntegerMap<BlossomData>* _blossom_data;
1.143 +
1.144 + NodeIntMap *_node_index;
1.145 + EdgeIntMap *_node_heap_index;
1.146 +
1.147 + struct NodeData {
1.148 +
1.149 + NodeData(EdgeIntMap& node_heap_index)
1.150 + : heap(node_heap_index) {}
1.151 +
1.152 + int blossom;
1.153 + Value pot;
1.154 + BinHeap<Value, EdgeIntMap> heap;
1.155 + std::map<int, Edge> heap_index;
1.156 +
1.157 + int tree;
1.158 + };
1.159 +
1.160 + IntegerMap<NodeData>* _node_data;
1.161 +
1.162 + typedef ExtendFindEnum<IntIntMap> TreeSet;
1.163 +
1.164 + IntIntMap *_tree_set_index;
1.165 + TreeSet *_tree_set;
1.166 +
1.167 + NodeIntMap *_delta1_index;
1.168 + BinHeap<Value, NodeIntMap> *_delta1;
1.169 +
1.170 + IntIntMap *_delta2_index;
1.171 + BinHeap<Value, IntIntMap> *_delta2;
1.172 +
1.173 + UEdgeIntMap *_delta3_index;
1.174 + BinHeap<Value, UEdgeIntMap> *_delta3;
1.175 +
1.176 + IntIntMap *_delta4_index;
1.177 + BinHeap<Value, IntIntMap> *_delta4;
1.178 +
1.179 + Value _delta_sum;
1.180 +
1.181 + void createStructures() {
1.182 + _node_num = countNodes(_ugraph);
1.183 + _blossom_num = _node_num * 3 / 2;
1.184 +
1.185 + if (!_matching) {
1.186 + _matching = new MatchingMap(_ugraph);
1.187 + }
1.188 + if (!_node_potential) {
1.189 + _node_potential = new NodePotential(_ugraph);
1.190 + }
1.191 + if (!_blossom_set) {
1.192 + _blossom_index = new NodeIntMap(_ugraph);
1.193 + _blossom_set = new BlossomSet(*_blossom_index);
1.194 + _blossom_data = new IntegerMap<BlossomData>(_blossom_num);
1.195 + }
1.196 +
1.197 + if (!_node_index) {
1.198 + _node_index = new NodeIntMap(_ugraph);
1.199 + _node_heap_index = new EdgeIntMap(_ugraph);
1.200 + _node_data = new IntegerMap<NodeData>(_node_num,
1.201 + NodeData(*_node_heap_index));
1.202 + }
1.203 +
1.204 + if (!_tree_set) {
1.205 + _tree_set_index = new IntIntMap(_blossom_num);
1.206 + _tree_set = new TreeSet(*_tree_set_index);
1.207 + }
1.208 + if (!_delta1) {
1.209 + _delta1_index = new NodeIntMap(_ugraph);
1.210 + _delta1 = new BinHeap<Value, NodeIntMap>(*_delta1_index);
1.211 + }
1.212 + if (!_delta2) {
1.213 + _delta2_index = new IntIntMap(_blossom_num);
1.214 + _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
1.215 + }
1.216 + if (!_delta3) {
1.217 + _delta3_index = new UEdgeIntMap(_ugraph);
1.218 + _delta3 = new BinHeap<Value, UEdgeIntMap>(*_delta3_index);
1.219 + }
1.220 + if (!_delta4) {
1.221 + _delta4_index = new IntIntMap(_blossom_num);
1.222 + _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
1.223 + }
1.224 + }
1.225 +
1.226 + void destroyStructures() {
1.227 + _node_num = countNodes(_ugraph);
1.228 + _blossom_num = _node_num * 3 / 2;
1.229 +
1.230 + if (_matching) {
1.231 + delete _matching;
1.232 + }
1.233 + if (_node_potential) {
1.234 + delete _node_potential;
1.235 + }
1.236 + if (_blossom_set) {
1.237 + delete _blossom_index;
1.238 + delete _blossom_set;
1.239 + delete _blossom_data;
1.240 + }
1.241 +
1.242 + if (_node_index) {
1.243 + delete _node_index;
1.244 + delete _node_heap_index;
1.245 + delete _node_data;
1.246 + }
1.247 +
1.248 + if (_tree_set) {
1.249 + delete _tree_set_index;
1.250 + delete _tree_set;
1.251 + }
1.252 + if (_delta1) {
1.253 + delete _delta1_index;
1.254 + delete _delta1;
1.255 + }
1.256 + if (_delta2) {
1.257 + delete _delta2_index;
1.258 + delete _delta2;
1.259 + }
1.260 + if (_delta3) {
1.261 + delete _delta3_index;
1.262 + delete _delta3;
1.263 + }
1.264 + if (_delta4) {
1.265 + delete _delta4_index;
1.266 + delete _delta4;
1.267 + }
1.268 + }
1.269 +
1.270 + void matchedToEven(int blossom, int tree) {
1.271 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.272 + _delta2->erase(blossom);
1.273 + }
1.274 +
1.275 + if (!_blossom_set->trivial(blossom)) {
1.276 + (*_blossom_data)[blossom].pot -=
1.277 + 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
1.278 + }
1.279 +
1.280 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.281 + n != INVALID; ++n) {
1.282 +
1.283 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
1.284 + int ni = (*_node_index)[n];
1.285 +
1.286 + (*_node_data)[ni].heap.clear();
1.287 + (*_node_data)[ni].heap_index.clear();
1.288 +
1.289 + (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
1.290 +
1.291 + _delta1->push(n, (*_node_data)[ni].pot);
1.292 +
1.293 + for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
1.294 + Node v = _ugraph.source(e);
1.295 + int vb = _blossom_set->find(v);
1.296 + int vi = (*_node_index)[v];
1.297 +
1.298 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.299 + dualScale * _weight[e];
1.300 +
1.301 + if ((*_blossom_data)[vb].status == EVEN) {
1.302 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1.303 + _delta3->push(e, rw / 2);
1.304 + }
1.305 + } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1.306 + if (_delta3->state(e) != _delta3->IN_HEAP) {
1.307 + _delta3->push(e, rw);
1.308 + }
1.309 + } else {
1.310 + typename std::map<int, Edge>::iterator it =
1.311 + (*_node_data)[vi].heap_index.find(tree);
1.312 +
1.313 + if (it != (*_node_data)[vi].heap_index.end()) {
1.314 + if ((*_node_data)[vi].heap[it->second] > rw) {
1.315 + (*_node_data)[vi].heap.replace(it->second, e);
1.316 + (*_node_data)[vi].heap.decrease(e, rw);
1.317 + it->second = e;
1.318 + }
1.319 + } else {
1.320 + (*_node_data)[vi].heap.push(e, rw);
1.321 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1.322 + }
1.323 +
1.324 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1.325 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1.326 +
1.327 + if ((*_blossom_data)[vb].status == MATCHED) {
1.328 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
1.329 + _delta2->push(vb, _blossom_set->classPrio(vb) -
1.330 + (*_blossom_data)[vb].offset);
1.331 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1.332 + (*_blossom_data)[vb].offset){
1.333 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1.334 + (*_blossom_data)[vb].offset);
1.335 + }
1.336 + }
1.337 + }
1.338 + }
1.339 + }
1.340 + }
1.341 + (*_blossom_data)[blossom].offset = 0;
1.342 + }
1.343 +
1.344 + void matchedToOdd(int blossom) {
1.345 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.346 + _delta2->erase(blossom);
1.347 + }
1.348 + (*_blossom_data)[blossom].offset += _delta_sum;
1.349 + if (!_blossom_set->trivial(blossom)) {
1.350 + _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
1.351 + (*_blossom_data)[blossom].offset);
1.352 + }
1.353 + }
1.354 +
1.355 + void evenToMatched(int blossom, int tree) {
1.356 + if (!_blossom_set->trivial(blossom)) {
1.357 + (*_blossom_data)[blossom].pot += 2 * _delta_sum;
1.358 + }
1.359 +
1.360 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.361 + n != INVALID; ++n) {
1.362 + int ni = (*_node_index)[n];
1.363 + (*_node_data)[ni].pot -= _delta_sum;
1.364 +
1.365 + _delta1->erase(n);
1.366 +
1.367 + for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
1.368 + Node v = _ugraph.source(e);
1.369 + int vb = _blossom_set->find(v);
1.370 + int vi = (*_node_index)[v];
1.371 +
1.372 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.373 + dualScale * _weight[e];
1.374 +
1.375 + if (vb == blossom) {
1.376 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.377 + _delta3->erase(e);
1.378 + }
1.379 + } else if ((*_blossom_data)[vb].status == EVEN) {
1.380 +
1.381 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.382 + _delta3->erase(e);
1.383 + }
1.384 +
1.385 + int vt = _tree_set->find(vb);
1.386 +
1.387 + if (vt != tree) {
1.388 +
1.389 + Edge r = _ugraph.oppositeEdge(e);
1.390 +
1.391 + typename std::map<int, Edge>::iterator it =
1.392 + (*_node_data)[ni].heap_index.find(vt);
1.393 +
1.394 + if (it != (*_node_data)[ni].heap_index.end()) {
1.395 + if ((*_node_data)[ni].heap[it->second] > rw) {
1.396 + (*_node_data)[ni].heap.replace(it->second, r);
1.397 + (*_node_data)[ni].heap.decrease(r, rw);
1.398 + it->second = r;
1.399 + }
1.400 + } else {
1.401 + (*_node_data)[ni].heap.push(r, rw);
1.402 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1.403 + }
1.404 +
1.405 + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1.406 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1.407 +
1.408 + if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1.409 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.410 + (*_blossom_data)[blossom].offset);
1.411 + } else if ((*_delta2)[blossom] >
1.412 + _blossom_set->classPrio(blossom) -
1.413 + (*_blossom_data)[blossom].offset){
1.414 + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1.415 + (*_blossom_data)[blossom].offset);
1.416 + }
1.417 + }
1.418 + }
1.419 +
1.420 + } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1.421 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.422 + _delta3->erase(e);
1.423 + }
1.424 + } else {
1.425 +
1.426 + typename std::map<int, Edge>::iterator it =
1.427 + (*_node_data)[vi].heap_index.find(tree);
1.428 +
1.429 + if (it != (*_node_data)[vi].heap_index.end()) {
1.430 + (*_node_data)[vi].heap.erase(it->second);
1.431 + (*_node_data)[vi].heap_index.erase(it);
1.432 + if ((*_node_data)[vi].heap.empty()) {
1.433 + _blossom_set->increase(v, std::numeric_limits<Value>::max());
1.434 + } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1.435 + _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1.436 + }
1.437 +
1.438 + if ((*_blossom_data)[vb].status == MATCHED) {
1.439 + if (_blossom_set->classPrio(vb) ==
1.440 + std::numeric_limits<Value>::max()) {
1.441 + _delta2->erase(vb);
1.442 + } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1.443 + (*_blossom_data)[vb].offset) {
1.444 + _delta2->increase(vb, _blossom_set->classPrio(vb) -
1.445 + (*_blossom_data)[vb].offset);
1.446 + }
1.447 + }
1.448 + }
1.449 + }
1.450 + }
1.451 + }
1.452 + }
1.453 +
1.454 + void oddToMatched(int blossom) {
1.455 + (*_blossom_data)[blossom].offset -= _delta_sum;
1.456 +
1.457 + if (_blossom_set->classPrio(blossom) !=
1.458 + std::numeric_limits<Value>::max()) {
1.459 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.460 + (*_blossom_data)[blossom].offset);
1.461 + }
1.462 +
1.463 + if (!_blossom_set->trivial(blossom)) {
1.464 + _delta4->erase(blossom);
1.465 + }
1.466 + }
1.467 +
1.468 + void oddToEven(int blossom, int tree) {
1.469 + if (!_blossom_set->trivial(blossom)) {
1.470 + _delta4->erase(blossom);
1.471 + (*_blossom_data)[blossom].pot -=
1.472 + 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1.473 + }
1.474 +
1.475 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.476 + n != INVALID; ++n) {
1.477 + int ni = (*_node_index)[n];
1.478 +
1.479 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
1.480 +
1.481 + (*_node_data)[ni].heap.clear();
1.482 + (*_node_data)[ni].heap_index.clear();
1.483 + (*_node_data)[ni].pot +=
1.484 + 2 * _delta_sum - (*_blossom_data)[blossom].offset;
1.485 +
1.486 + _delta1->push(n, (*_node_data)[ni].pot);
1.487 +
1.488 + for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
1.489 + Node v = _ugraph.source(e);
1.490 + int vb = _blossom_set->find(v);
1.491 + int vi = (*_node_index)[v];
1.492 +
1.493 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.494 + dualScale * _weight[e];
1.495 +
1.496 + if ((*_blossom_data)[vb].status == EVEN) {
1.497 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1.498 + _delta3->push(e, rw / 2);
1.499 + }
1.500 + } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1.501 + if (_delta3->state(e) != _delta3->IN_HEAP) {
1.502 + _delta3->push(e, rw);
1.503 + }
1.504 + } else {
1.505 +
1.506 + typename std::map<int, Edge>::iterator it =
1.507 + (*_node_data)[vi].heap_index.find(tree);
1.508 +
1.509 + if (it != (*_node_data)[vi].heap_index.end()) {
1.510 + if ((*_node_data)[vi].heap[it->second] > rw) {
1.511 + (*_node_data)[vi].heap.replace(it->second, e);
1.512 + (*_node_data)[vi].heap.decrease(e, rw);
1.513 + it->second = e;
1.514 + }
1.515 + } else {
1.516 + (*_node_data)[vi].heap.push(e, rw);
1.517 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1.518 + }
1.519 +
1.520 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1.521 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1.522 +
1.523 + if ((*_blossom_data)[vb].status == MATCHED) {
1.524 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
1.525 + _delta2->push(vb, _blossom_set->classPrio(vb) -
1.526 + (*_blossom_data)[vb].offset);
1.527 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1.528 + (*_blossom_data)[vb].offset) {
1.529 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1.530 + (*_blossom_data)[vb].offset);
1.531 + }
1.532 + }
1.533 + }
1.534 + }
1.535 + }
1.536 + }
1.537 + (*_blossom_data)[blossom].offset = 0;
1.538 + }
1.539 +
1.540 +
1.541 + void matchedToUnmatched(int blossom) {
1.542 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.543 + _delta2->erase(blossom);
1.544 + }
1.545 +
1.546 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.547 + n != INVALID; ++n) {
1.548 + int ni = (*_node_index)[n];
1.549 +
1.550 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
1.551 +
1.552 + (*_node_data)[ni].heap.clear();
1.553 + (*_node_data)[ni].heap_index.clear();
1.554 +
1.555 + for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
1.556 + Node v = _ugraph.target(e);
1.557 + int vb = _blossom_set->find(v);
1.558 + int vi = (*_node_index)[v];
1.559 +
1.560 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.561 + dualScale * _weight[e];
1.562 +
1.563 + if ((*_blossom_data)[vb].status == EVEN) {
1.564 + if (_delta3->state(e) != _delta3->IN_HEAP) {
1.565 + _delta3->push(e, rw);
1.566 + }
1.567 + }
1.568 + }
1.569 + }
1.570 + }
1.571 +
1.572 + void unmatchedToMatched(int blossom) {
1.573 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.574 + n != INVALID; ++n) {
1.575 + int ni = (*_node_index)[n];
1.576 +
1.577 + for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
1.578 + Node v = _ugraph.source(e);
1.579 + int vb = _blossom_set->find(v);
1.580 + int vi = (*_node_index)[v];
1.581 +
1.582 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.583 + dualScale * _weight[e];
1.584 +
1.585 + if (vb == blossom) {
1.586 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.587 + _delta3->erase(e);
1.588 + }
1.589 + } else if ((*_blossom_data)[vb].status == EVEN) {
1.590 +
1.591 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.592 + _delta3->erase(e);
1.593 + }
1.594 +
1.595 + int vt = _tree_set->find(vb);
1.596 +
1.597 + Edge r = _ugraph.oppositeEdge(e);
1.598 +
1.599 + typename std::map<int, Edge>::iterator it =
1.600 + (*_node_data)[ni].heap_index.find(vt);
1.601 +
1.602 + if (it != (*_node_data)[ni].heap_index.end()) {
1.603 + if ((*_node_data)[ni].heap[it->second] > rw) {
1.604 + (*_node_data)[ni].heap.replace(it->second, r);
1.605 + (*_node_data)[ni].heap.decrease(r, rw);
1.606 + it->second = r;
1.607 + }
1.608 + } else {
1.609 + (*_node_data)[ni].heap.push(r, rw);
1.610 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1.611 + }
1.612 +
1.613 + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1.614 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1.615 +
1.616 + if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1.617 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.618 + (*_blossom_data)[blossom].offset);
1.619 + } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
1.620 + (*_blossom_data)[blossom].offset){
1.621 + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1.622 + (*_blossom_data)[blossom].offset);
1.623 + }
1.624 + }
1.625 +
1.626 + } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1.627 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.628 + _delta3->erase(e);
1.629 + }
1.630 + }
1.631 + }
1.632 + }
1.633 + }
1.634 +
1.635 + void alternatePath(int even, int tree) {
1.636 + int odd;
1.637 +
1.638 + evenToMatched(even, tree);
1.639 + (*_blossom_data)[even].status = MATCHED;
1.640 +
1.641 + while ((*_blossom_data)[even].pred != INVALID) {
1.642 + odd = _blossom_set->find(_ugraph.target((*_blossom_data)[even].pred));
1.643 + (*_blossom_data)[odd].status = MATCHED;
1.644 + oddToMatched(odd);
1.645 + (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
1.646 +
1.647 + even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].pred));
1.648 + (*_blossom_data)[even].status = MATCHED;
1.649 + evenToMatched(even, tree);
1.650 + (*_blossom_data)[even].next =
1.651 + _ugraph.oppositeEdge((*_blossom_data)[odd].pred);
1.652 + }
1.653 +
1.654 + }
1.655 +
1.656 + void destroyTree(int tree) {
1.657 + for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
1.658 + if ((*_blossom_data)[b].status == EVEN) {
1.659 + (*_blossom_data)[b].status = MATCHED;
1.660 + evenToMatched(b, tree);
1.661 + } else if ((*_blossom_data)[b].status == ODD) {
1.662 + (*_blossom_data)[b].status = MATCHED;
1.663 + oddToMatched(b);
1.664 + }
1.665 + }
1.666 + _tree_set->eraseClass(tree);
1.667 + }
1.668 +
1.669 +
1.670 + void unmatchNode(const Node& node) {
1.671 + int blossom = _blossom_set->find(node);
1.672 + int tree = _tree_set->find(blossom);
1.673 +
1.674 + alternatePath(blossom, tree);
1.675 + destroyTree(tree);
1.676 +
1.677 + (*_blossom_data)[blossom].status = UNMATCHED;
1.678 + (*_blossom_data)[blossom].base = node;
1.679 + matchedToUnmatched(blossom);
1.680 + }
1.681 +
1.682 +
1.683 + void augmentOnEdge(const UEdge& edge) {
1.684 +
1.685 + int left = _blossom_set->find(_ugraph.source(edge));
1.686 + int right = _blossom_set->find(_ugraph.target(edge));
1.687 +
1.688 + if ((*_blossom_data)[left].status == EVEN) {
1.689 + int left_tree = _tree_set->find(left);
1.690 + alternatePath(left, left_tree);
1.691 + destroyTree(left_tree);
1.692 + } else {
1.693 + (*_blossom_data)[left].status = MATCHED;
1.694 + unmatchedToMatched(left);
1.695 + }
1.696 +
1.697 + if ((*_blossom_data)[right].status == EVEN) {
1.698 + int right_tree = _tree_set->find(right);
1.699 + alternatePath(right, right_tree);
1.700 + destroyTree(right_tree);
1.701 + } else {
1.702 + (*_blossom_data)[right].status = MATCHED;
1.703 + unmatchedToMatched(right);
1.704 + }
1.705 +
1.706 + (*_blossom_data)[left].next = _ugraph.direct(edge, true);
1.707 + (*_blossom_data)[right].next = _ugraph.direct(edge, false);
1.708 + }
1.709 +
1.710 + void extendOnEdge(const Edge& edge) {
1.711 + int base = _blossom_set->find(_ugraph.target(edge));
1.712 + int tree = _tree_set->find(base);
1.713 +
1.714 + int odd = _blossom_set->find(_ugraph.source(edge));
1.715 + _tree_set->insert(odd, tree);
1.716 + (*_blossom_data)[odd].status = ODD;
1.717 + matchedToOdd(odd);
1.718 + (*_blossom_data)[odd].pred = edge;
1.719 +
1.720 + int even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].next));
1.721 + (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
1.722 + _tree_set->insert(even, tree);
1.723 + (*_blossom_data)[even].status = EVEN;
1.724 + matchedToEven(even, tree);
1.725 + }
1.726 +
1.727 + void shrinkOnEdge(const UEdge& uedge, int tree) {
1.728 + int nca = -1;
1.729 + std::vector<int> left_path, right_path;
1.730 +
1.731 + {
1.732 + std::set<int> left_set, right_set;
1.733 + int left = _blossom_set->find(_ugraph.source(uedge));
1.734 + left_path.push_back(left);
1.735 + left_set.insert(left);
1.736 +
1.737 + int right = _blossom_set->find(_ugraph.target(uedge));
1.738 + right_path.push_back(right);
1.739 + right_set.insert(right);
1.740 +
1.741 + while (true) {
1.742 +
1.743 + if ((*_blossom_data)[left].pred == INVALID) break;
1.744 +
1.745 + left =
1.746 + _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
1.747 + left_path.push_back(left);
1.748 + left =
1.749 + _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
1.750 + left_path.push_back(left);
1.751 +
1.752 + left_set.insert(left);
1.753 +
1.754 + if (right_set.find(left) != right_set.end()) {
1.755 + nca = left;
1.756 + break;
1.757 + }
1.758 +
1.759 + if ((*_blossom_data)[right].pred == INVALID) break;
1.760 +
1.761 + right =
1.762 + _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
1.763 + right_path.push_back(right);
1.764 + right =
1.765 + _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
1.766 + right_path.push_back(right);
1.767 +
1.768 + right_set.insert(right);
1.769 +
1.770 + if (left_set.find(right) != left_set.end()) {
1.771 + nca = right;
1.772 + break;
1.773 + }
1.774 +
1.775 + }
1.776 +
1.777 + if (nca == -1) {
1.778 + if ((*_blossom_data)[left].pred == INVALID) {
1.779 + nca = right;
1.780 + while (left_set.find(nca) == left_set.end()) {
1.781 + nca =
1.782 + _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
1.783 + right_path.push_back(nca);
1.784 + nca =
1.785 + _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
1.786 + right_path.push_back(nca);
1.787 + }
1.788 + } else {
1.789 + nca = left;
1.790 + while (right_set.find(nca) == right_set.end()) {
1.791 + nca =
1.792 + _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
1.793 + left_path.push_back(nca);
1.794 + nca =
1.795 + _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
1.796 + left_path.push_back(nca);
1.797 + }
1.798 + }
1.799 + }
1.800 + }
1.801 +
1.802 + std::vector<int> subblossoms;
1.803 + Edge prev;
1.804 +
1.805 + prev = _ugraph.direct(uedge, true);
1.806 + for (int i = 0; left_path[i] != nca; i += 2) {
1.807 + subblossoms.push_back(left_path[i]);
1.808 + (*_blossom_data)[left_path[i]].next = prev;
1.809 + _tree_set->erase(left_path[i]);
1.810 +
1.811 + subblossoms.push_back(left_path[i + 1]);
1.812 + (*_blossom_data)[left_path[i + 1]].status = EVEN;
1.813 + oddToEven(left_path[i + 1], tree);
1.814 + _tree_set->erase(left_path[i + 1]);
1.815 + prev = _ugraph.oppositeEdge((*_blossom_data)[left_path[i + 1]].pred);
1.816 + }
1.817 +
1.818 + int k = 0;
1.819 + while (right_path[k] != nca) ++k;
1.820 +
1.821 + subblossoms.push_back(nca);
1.822 + (*_blossom_data)[nca].next = prev;
1.823 +
1.824 + for (int i = k - 2; i >= 0; i -= 2) {
1.825 + subblossoms.push_back(right_path[i + 1]);
1.826 + (*_blossom_data)[right_path[i + 1]].status = EVEN;
1.827 + oddToEven(right_path[i + 1], tree);
1.828 + _tree_set->erase(right_path[i + 1]);
1.829 +
1.830 + (*_blossom_data)[right_path[i + 1]].next =
1.831 + (*_blossom_data)[right_path[i + 1]].pred;
1.832 +
1.833 + subblossoms.push_back(right_path[i]);
1.834 + _tree_set->erase(right_path[i]);
1.835 + }
1.836 +
1.837 + int surface =
1.838 + _blossom_set->join(subblossoms.begin(), subblossoms.end());
1.839 +
1.840 + for (int i = 0; i < int(subblossoms.size()); ++i) {
1.841 + if (!_blossom_set->trivial(subblossoms[i])) {
1.842 + (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
1.843 + }
1.844 + (*_blossom_data)[subblossoms[i]].status = MATCHED;
1.845 + }
1.846 +
1.847 + (*_blossom_data)[surface].pot = -2 * _delta_sum;
1.848 + (*_blossom_data)[surface].offset = 0;
1.849 + (*_blossom_data)[surface].status = EVEN;
1.850 + (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
1.851 + (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
1.852 +
1.853 + _tree_set->insert(surface, tree);
1.854 + _tree_set->erase(nca);
1.855 + }
1.856 +
1.857 + void splitBlossom(int blossom) {
1.858 + Edge next = (*_blossom_data)[blossom].next;
1.859 + Edge pred = (*_blossom_data)[blossom].pred;
1.860 +
1.861 + int tree = _tree_set->find(blossom);
1.862 +
1.863 + (*_blossom_data)[blossom].status = MATCHED;
1.864 + oddToMatched(blossom);
1.865 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.866 + _delta2->erase(blossom);
1.867 + }
1.868 +
1.869 + std::vector<int> subblossoms;
1.870 + _blossom_set->split(blossom, std::back_inserter(subblossoms));
1.871 +
1.872 + Value offset = (*_blossom_data)[blossom].offset;
1.873 + int b = _blossom_set->find(_ugraph.source(pred));
1.874 + int d = _blossom_set->find(_ugraph.source(next));
1.875 +
1.876 + int ib, id;
1.877 + for (int i = 0; i < int(subblossoms.size()); ++i) {
1.878 + if (subblossoms[i] == b) ib = i;
1.879 + if (subblossoms[i] == d) id = i;
1.880 +
1.881 + (*_blossom_data)[subblossoms[i]].offset = offset;
1.882 + if (!_blossom_set->trivial(subblossoms[i])) {
1.883 + (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
1.884 + }
1.885 + if (_blossom_set->classPrio(subblossoms[i]) !=
1.886 + std::numeric_limits<Value>::max()) {
1.887 + _delta2->push(subblossoms[i],
1.888 + _blossom_set->classPrio(subblossoms[i]) -
1.889 + (*_blossom_data)[subblossoms[i]].offset);
1.890 + }
1.891 + }
1.892 +
1.893 + if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
1.894 + for (int i = (id + 1) % subblossoms.size();
1.895 + i != ib; i = (i + 2) % subblossoms.size()) {
1.896 + int sb = subblossoms[i];
1.897 + int tb = subblossoms[(i + 1) % subblossoms.size()];
1.898 + (*_blossom_data)[sb].next =
1.899 + _ugraph.oppositeEdge((*_blossom_data)[tb].next);
1.900 + }
1.901 +
1.902 + for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
1.903 + int sb = subblossoms[i];
1.904 + int tb = subblossoms[(i + 1) % subblossoms.size()];
1.905 + int ub = subblossoms[(i + 2) % subblossoms.size()];
1.906 +
1.907 + (*_blossom_data)[sb].status = ODD;
1.908 + matchedToOdd(sb);
1.909 + _tree_set->insert(sb, tree);
1.910 + (*_blossom_data)[sb].pred = pred;
1.911 + (*_blossom_data)[sb].next =
1.912 + _ugraph.oppositeEdge((*_blossom_data)[tb].next);
1.913 +
1.914 + pred = (*_blossom_data)[ub].next;
1.915 +
1.916 + (*_blossom_data)[tb].status = EVEN;
1.917 + matchedToEven(tb, tree);
1.918 + _tree_set->insert(tb, tree);
1.919 + (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
1.920 + }
1.921 +
1.922 + (*_blossom_data)[subblossoms[id]].status = ODD;
1.923 + matchedToOdd(subblossoms[id]);
1.924 + _tree_set->insert(subblossoms[id], tree);
1.925 + (*_blossom_data)[subblossoms[id]].next = next;
1.926 + (*_blossom_data)[subblossoms[id]].pred = pred;
1.927 +
1.928 + } else {
1.929 +
1.930 + for (int i = (ib + 1) % subblossoms.size();
1.931 + i != id; i = (i + 2) % subblossoms.size()) {
1.932 + int sb = subblossoms[i];
1.933 + int tb = subblossoms[(i + 1) % subblossoms.size()];
1.934 + (*_blossom_data)[sb].next =
1.935 + _ugraph.oppositeEdge((*_blossom_data)[tb].next);
1.936 + }
1.937 +
1.938 + for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
1.939 + int sb = subblossoms[i];
1.940 + int tb = subblossoms[(i + 1) % subblossoms.size()];
1.941 + int ub = subblossoms[(i + 2) % subblossoms.size()];
1.942 +
1.943 + (*_blossom_data)[sb].status = ODD;
1.944 + matchedToOdd(sb);
1.945 + _tree_set->insert(sb, tree);
1.946 + (*_blossom_data)[sb].next = next;
1.947 + (*_blossom_data)[sb].pred =
1.948 + _ugraph.oppositeEdge((*_blossom_data)[tb].next);
1.949 +
1.950 + (*_blossom_data)[tb].status = EVEN;
1.951 + matchedToEven(tb, tree);
1.952 + _tree_set->insert(tb, tree);
1.953 + (*_blossom_data)[tb].pred =
1.954 + (*_blossom_data)[tb].next =
1.955 + _ugraph.oppositeEdge((*_blossom_data)[ub].next);
1.956 + next = (*_blossom_data)[ub].next;
1.957 + }
1.958 +
1.959 + (*_blossom_data)[subblossoms[ib]].status = ODD;
1.960 + matchedToOdd(subblossoms[ib]);
1.961 + _tree_set->insert(subblossoms[ib], tree);
1.962 + (*_blossom_data)[subblossoms[ib]].next = next;
1.963 + (*_blossom_data)[subblossoms[ib]].pred = pred;
1.964 + }
1.965 + _tree_set->erase(blossom);
1.966 + }
1.967 +
1.968 + void extractBlossom(int blossom, const Node& base, const Edge& matching) {
1.969 + if (_blossom_set->trivial(blossom)) {
1.970 + int bi = (*_node_index)[base];
1.971 + Value pot = (*_node_data)[bi].pot;
1.972 +
1.973 + _matching->set(base, matching);
1.974 + _blossom_node_list.push_back(base);
1.975 + _node_potential->set(base, pot);
1.976 + } else {
1.977 +
1.978 + Value pot = (*_blossom_data)[blossom].pot;
1.979 + int bn = _blossom_node_list.size();
1.980 +
1.981 + std::vector<int> subblossoms;
1.982 + _blossom_set->split(blossom, std::back_inserter(subblossoms));
1.983 + int b = _blossom_set->find(base);
1.984 + int ib = -1;
1.985 + for (int i = 0; i < int(subblossoms.size()); ++i) {
1.986 + if (subblossoms[i] == b) { ib = i; break; }
1.987 + }
1.988 +
1.989 + for (int i = 1; i < int(subblossoms.size()); i += 2) {
1.990 + int sb = subblossoms[(ib + i) % subblossoms.size()];
1.991 + int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
1.992 +
1.993 + Edge m = (*_blossom_data)[tb].next;
1.994 + extractBlossom(sb, _ugraph.target(m), _ugraph.oppositeEdge(m));
1.995 + extractBlossom(tb, _ugraph.source(m), m);
1.996 + }
1.997 + extractBlossom(subblossoms[ib], base, matching);
1.998 +
1.999 + int en = _blossom_node_list.size();
1.1000 +
1.1001 + _blossom_potential.push_back(BlossomVariable(bn, en, pot));
1.1002 + }
1.1003 + }
1.1004 +
1.1005 + void extractMatching() {
1.1006 + std::vector<int> blossoms;
1.1007 + for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1.1008 + blossoms.push_back(c);
1.1009 + }
1.1010 +
1.1011 + for (int i = 0; i < int(blossoms.size()); ++i) {
1.1012 + if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
1.1013 +
1.1014 + Value offset = (*_blossom_data)[blossoms[i]].offset;
1.1015 + (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1.1016 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1.1017 + n != INVALID; ++n) {
1.1018 + (*_node_data)[(*_node_index)[n]].pot -= offset;
1.1019 + }
1.1020 +
1.1021 + Edge matching = (*_blossom_data)[blossoms[i]].next;
1.1022 + Node base = _ugraph.source(matching);
1.1023 + extractBlossom(blossoms[i], base, matching);
1.1024 + } else {
1.1025 + Node base = (*_blossom_data)[blossoms[i]].base;
1.1026 + extractBlossom(blossoms[i], base, INVALID);
1.1027 + }
1.1028 + }
1.1029 + }
1.1030 +
1.1031 + public:
1.1032 +
1.1033 + /// \brief Constructor
1.1034 + ///
1.1035 + /// Constructor.
1.1036 + MaxWeightedMatching(const UGraph& ugraph, const WeightMap& weight)
1.1037 + : _ugraph(ugraph), _weight(weight), _matching(0),
1.1038 + _node_potential(0), _blossom_potential(), _blossom_node_list(),
1.1039 + _node_num(0), _blossom_num(0),
1.1040 +
1.1041 + _blossom_index(0), _blossom_set(0), _blossom_data(0),
1.1042 + _node_index(0), _node_heap_index(0), _node_data(0),
1.1043 + _tree_set_index(0), _tree_set(0),
1.1044 +
1.1045 + _delta1_index(0), _delta1(0),
1.1046 + _delta2_index(0), _delta2(0),
1.1047 + _delta3_index(0), _delta3(0),
1.1048 + _delta4_index(0), _delta4(0),
1.1049 +
1.1050 + _delta_sum() {}
1.1051 +
1.1052 + ~MaxWeightedMatching() {
1.1053 + destroyStructures();
1.1054 + }
1.1055 +
1.1056 + /// \name Execution control
1.1057 + /// The simplest way to execute the algorithm is to use the member
1.1058 + /// \c run() member function.
1.1059 +
1.1060 + ///@{
1.1061 +
1.1062 + /// \brief Initialize the algorithm
1.1063 + ///
1.1064 + /// Initialize the algorithm
1.1065 + void init() {
1.1066 + createStructures();
1.1067 +
1.1068 + for (EdgeIt e(_ugraph); e != INVALID; ++e) {
1.1069 + _node_heap_index->set(e, BinHeap<Value, EdgeIntMap>::PRE_HEAP);
1.1070 + }
1.1071 + for (NodeIt n(_ugraph); n != INVALID; ++n) {
1.1072 + _delta1_index->set(n, _delta1->PRE_HEAP);
1.1073 + }
1.1074 + for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
1.1075 + _delta3_index->set(e, _delta3->PRE_HEAP);
1.1076 + }
1.1077 + for (int i = 0; i < _blossom_num; ++i) {
1.1078 + _delta2_index->set(i, _delta2->PRE_HEAP);
1.1079 + _delta4_index->set(i, _delta4->PRE_HEAP);
1.1080 + }
1.1081 +
1.1082 + int index = 0;
1.1083 + for (NodeIt n(_ugraph); n != INVALID; ++n) {
1.1084 + Value max = 0;
1.1085 + for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
1.1086 + if (_ugraph.target(e) == n) continue;
1.1087 + if ((dualScale * _weight[e]) / 2 > max) {
1.1088 + max = (dualScale * _weight[e]) / 2;
1.1089 + }
1.1090 + }
1.1091 + _node_index->set(n, index);
1.1092 + (*_node_data)[index].pot = max;
1.1093 + _delta1->push(n, max);
1.1094 + int blossom =
1.1095 + _blossom_set->insert(n, std::numeric_limits<Value>::max());
1.1096 +
1.1097 + _tree_set->insert(blossom);
1.1098 +
1.1099 + (*_blossom_data)[blossom].status = EVEN;
1.1100 + (*_blossom_data)[blossom].pred = INVALID;
1.1101 + (*_blossom_data)[blossom].next = INVALID;
1.1102 + (*_blossom_data)[blossom].pot = 0;
1.1103 + (*_blossom_data)[blossom].offset = 0;
1.1104 + ++index;
1.1105 + }
1.1106 + for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
1.1107 + int si = (*_node_index)[_ugraph.source(e)];
1.1108 + int ti = (*_node_index)[_ugraph.target(e)];
1.1109 + if (_ugraph.source(e) != _ugraph.target(e)) {
1.1110 + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1.1111 + dualScale * _weight[e]) / 2);
1.1112 + }
1.1113 + }
1.1114 + }
1.1115 +
1.1116 + /// \brief Starts the algorithm
1.1117 + ///
1.1118 + /// Starts the algorithm
1.1119 + void start() {
1.1120 + enum OpType {
1.1121 + D1, D2, D3, D4
1.1122 + };
1.1123 +
1.1124 + int unmatched = _node_num;
1.1125 + while (unmatched > 0) {
1.1126 + Value d1 = !_delta1->empty() ?
1.1127 + _delta1->prio() : std::numeric_limits<Value>::max();
1.1128 +
1.1129 + Value d2 = !_delta2->empty() ?
1.1130 + _delta2->prio() : std::numeric_limits<Value>::max();
1.1131 +
1.1132 + Value d3 = !_delta3->empty() ?
1.1133 + _delta3->prio() : std::numeric_limits<Value>::max();
1.1134 +
1.1135 + Value d4 = !_delta4->empty() ?
1.1136 + _delta4->prio() : std::numeric_limits<Value>::max();
1.1137 +
1.1138 + _delta_sum = d1; OpType ot = D1;
1.1139 + if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1.1140 + if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
1.1141 + if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1.1142 +
1.1143 +
1.1144 + switch (ot) {
1.1145 + case D1:
1.1146 + {
1.1147 + Node n = _delta1->top();
1.1148 + unmatchNode(n);
1.1149 + --unmatched;
1.1150 + }
1.1151 + break;
1.1152 + case D2:
1.1153 + {
1.1154 + int blossom = _delta2->top();
1.1155 + Node n = _blossom_set->classTop(blossom);
1.1156 + Edge e = (*_node_data)[(*_node_index)[n]].heap.top();
1.1157 + extendOnEdge(e);
1.1158 + }
1.1159 + break;
1.1160 + case D3:
1.1161 + {
1.1162 + UEdge e = _delta3->top();
1.1163 +
1.1164 + int left_blossom = _blossom_set->find(_ugraph.source(e));
1.1165 + int right_blossom = _blossom_set->find(_ugraph.target(e));
1.1166 +
1.1167 + if (left_blossom == right_blossom) {
1.1168 + _delta3->pop();
1.1169 + } else {
1.1170 + int left_tree;
1.1171 + if ((*_blossom_data)[left_blossom].status == EVEN) {
1.1172 + left_tree = _tree_set->find(left_blossom);
1.1173 + } else {
1.1174 + left_tree = -1;
1.1175 + ++unmatched;
1.1176 + }
1.1177 + int right_tree;
1.1178 + if ((*_blossom_data)[right_blossom].status == EVEN) {
1.1179 + right_tree = _tree_set->find(right_blossom);
1.1180 + } else {
1.1181 + right_tree = -1;
1.1182 + ++unmatched;
1.1183 + }
1.1184 +
1.1185 + if (left_tree == right_tree) {
1.1186 + shrinkOnEdge(e, left_tree);
1.1187 + } else {
1.1188 + augmentOnEdge(e);
1.1189 + unmatched -= 2;
1.1190 + }
1.1191 + }
1.1192 + } break;
1.1193 + case D4:
1.1194 + splitBlossom(_delta4->top());
1.1195 + break;
1.1196 + }
1.1197 + }
1.1198 + extractMatching();
1.1199 + }
1.1200 +
1.1201 + /// \brief Runs %MaxWeightedMatching algorithm.
1.1202 + ///
1.1203 + /// This method runs the %MaxWeightedMatching algorithm.
1.1204 + ///
1.1205 + /// \note mwm.run() is just a shortcut of the following code.
1.1206 + /// \code
1.1207 + /// mwm.init();
1.1208 + /// mwm.start();
1.1209 + /// \endcode
1.1210 + void run() {
1.1211 + init();
1.1212 + start();
1.1213 + }
1.1214 +
1.1215 + /// @}
1.1216 +
1.1217 + /// \name Primal solution
1.1218 + /// Functions for get the primal solution, ie. the matching.
1.1219 +
1.1220 + /// @{
1.1221 +
1.1222 + /// \brief Returns the matching value.
1.1223 + ///
1.1224 + /// Returns the matching value.
1.1225 + Value matchingValue() const {
1.1226 + Value sum = 0;
1.1227 + for (NodeIt n(_ugraph); n != INVALID; ++n) {
1.1228 + if ((*_matching)[n] != INVALID) {
1.1229 + sum += _weight[(*_matching)[n]];
1.1230 + }
1.1231 + }
1.1232 + return sum /= 2;
1.1233 + }
1.1234 +
1.1235 + /// \brief Returns true when the edge is in the matching.
1.1236 + ///
1.1237 + /// Returns true when the edge is in the matching.
1.1238 + bool matching(const UEdge& edge) const {
1.1239 + return (*_matching)[_ugraph.source(edge)] == _ugraph.direct(edge, true);
1.1240 + }
1.1241 +
1.1242 + /// \brief Returns the incident matching edge.
1.1243 + ///
1.1244 + /// Returns the incident matching edge from given node. If the
1.1245 + /// node is not matched then it gives back \c INVALID.
1.1246 + Edge matching(const Node& node) const {
1.1247 + return (*_matching)[node];
1.1248 + }
1.1249 +
1.1250 + /// \brief Returns the mate of the node.
1.1251 + ///
1.1252 + /// Returns the adjancent node in a mathcing edge. If the node is
1.1253 + /// not matched then it gives back \c INVALID.
1.1254 + Node mate(const Node& node) const {
1.1255 + return (*_matching)[node] != INVALID ?
1.1256 + _ugraph.target((*_matching)[node]) : INVALID;
1.1257 + }
1.1258 +
1.1259 + /// @}
1.1260 +
1.1261 + /// \name Dual solution
1.1262 + /// Functions for get the dual solution.
1.1263 +
1.1264 + /// @{
1.1265 +
1.1266 + /// \brief Returns the value of the dual solution.
1.1267 + ///
1.1268 + /// Returns the value of the dual solution. It should be equal to
1.1269 + /// the primal value scaled by \ref dualScale "dual scale".
1.1270 + Value dualValue() const {
1.1271 + Value sum = 0;
1.1272 + for (NodeIt n(_ugraph); n != INVALID; ++n) {
1.1273 + sum += nodeValue(n);
1.1274 + }
1.1275 + for (int i = 0; i < blossomNum(); ++i) {
1.1276 + sum += blossomValue(i) * (blossomSize(i) / 2);
1.1277 + }
1.1278 + return sum;
1.1279 + }
1.1280 +
1.1281 + /// \brief Returns the value of the node.
1.1282 + ///
1.1283 + /// Returns the the value of the node.
1.1284 + Value nodeValue(const Node& n) const {
1.1285 + return (*_node_potential)[n];
1.1286 + }
1.1287 +
1.1288 + /// \brief Returns the number of the blossoms in the basis.
1.1289 + ///
1.1290 + /// Returns the number of the blossoms in the basis.
1.1291 + /// \see BlossomIt
1.1292 + int blossomNum() const {
1.1293 + return _blossom_potential.size();
1.1294 + }
1.1295 +
1.1296 +
1.1297 + /// \brief Returns the number of the nodes in the blossom.
1.1298 + ///
1.1299 + /// Returns the number of the nodes in the blossom.
1.1300 + int blossomSize(int k) const {
1.1301 + return _blossom_potential[k].end - _blossom_potential[k].begin;
1.1302 + }
1.1303 +
1.1304 + /// \brief Returns the value of the blossom.
1.1305 + ///
1.1306 + /// Returns the the value of the blossom.
1.1307 + /// \see BlossomIt
1.1308 + Value blossomValue(int k) const {
1.1309 + return _blossom_potential[k].value;
1.1310 + }
1.1311 +
1.1312 + /// \brief Lemon iterator for get the items of the blossom.
1.1313 + ///
1.1314 + /// Lemon iterator for get the nodes of the blossom. This class
1.1315 + /// provides a common style lemon iterator which gives back a
1.1316 + /// subset of the nodes.
1.1317 + class BlossomIt {
1.1318 + public:
1.1319 +
1.1320 + /// \brief Constructor.
1.1321 + ///
1.1322 + /// Constructor for get the nodes of the variable.
1.1323 + BlossomIt(const MaxWeightedMatching& algorithm, int variable)
1.1324 + : _algorithm(&algorithm)
1.1325 + {
1.1326 + _index = _algorithm->_blossom_potential[variable].begin;
1.1327 + _last = _algorithm->_blossom_potential[variable].end;
1.1328 + }
1.1329 +
1.1330 + /// \brief Invalid constructor.
1.1331 + ///
1.1332 + /// Invalid constructor.
1.1333 + BlossomIt(Invalid) : _index(-1) {}
1.1334 +
1.1335 + /// \brief Conversion to node.
1.1336 + ///
1.1337 + /// Conversion to node.
1.1338 + operator Node() const {
1.1339 + return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
1.1340 + }
1.1341 +
1.1342 + /// \brief Increment operator.
1.1343 + ///
1.1344 + /// Increment operator.
1.1345 + BlossomIt& operator++() {
1.1346 + ++_index;
1.1347 + if (_index == _last) {
1.1348 + _index = -1;
1.1349 + }
1.1350 + return *this;
1.1351 + }
1.1352 +
1.1353 + bool operator==(const BlossomIt& it) const {
1.1354 + return _index == it._index;
1.1355 + }
1.1356 + bool operator!=(const BlossomIt& it) const {
1.1357 + return _index != it._index;
1.1358 + }
1.1359 +
1.1360 + private:
1.1361 + const MaxWeightedMatching* _algorithm;
1.1362 + int _last;
1.1363 + int _index;
1.1364 + };
1.1365 +
1.1366 + /// @}
1.1367 +
1.1368 + };
1.1369 +
1.1370 + /// \ingroup matching
1.1371 + ///
1.1372 + /// \brief Weighted perfect matching in general undirected graphs
1.1373 + ///
1.1374 + /// This class provides an efficient implementation of Edmond's
1.1375 + /// maximum weighted perfecr matching algorithm. The implementation
1.1376 + /// is based on extensive use of priority queues and provides
1.1377 + /// \f$O(nm\log(n))\f$ time complexity.
1.1378 + ///
1.1379 + /// The maximum weighted matching problem is to find undirected
1.1380 + /// edges in the graph with maximum overall weight and no two of
1.1381 + /// them shares their endpoints and covers all nodes. The problem
1.1382 + /// can be formulated with the next linear program:
1.1383 + /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1.1384 + ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
1.1385 + /// \f[x_e \ge 0\quad \forall e\in E\f]
1.1386 + /// \f[\max \sum_{e\in E}x_ew_e\f]
1.1387 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1.1388 + /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both endpoints in
1.1389 + /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
1.1390 + /// the nodes.
1.1391 + ///
1.1392 + /// The algorithm calculates an optimal matching and a proof of the
1.1393 + /// optimality. The solution of the dual problem can be used to check
1.1394 + /// the result of the algorithm. The dual linear problem is the next:
1.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]
1.1396 + /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
1.1397 + /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
1.1398 + ///
1.1399 + /// The algorithm can be executed with \c run() or the \c init() and
1.1400 + /// then the \c start() member functions. After it the matching can
1.1401 + /// be asked with \c matching() or mate() functions. The dual
1.1402 + /// solution can be get with \c nodeValue(), \c blossomNum() and \c
1.1403 + /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
1.1404 + /// "BlossomIt" nested class which is able to iterate on the nodes
1.1405 + /// of a blossom. If the value type is integral then the dual
1.1406 + /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
1.1407 + ///
1.1408 + /// \author Balazs Dezso
1.1409 + template <typename _UGraph,
1.1410 + typename _WeightMap = typename _UGraph::template UEdgeMap<int> >
1.1411 + class MaxWeightedPerfectMatching {
1.1412 + public:
1.1413 +
1.1414 + typedef _UGraph UGraph;
1.1415 + typedef _WeightMap WeightMap;
1.1416 + typedef typename WeightMap::Value Value;
1.1417 +
1.1418 + /// \brief Scaling factor for dual solution
1.1419 + ///
1.1420 + /// Scaling factor for dual solution, it is equal to 4 or 1
1.1421 + /// according to the value type.
1.1422 + static const int dualScale =
1.1423 + std::numeric_limits<Value>::is_integer ? 4 : 1;
1.1424 +
1.1425 + typedef typename UGraph::template NodeMap<typename UGraph::Edge>
1.1426 + MatchingMap;
1.1427 +
1.1428 + private:
1.1429 +
1.1430 + UGRAPH_TYPEDEFS(typename UGraph);
1.1431 +
1.1432 + typedef typename UGraph::template NodeMap<Value> NodePotential;
1.1433 + typedef std::vector<Node> BlossomNodeList;
1.1434 +
1.1435 + struct BlossomVariable {
1.1436 + int begin, end;
1.1437 + Value value;
1.1438 +
1.1439 + BlossomVariable(int _begin, int _end, Value _value)
1.1440 + : begin(_begin), end(_end), value(_value) {}
1.1441 +
1.1442 + };
1.1443 +
1.1444 + typedef std::vector<BlossomVariable> BlossomPotential;
1.1445 +
1.1446 + const UGraph& _ugraph;
1.1447 + const WeightMap& _weight;
1.1448 +
1.1449 + MatchingMap* _matching;
1.1450 +
1.1451 + NodePotential* _node_potential;
1.1452 +
1.1453 + BlossomPotential _blossom_potential;
1.1454 + BlossomNodeList _blossom_node_list;
1.1455 +
1.1456 + int _node_num;
1.1457 + int _blossom_num;
1.1458 +
1.1459 + typedef typename UGraph::template NodeMap<int> NodeIntMap;
1.1460 + typedef typename UGraph::template EdgeMap<int> EdgeIntMap;
1.1461 + typedef typename UGraph::template UEdgeMap<int> UEdgeIntMap;
1.1462 + typedef IntegerMap<int> IntIntMap;
1.1463 +
1.1464 + enum Status {
1.1465 + EVEN = -1, MATCHED = 0, ODD = 1
1.1466 + };
1.1467 +
1.1468 + typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
1.1469 + struct BlossomData {
1.1470 + int tree;
1.1471 + Status status;
1.1472 + Edge pred, next;
1.1473 + Value pot, offset;
1.1474 + };
1.1475 +
1.1476 + NodeIntMap *_blossom_index;
1.1477 + BlossomSet *_blossom_set;
1.1478 + IntegerMap<BlossomData>* _blossom_data;
1.1479 +
1.1480 + NodeIntMap *_node_index;
1.1481 + EdgeIntMap *_node_heap_index;
1.1482 +
1.1483 + struct NodeData {
1.1484 +
1.1485 + NodeData(EdgeIntMap& node_heap_index)
1.1486 + : heap(node_heap_index) {}
1.1487 +
1.1488 + int blossom;
1.1489 + Value pot;
1.1490 + BinHeap<Value, EdgeIntMap> heap;
1.1491 + std::map<int, Edge> heap_index;
1.1492 +
1.1493 + int tree;
1.1494 + };
1.1495 +
1.1496 + IntegerMap<NodeData>* _node_data;
1.1497 +
1.1498 + typedef ExtendFindEnum<IntIntMap> TreeSet;
1.1499 +
1.1500 + IntIntMap *_tree_set_index;
1.1501 + TreeSet *_tree_set;
1.1502 +
1.1503 + IntIntMap *_delta2_index;
1.1504 + BinHeap<Value, IntIntMap> *_delta2;
1.1505 +
1.1506 + UEdgeIntMap *_delta3_index;
1.1507 + BinHeap<Value, UEdgeIntMap> *_delta3;
1.1508 +
1.1509 + IntIntMap *_delta4_index;
1.1510 + BinHeap<Value, IntIntMap> *_delta4;
1.1511 +
1.1512 + Value _delta_sum;
1.1513 +
1.1514 + void createStructures() {
1.1515 + _node_num = countNodes(_ugraph);
1.1516 + _blossom_num = _node_num * 3 / 2;
1.1517 +
1.1518 + if (!_matching) {
1.1519 + _matching = new MatchingMap(_ugraph);
1.1520 + }
1.1521 + if (!_node_potential) {
1.1522 + _node_potential = new NodePotential(_ugraph);
1.1523 + }
1.1524 + if (!_blossom_set) {
1.1525 + _blossom_index = new NodeIntMap(_ugraph);
1.1526 + _blossom_set = new BlossomSet(*_blossom_index);
1.1527 + _blossom_data = new IntegerMap<BlossomData>(_blossom_num);
1.1528 + }
1.1529 +
1.1530 + if (!_node_index) {
1.1531 + _node_index = new NodeIntMap(_ugraph);
1.1532 + _node_heap_index = new EdgeIntMap(_ugraph);
1.1533 + _node_data = new IntegerMap<NodeData>(_node_num,
1.1534 + NodeData(*_node_heap_index));
1.1535 + }
1.1536 +
1.1537 + if (!_tree_set) {
1.1538 + _tree_set_index = new IntIntMap(_blossom_num);
1.1539 + _tree_set = new TreeSet(*_tree_set_index);
1.1540 + }
1.1541 + if (!_delta2) {
1.1542 + _delta2_index = new IntIntMap(_blossom_num);
1.1543 + _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
1.1544 + }
1.1545 + if (!_delta3) {
1.1546 + _delta3_index = new UEdgeIntMap(_ugraph);
1.1547 + _delta3 = new BinHeap<Value, UEdgeIntMap>(*_delta3_index);
1.1548 + }
1.1549 + if (!_delta4) {
1.1550 + _delta4_index = new IntIntMap(_blossom_num);
1.1551 + _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
1.1552 + }
1.1553 + }
1.1554 +
1.1555 + void destroyStructures() {
1.1556 + _node_num = countNodes(_ugraph);
1.1557 + _blossom_num = _node_num * 3 / 2;
1.1558 +
1.1559 + if (_matching) {
1.1560 + delete _matching;
1.1561 + }
1.1562 + if (_node_potential) {
1.1563 + delete _node_potential;
1.1564 + }
1.1565 + if (_blossom_set) {
1.1566 + delete _blossom_index;
1.1567 + delete _blossom_set;
1.1568 + delete _blossom_data;
1.1569 + }
1.1570 +
1.1571 + if (_node_index) {
1.1572 + delete _node_index;
1.1573 + delete _node_heap_index;
1.1574 + delete _node_data;
1.1575 + }
1.1576 +
1.1577 + if (_tree_set) {
1.1578 + delete _tree_set_index;
1.1579 + delete _tree_set;
1.1580 + }
1.1581 + if (_delta2) {
1.1582 + delete _delta2_index;
1.1583 + delete _delta2;
1.1584 + }
1.1585 + if (_delta3) {
1.1586 + delete _delta3_index;
1.1587 + delete _delta3;
1.1588 + }
1.1589 + if (_delta4) {
1.1590 + delete _delta4_index;
1.1591 + delete _delta4;
1.1592 + }
1.1593 + }
1.1594 +
1.1595 + void matchedToEven(int blossom, int tree) {
1.1596 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.1597 + _delta2->erase(blossom);
1.1598 + }
1.1599 +
1.1600 + if (!_blossom_set->trivial(blossom)) {
1.1601 + (*_blossom_data)[blossom].pot -=
1.1602 + 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
1.1603 + }
1.1604 +
1.1605 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.1606 + n != INVALID; ++n) {
1.1607 +
1.1608 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
1.1609 + int ni = (*_node_index)[n];
1.1610 +
1.1611 + (*_node_data)[ni].heap.clear();
1.1612 + (*_node_data)[ni].heap_index.clear();
1.1613 +
1.1614 + (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
1.1615 +
1.1616 + for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
1.1617 + Node v = _ugraph.source(e);
1.1618 + int vb = _blossom_set->find(v);
1.1619 + int vi = (*_node_index)[v];
1.1620 +
1.1621 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.1622 + dualScale * _weight[e];
1.1623 +
1.1624 + if ((*_blossom_data)[vb].status == EVEN) {
1.1625 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1.1626 + _delta3->push(e, rw / 2);
1.1627 + }
1.1628 + } else {
1.1629 + typename std::map<int, Edge>::iterator it =
1.1630 + (*_node_data)[vi].heap_index.find(tree);
1.1631 +
1.1632 + if (it != (*_node_data)[vi].heap_index.end()) {
1.1633 + if ((*_node_data)[vi].heap[it->second] > rw) {
1.1634 + (*_node_data)[vi].heap.replace(it->second, e);
1.1635 + (*_node_data)[vi].heap.decrease(e, rw);
1.1636 + it->second = e;
1.1637 + }
1.1638 + } else {
1.1639 + (*_node_data)[vi].heap.push(e, rw);
1.1640 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1.1641 + }
1.1642 +
1.1643 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1.1644 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1.1645 +
1.1646 + if ((*_blossom_data)[vb].status == MATCHED) {
1.1647 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
1.1648 + _delta2->push(vb, _blossom_set->classPrio(vb) -
1.1649 + (*_blossom_data)[vb].offset);
1.1650 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1.1651 + (*_blossom_data)[vb].offset){
1.1652 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1.1653 + (*_blossom_data)[vb].offset);
1.1654 + }
1.1655 + }
1.1656 + }
1.1657 + }
1.1658 + }
1.1659 + }
1.1660 + (*_blossom_data)[blossom].offset = 0;
1.1661 + }
1.1662 +
1.1663 + void matchedToOdd(int blossom) {
1.1664 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.1665 + _delta2->erase(blossom);
1.1666 + }
1.1667 + (*_blossom_data)[blossom].offset += _delta_sum;
1.1668 + if (!_blossom_set->trivial(blossom)) {
1.1669 + _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
1.1670 + (*_blossom_data)[blossom].offset);
1.1671 + }
1.1672 + }
1.1673 +
1.1674 + void evenToMatched(int blossom, int tree) {
1.1675 + if (!_blossom_set->trivial(blossom)) {
1.1676 + (*_blossom_data)[blossom].pot += 2 * _delta_sum;
1.1677 + }
1.1678 +
1.1679 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.1680 + n != INVALID; ++n) {
1.1681 + int ni = (*_node_index)[n];
1.1682 + (*_node_data)[ni].pot -= _delta_sum;
1.1683 +
1.1684 + for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
1.1685 + Node v = _ugraph.source(e);
1.1686 + int vb = _blossom_set->find(v);
1.1687 + int vi = (*_node_index)[v];
1.1688 +
1.1689 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.1690 + dualScale * _weight[e];
1.1691 +
1.1692 + if (vb == blossom) {
1.1693 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.1694 + _delta3->erase(e);
1.1695 + }
1.1696 + } else if ((*_blossom_data)[vb].status == EVEN) {
1.1697 +
1.1698 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.1699 + _delta3->erase(e);
1.1700 + }
1.1701 +
1.1702 + int vt = _tree_set->find(vb);
1.1703 +
1.1704 + if (vt != tree) {
1.1705 +
1.1706 + Edge r = _ugraph.oppositeEdge(e);
1.1707 +
1.1708 + typename std::map<int, Edge>::iterator it =
1.1709 + (*_node_data)[ni].heap_index.find(vt);
1.1710 +
1.1711 + if (it != (*_node_data)[ni].heap_index.end()) {
1.1712 + if ((*_node_data)[ni].heap[it->second] > rw) {
1.1713 + (*_node_data)[ni].heap.replace(it->second, r);
1.1714 + (*_node_data)[ni].heap.decrease(r, rw);
1.1715 + it->second = r;
1.1716 + }
1.1717 + } else {
1.1718 + (*_node_data)[ni].heap.push(r, rw);
1.1719 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1.1720 + }
1.1721 +
1.1722 + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1.1723 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1.1724 +
1.1725 + if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1.1726 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.1727 + (*_blossom_data)[blossom].offset);
1.1728 + } else if ((*_delta2)[blossom] >
1.1729 + _blossom_set->classPrio(blossom) -
1.1730 + (*_blossom_data)[blossom].offset){
1.1731 + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1.1732 + (*_blossom_data)[blossom].offset);
1.1733 + }
1.1734 + }
1.1735 + }
1.1736 + } else {
1.1737 +
1.1738 + typename std::map<int, Edge>::iterator it =
1.1739 + (*_node_data)[vi].heap_index.find(tree);
1.1740 +
1.1741 + if (it != (*_node_data)[vi].heap_index.end()) {
1.1742 + (*_node_data)[vi].heap.erase(it->second);
1.1743 + (*_node_data)[vi].heap_index.erase(it);
1.1744 + if ((*_node_data)[vi].heap.empty()) {
1.1745 + _blossom_set->increase(v, std::numeric_limits<Value>::max());
1.1746 + } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1.1747 + _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1.1748 + }
1.1749 +
1.1750 + if ((*_blossom_data)[vb].status == MATCHED) {
1.1751 + if (_blossom_set->classPrio(vb) ==
1.1752 + std::numeric_limits<Value>::max()) {
1.1753 + _delta2->erase(vb);
1.1754 + } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1.1755 + (*_blossom_data)[vb].offset) {
1.1756 + _delta2->increase(vb, _blossom_set->classPrio(vb) -
1.1757 + (*_blossom_data)[vb].offset);
1.1758 + }
1.1759 + }
1.1760 + }
1.1761 + }
1.1762 + }
1.1763 + }
1.1764 + }
1.1765 +
1.1766 + void oddToMatched(int blossom) {
1.1767 + (*_blossom_data)[blossom].offset -= _delta_sum;
1.1768 +
1.1769 + if (_blossom_set->classPrio(blossom) !=
1.1770 + std::numeric_limits<Value>::max()) {
1.1771 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.1772 + (*_blossom_data)[blossom].offset);
1.1773 + }
1.1774 +
1.1775 + if (!_blossom_set->trivial(blossom)) {
1.1776 + _delta4->erase(blossom);
1.1777 + }
1.1778 + }
1.1779 +
1.1780 + void oddToEven(int blossom, int tree) {
1.1781 + if (!_blossom_set->trivial(blossom)) {
1.1782 + _delta4->erase(blossom);
1.1783 + (*_blossom_data)[blossom].pot -=
1.1784 + 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1.1785 + }
1.1786 +
1.1787 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.1788 + n != INVALID; ++n) {
1.1789 + int ni = (*_node_index)[n];
1.1790 +
1.1791 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
1.1792 +
1.1793 + (*_node_data)[ni].heap.clear();
1.1794 + (*_node_data)[ni].heap_index.clear();
1.1795 + (*_node_data)[ni].pot +=
1.1796 + 2 * _delta_sum - (*_blossom_data)[blossom].offset;
1.1797 +
1.1798 + for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
1.1799 + Node v = _ugraph.source(e);
1.1800 + int vb = _blossom_set->find(v);
1.1801 + int vi = (*_node_index)[v];
1.1802 +
1.1803 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.1804 + dualScale * _weight[e];
1.1805 +
1.1806 + if ((*_blossom_data)[vb].status == EVEN) {
1.1807 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1.1808 + _delta3->push(e, rw / 2);
1.1809 + }
1.1810 + } else {
1.1811 +
1.1812 + typename std::map<int, Edge>::iterator it =
1.1813 + (*_node_data)[vi].heap_index.find(tree);
1.1814 +
1.1815 + if (it != (*_node_data)[vi].heap_index.end()) {
1.1816 + if ((*_node_data)[vi].heap[it->second] > rw) {
1.1817 + (*_node_data)[vi].heap.replace(it->second, e);
1.1818 + (*_node_data)[vi].heap.decrease(e, rw);
1.1819 + it->second = e;
1.1820 + }
1.1821 + } else {
1.1822 + (*_node_data)[vi].heap.push(e, rw);
1.1823 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1.1824 + }
1.1825 +
1.1826 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1.1827 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1.1828 +
1.1829 + if ((*_blossom_data)[vb].status == MATCHED) {
1.1830 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
1.1831 + _delta2->push(vb, _blossom_set->classPrio(vb) -
1.1832 + (*_blossom_data)[vb].offset);
1.1833 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1.1834 + (*_blossom_data)[vb].offset) {
1.1835 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1.1836 + (*_blossom_data)[vb].offset);
1.1837 + }
1.1838 + }
1.1839 + }
1.1840 + }
1.1841 + }
1.1842 + }
1.1843 + (*_blossom_data)[blossom].offset = 0;
1.1844 + }
1.1845 +
1.1846 + void alternatePath(int even, int tree) {
1.1847 + int odd;
1.1848 +
1.1849 + evenToMatched(even, tree);
1.1850 + (*_blossom_data)[even].status = MATCHED;
1.1851 +
1.1852 + while ((*_blossom_data)[even].pred != INVALID) {
1.1853 + odd = _blossom_set->find(_ugraph.target((*_blossom_data)[even].pred));
1.1854 + (*_blossom_data)[odd].status = MATCHED;
1.1855 + oddToMatched(odd);
1.1856 + (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
1.1857 +
1.1858 + even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].pred));
1.1859 + (*_blossom_data)[even].status = MATCHED;
1.1860 + evenToMatched(even, tree);
1.1861 + (*_blossom_data)[even].next =
1.1862 + _ugraph.oppositeEdge((*_blossom_data)[odd].pred);
1.1863 + }
1.1864 +
1.1865 + }
1.1866 +
1.1867 + void destroyTree(int tree) {
1.1868 + for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
1.1869 + if ((*_blossom_data)[b].status == EVEN) {
1.1870 + (*_blossom_data)[b].status = MATCHED;
1.1871 + evenToMatched(b, tree);
1.1872 + } else if ((*_blossom_data)[b].status == ODD) {
1.1873 + (*_blossom_data)[b].status = MATCHED;
1.1874 + oddToMatched(b);
1.1875 + }
1.1876 + }
1.1877 + _tree_set->eraseClass(tree);
1.1878 + }
1.1879 +
1.1880 + void augmentOnEdge(const UEdge& edge) {
1.1881 +
1.1882 + int left = _blossom_set->find(_ugraph.source(edge));
1.1883 + int right = _blossom_set->find(_ugraph.target(edge));
1.1884 +
1.1885 + int left_tree = _tree_set->find(left);
1.1886 + alternatePath(left, left_tree);
1.1887 + destroyTree(left_tree);
1.1888 +
1.1889 + int right_tree = _tree_set->find(right);
1.1890 + alternatePath(right, right_tree);
1.1891 + destroyTree(right_tree);
1.1892 +
1.1893 + (*_blossom_data)[left].next = _ugraph.direct(edge, true);
1.1894 + (*_blossom_data)[right].next = _ugraph.direct(edge, false);
1.1895 + }
1.1896 +
1.1897 + void extendOnEdge(const Edge& edge) {
1.1898 + int base = _blossom_set->find(_ugraph.target(edge));
1.1899 + int tree = _tree_set->find(base);
1.1900 +
1.1901 + int odd = _blossom_set->find(_ugraph.source(edge));
1.1902 + _tree_set->insert(odd, tree);
1.1903 + (*_blossom_data)[odd].status = ODD;
1.1904 + matchedToOdd(odd);
1.1905 + (*_blossom_data)[odd].pred = edge;
1.1906 +
1.1907 + int even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].next));
1.1908 + (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
1.1909 + _tree_set->insert(even, tree);
1.1910 + (*_blossom_data)[even].status = EVEN;
1.1911 + matchedToEven(even, tree);
1.1912 + }
1.1913 +
1.1914 + void shrinkOnEdge(const UEdge& uedge, int tree) {
1.1915 + int nca = -1;
1.1916 + std::vector<int> left_path, right_path;
1.1917 +
1.1918 + {
1.1919 + std::set<int> left_set, right_set;
1.1920 + int left = _blossom_set->find(_ugraph.source(uedge));
1.1921 + left_path.push_back(left);
1.1922 + left_set.insert(left);
1.1923 +
1.1924 + int right = _blossom_set->find(_ugraph.target(uedge));
1.1925 + right_path.push_back(right);
1.1926 + right_set.insert(right);
1.1927 +
1.1928 + while (true) {
1.1929 +
1.1930 + if ((*_blossom_data)[left].pred == INVALID) break;
1.1931 +
1.1932 + left =
1.1933 + _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
1.1934 + left_path.push_back(left);
1.1935 + left =
1.1936 + _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
1.1937 + left_path.push_back(left);
1.1938 +
1.1939 + left_set.insert(left);
1.1940 +
1.1941 + if (right_set.find(left) != right_set.end()) {
1.1942 + nca = left;
1.1943 + break;
1.1944 + }
1.1945 +
1.1946 + if ((*_blossom_data)[right].pred == INVALID) break;
1.1947 +
1.1948 + right =
1.1949 + _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
1.1950 + right_path.push_back(right);
1.1951 + right =
1.1952 + _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
1.1953 + right_path.push_back(right);
1.1954 +
1.1955 + right_set.insert(right);
1.1956 +
1.1957 + if (left_set.find(right) != left_set.end()) {
1.1958 + nca = right;
1.1959 + break;
1.1960 + }
1.1961 +
1.1962 + }
1.1963 +
1.1964 + if (nca == -1) {
1.1965 + if ((*_blossom_data)[left].pred == INVALID) {
1.1966 + nca = right;
1.1967 + while (left_set.find(nca) == left_set.end()) {
1.1968 + nca =
1.1969 + _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
1.1970 + right_path.push_back(nca);
1.1971 + nca =
1.1972 + _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
1.1973 + right_path.push_back(nca);
1.1974 + }
1.1975 + } else {
1.1976 + nca = left;
1.1977 + while (right_set.find(nca) == right_set.end()) {
1.1978 + nca =
1.1979 + _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
1.1980 + left_path.push_back(nca);
1.1981 + nca =
1.1982 + _blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
1.1983 + left_path.push_back(nca);
1.1984 + }
1.1985 + }
1.1986 + }
1.1987 + }
1.1988 +
1.1989 + std::vector<int> subblossoms;
1.1990 + Edge prev;
1.1991 +
1.1992 + prev = _ugraph.direct(uedge, true);
1.1993 + for (int i = 0; left_path[i] != nca; i += 2) {
1.1994 + subblossoms.push_back(left_path[i]);
1.1995 + (*_blossom_data)[left_path[i]].next = prev;
1.1996 + _tree_set->erase(left_path[i]);
1.1997 +
1.1998 + subblossoms.push_back(left_path[i + 1]);
1.1999 + (*_blossom_data)[left_path[i + 1]].status = EVEN;
1.2000 + oddToEven(left_path[i + 1], tree);
1.2001 + _tree_set->erase(left_path[i + 1]);
1.2002 + prev = _ugraph.oppositeEdge((*_blossom_data)[left_path[i + 1]].pred);
1.2003 + }
1.2004 +
1.2005 + int k = 0;
1.2006 + while (right_path[k] != nca) ++k;
1.2007 +
1.2008 + subblossoms.push_back(nca);
1.2009 + (*_blossom_data)[nca].next = prev;
1.2010 +
1.2011 + for (int i = k - 2; i >= 0; i -= 2) {
1.2012 + subblossoms.push_back(right_path[i + 1]);
1.2013 + (*_blossom_data)[right_path[i + 1]].status = EVEN;
1.2014 + oddToEven(right_path[i + 1], tree);
1.2015 + _tree_set->erase(right_path[i + 1]);
1.2016 +
1.2017 + (*_blossom_data)[right_path[i + 1]].next =
1.2018 + (*_blossom_data)[right_path[i + 1]].pred;
1.2019 +
1.2020 + subblossoms.push_back(right_path[i]);
1.2021 + _tree_set->erase(right_path[i]);
1.2022 + }
1.2023 +
1.2024 + int surface =
1.2025 + _blossom_set->join(subblossoms.begin(), subblossoms.end());
1.2026 +
1.2027 + for (int i = 0; i < int(subblossoms.size()); ++i) {
1.2028 + if (!_blossom_set->trivial(subblossoms[i])) {
1.2029 + (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
1.2030 + }
1.2031 + (*_blossom_data)[subblossoms[i]].status = MATCHED;
1.2032 + }
1.2033 +
1.2034 + (*_blossom_data)[surface].pot = -2 * _delta_sum;
1.2035 + (*_blossom_data)[surface].offset = 0;
1.2036 + (*_blossom_data)[surface].status = EVEN;
1.2037 + (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
1.2038 + (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
1.2039 +
1.2040 + _tree_set->insert(surface, tree);
1.2041 + _tree_set->erase(nca);
1.2042 + }
1.2043 +
1.2044 + void splitBlossom(int blossom) {
1.2045 + Edge next = (*_blossom_data)[blossom].next;
1.2046 + Edge pred = (*_blossom_data)[blossom].pred;
1.2047 +
1.2048 + int tree = _tree_set->find(blossom);
1.2049 +
1.2050 + (*_blossom_data)[blossom].status = MATCHED;
1.2051 + oddToMatched(blossom);
1.2052 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.2053 + _delta2->erase(blossom);
1.2054 + }
1.2055 +
1.2056 + std::vector<int> subblossoms;
1.2057 + _blossom_set->split(blossom, std::back_inserter(subblossoms));
1.2058 +
1.2059 + Value offset = (*_blossom_data)[blossom].offset;
1.2060 + int b = _blossom_set->find(_ugraph.source(pred));
1.2061 + int d = _blossom_set->find(_ugraph.source(next));
1.2062 +
1.2063 + int ib, id;
1.2064 + for (int i = 0; i < int(subblossoms.size()); ++i) {
1.2065 + if (subblossoms[i] == b) ib = i;
1.2066 + if (subblossoms[i] == d) id = i;
1.2067 +
1.2068 + (*_blossom_data)[subblossoms[i]].offset = offset;
1.2069 + if (!_blossom_set->trivial(subblossoms[i])) {
1.2070 + (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
1.2071 + }
1.2072 + if (_blossom_set->classPrio(subblossoms[i]) !=
1.2073 + std::numeric_limits<Value>::max()) {
1.2074 + _delta2->push(subblossoms[i],
1.2075 + _blossom_set->classPrio(subblossoms[i]) -
1.2076 + (*_blossom_data)[subblossoms[i]].offset);
1.2077 + }
1.2078 + }
1.2079 +
1.2080 + if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
1.2081 + for (int i = (id + 1) % subblossoms.size();
1.2082 + i != ib; i = (i + 2) % subblossoms.size()) {
1.2083 + int sb = subblossoms[i];
1.2084 + int tb = subblossoms[(i + 1) % subblossoms.size()];
1.2085 + (*_blossom_data)[sb].next =
1.2086 + _ugraph.oppositeEdge((*_blossom_data)[tb].next);
1.2087 + }
1.2088 +
1.2089 + for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
1.2090 + int sb = subblossoms[i];
1.2091 + int tb = subblossoms[(i + 1) % subblossoms.size()];
1.2092 + int ub = subblossoms[(i + 2) % subblossoms.size()];
1.2093 +
1.2094 + (*_blossom_data)[sb].status = ODD;
1.2095 + matchedToOdd(sb);
1.2096 + _tree_set->insert(sb, tree);
1.2097 + (*_blossom_data)[sb].pred = pred;
1.2098 + (*_blossom_data)[sb].next =
1.2099 + _ugraph.oppositeEdge((*_blossom_data)[tb].next);
1.2100 +
1.2101 + pred = (*_blossom_data)[ub].next;
1.2102 +
1.2103 + (*_blossom_data)[tb].status = EVEN;
1.2104 + matchedToEven(tb, tree);
1.2105 + _tree_set->insert(tb, tree);
1.2106 + (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
1.2107 + }
1.2108 +
1.2109 + (*_blossom_data)[subblossoms[id]].status = ODD;
1.2110 + matchedToOdd(subblossoms[id]);
1.2111 + _tree_set->insert(subblossoms[id], tree);
1.2112 + (*_blossom_data)[subblossoms[id]].next = next;
1.2113 + (*_blossom_data)[subblossoms[id]].pred = pred;
1.2114 +
1.2115 + } else {
1.2116 +
1.2117 + for (int i = (ib + 1) % subblossoms.size();
1.2118 + i != id; i = (i + 2) % subblossoms.size()) {
1.2119 + int sb = subblossoms[i];
1.2120 + int tb = subblossoms[(i + 1) % subblossoms.size()];
1.2121 + (*_blossom_data)[sb].next =
1.2122 + _ugraph.oppositeEdge((*_blossom_data)[tb].next);
1.2123 + }
1.2124 +
1.2125 + for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
1.2126 + int sb = subblossoms[i];
1.2127 + int tb = subblossoms[(i + 1) % subblossoms.size()];
1.2128 + int ub = subblossoms[(i + 2) % subblossoms.size()];
1.2129 +
1.2130 + (*_blossom_data)[sb].status = ODD;
1.2131 + matchedToOdd(sb);
1.2132 + _tree_set->insert(sb, tree);
1.2133 + (*_blossom_data)[sb].next = next;
1.2134 + (*_blossom_data)[sb].pred =
1.2135 + _ugraph.oppositeEdge((*_blossom_data)[tb].next);
1.2136 +
1.2137 + (*_blossom_data)[tb].status = EVEN;
1.2138 + matchedToEven(tb, tree);
1.2139 + _tree_set->insert(tb, tree);
1.2140 + (*_blossom_data)[tb].pred =
1.2141 + (*_blossom_data)[tb].next =
1.2142 + _ugraph.oppositeEdge((*_blossom_data)[ub].next);
1.2143 + next = (*_blossom_data)[ub].next;
1.2144 + }
1.2145 +
1.2146 + (*_blossom_data)[subblossoms[ib]].status = ODD;
1.2147 + matchedToOdd(subblossoms[ib]);
1.2148 + _tree_set->insert(subblossoms[ib], tree);
1.2149 + (*_blossom_data)[subblossoms[ib]].next = next;
1.2150 + (*_blossom_data)[subblossoms[ib]].pred = pred;
1.2151 + }
1.2152 + _tree_set->erase(blossom);
1.2153 + }
1.2154 +
1.2155 + void extractBlossom(int blossom, const Node& base, const Edge& matching) {
1.2156 + if (_blossom_set->trivial(blossom)) {
1.2157 + int bi = (*_node_index)[base];
1.2158 + Value pot = (*_node_data)[bi].pot;
1.2159 +
1.2160 + _matching->set(base, matching);
1.2161 + _blossom_node_list.push_back(base);
1.2162 + _node_potential->set(base, pot);
1.2163 + } else {
1.2164 +
1.2165 + Value pot = (*_blossom_data)[blossom].pot;
1.2166 + int bn = _blossom_node_list.size();
1.2167 +
1.2168 + std::vector<int> subblossoms;
1.2169 + _blossom_set->split(blossom, std::back_inserter(subblossoms));
1.2170 + int b = _blossom_set->find(base);
1.2171 + int ib = -1;
1.2172 + for (int i = 0; i < int(subblossoms.size()); ++i) {
1.2173 + if (subblossoms[i] == b) { ib = i; break; }
1.2174 + }
1.2175 +
1.2176 + for (int i = 1; i < int(subblossoms.size()); i += 2) {
1.2177 + int sb = subblossoms[(ib + i) % subblossoms.size()];
1.2178 + int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
1.2179 +
1.2180 + Edge m = (*_blossom_data)[tb].next;
1.2181 + extractBlossom(sb, _ugraph.target(m), _ugraph.oppositeEdge(m));
1.2182 + extractBlossom(tb, _ugraph.source(m), m);
1.2183 + }
1.2184 + extractBlossom(subblossoms[ib], base, matching);
1.2185 +
1.2186 + int en = _blossom_node_list.size();
1.2187 +
1.2188 + _blossom_potential.push_back(BlossomVariable(bn, en, pot));
1.2189 + }
1.2190 + }
1.2191 +
1.2192 + void extractMatching() {
1.2193 + std::vector<int> blossoms;
1.2194 + for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1.2195 + blossoms.push_back(c);
1.2196 + }
1.2197 +
1.2198 + for (int i = 0; i < int(blossoms.size()); ++i) {
1.2199 +
1.2200 + Value offset = (*_blossom_data)[blossoms[i]].offset;
1.2201 + (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1.2202 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1.2203 + n != INVALID; ++n) {
1.2204 + (*_node_data)[(*_node_index)[n]].pot -= offset;
1.2205 + }
1.2206 +
1.2207 + Edge matching = (*_blossom_data)[blossoms[i]].next;
1.2208 + Node base = _ugraph.source(matching);
1.2209 + extractBlossom(blossoms[i], base, matching);
1.2210 + }
1.2211 + }
1.2212 +
1.2213 + public:
1.2214 +
1.2215 + /// \brief Constructor
1.2216 + ///
1.2217 + /// Constructor.
1.2218 + MaxWeightedPerfectMatching(const UGraph& ugraph, const WeightMap& weight)
1.2219 + : _ugraph(ugraph), _weight(weight), _matching(0),
1.2220 + _node_potential(0), _blossom_potential(), _blossom_node_list(),
1.2221 + _node_num(0), _blossom_num(0),
1.2222 +
1.2223 + _blossom_index(0), _blossom_set(0), _blossom_data(0),
1.2224 + _node_index(0), _node_heap_index(0), _node_data(0),
1.2225 + _tree_set_index(0), _tree_set(0),
1.2226 +
1.2227 + _delta2_index(0), _delta2(0),
1.2228 + _delta3_index(0), _delta3(0),
1.2229 + _delta4_index(0), _delta4(0),
1.2230 +
1.2231 + _delta_sum() {}
1.2232 +
1.2233 + ~MaxWeightedPerfectMatching() {
1.2234 + destroyStructures();
1.2235 + }
1.2236 +
1.2237 + /// \name Execution control
1.2238 + /// The simplest way to execute the algorithm is to use the member
1.2239 + /// \c run() member function.
1.2240 +
1.2241 + ///@{
1.2242 +
1.2243 + /// \brief Initialize the algorithm
1.2244 + ///
1.2245 + /// Initialize the algorithm
1.2246 + void init() {
1.2247 + createStructures();
1.2248 +
1.2249 + for (EdgeIt e(_ugraph); e != INVALID; ++e) {
1.2250 + _node_heap_index->set(e, BinHeap<Value, EdgeIntMap>::PRE_HEAP);
1.2251 + }
1.2252 + for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
1.2253 + _delta3_index->set(e, _delta3->PRE_HEAP);
1.2254 + }
1.2255 + for (int i = 0; i < _blossom_num; ++i) {
1.2256 + _delta2_index->set(i, _delta2->PRE_HEAP);
1.2257 + _delta4_index->set(i, _delta4->PRE_HEAP);
1.2258 + }
1.2259 +
1.2260 + int index = 0;
1.2261 + for (NodeIt n(_ugraph); n != INVALID; ++n) {
1.2262 + Value max = std::numeric_limits<Value>::min();
1.2263 + for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
1.2264 + if (_ugraph.target(e) == n) continue;
1.2265 + if ((dualScale * _weight[e]) / 2 > max) {
1.2266 + max = (dualScale * _weight[e]) / 2;
1.2267 + }
1.2268 + }
1.2269 + _node_index->set(n, index);
1.2270 + (*_node_data)[index].pot = max;
1.2271 + int blossom =
1.2272 + _blossom_set->insert(n, std::numeric_limits<Value>::max());
1.2273 +
1.2274 + _tree_set->insert(blossom);
1.2275 +
1.2276 + (*_blossom_data)[blossom].status = EVEN;
1.2277 + (*_blossom_data)[blossom].pred = INVALID;
1.2278 + (*_blossom_data)[blossom].next = INVALID;
1.2279 + (*_blossom_data)[blossom].pot = 0;
1.2280 + (*_blossom_data)[blossom].offset = 0;
1.2281 + ++index;
1.2282 + }
1.2283 + for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
1.2284 + int si = (*_node_index)[_ugraph.source(e)];
1.2285 + int ti = (*_node_index)[_ugraph.target(e)];
1.2286 + if (_ugraph.source(e) != _ugraph.target(e)) {
1.2287 + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1.2288 + dualScale * _weight[e]) / 2);
1.2289 + }
1.2290 + }
1.2291 + }
1.2292 +
1.2293 + /// \brief Starts the algorithm
1.2294 + ///
1.2295 + /// Starts the algorithm
1.2296 + bool start() {
1.2297 + enum OpType {
1.2298 + D2, D3, D4
1.2299 + };
1.2300 +
1.2301 + int unmatched = _node_num;
1.2302 + while (unmatched > 0) {
1.2303 + Value d2 = !_delta2->empty() ?
1.2304 + _delta2->prio() : std::numeric_limits<Value>::max();
1.2305 +
1.2306 + Value d3 = !_delta3->empty() ?
1.2307 + _delta3->prio() : std::numeric_limits<Value>::max();
1.2308 +
1.2309 + Value d4 = !_delta4->empty() ?
1.2310 + _delta4->prio() : std::numeric_limits<Value>::max();
1.2311 +
1.2312 + _delta_sum = d2; OpType ot = D2;
1.2313 + if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
1.2314 + if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1.2315 +
1.2316 + if (_delta_sum == std::numeric_limits<Value>::max()) {
1.2317 + return false;
1.2318 + }
1.2319 +
1.2320 + switch (ot) {
1.2321 + case D2:
1.2322 + {
1.2323 + int blossom = _delta2->top();
1.2324 + Node n = _blossom_set->classTop(blossom);
1.2325 + Edge e = (*_node_data)[(*_node_index)[n]].heap.top();
1.2326 + extendOnEdge(e);
1.2327 + }
1.2328 + break;
1.2329 + case D3:
1.2330 + {
1.2331 + UEdge e = _delta3->top();
1.2332 +
1.2333 + int left_blossom = _blossom_set->find(_ugraph.source(e));
1.2334 + int right_blossom = _blossom_set->find(_ugraph.target(e));
1.2335 +
1.2336 + if (left_blossom == right_blossom) {
1.2337 + _delta3->pop();
1.2338 + } else {
1.2339 + int left_tree = _tree_set->find(left_blossom);
1.2340 + int right_tree = _tree_set->find(right_blossom);
1.2341 +
1.2342 + if (left_tree == right_tree) {
1.2343 + shrinkOnEdge(e, left_tree);
1.2344 + } else {
1.2345 + augmentOnEdge(e);
1.2346 + unmatched -= 2;
1.2347 + }
1.2348 + }
1.2349 + } break;
1.2350 + case D4:
1.2351 + splitBlossom(_delta4->top());
1.2352 + break;
1.2353 + }
1.2354 + }
1.2355 + extractMatching();
1.2356 + return true;
1.2357 + }
1.2358 +
1.2359 + /// \brief Runs %MaxWeightedPerfectMatching algorithm.
1.2360 + ///
1.2361 + /// This method runs the %MaxWeightedPerfectMatching algorithm.
1.2362 + ///
1.2363 + /// \note mwm.run() is just a shortcut of the following code.
1.2364 + /// \code
1.2365 + /// mwm.init();
1.2366 + /// mwm.start();
1.2367 + /// \endcode
1.2368 + bool run() {
1.2369 + init();
1.2370 + return start();
1.2371 + }
1.2372 +
1.2373 + /// @}
1.2374 +
1.2375 + /// \name Primal solution
1.2376 + /// Functions for get the primal solution, ie. the matching.
1.2377 +
1.2378 + /// @{
1.2379 +
1.2380 + /// \brief Returns the matching value.
1.2381 + ///
1.2382 + /// Returns the matching value.
1.2383 + Value matchingValue() const {
1.2384 + Value sum = 0;
1.2385 + for (NodeIt n(_ugraph); n != INVALID; ++n) {
1.2386 + if ((*_matching)[n] != INVALID) {
1.2387 + sum += _weight[(*_matching)[n]];
1.2388 + }
1.2389 + }
1.2390 + return sum /= 2;
1.2391 + }
1.2392 +
1.2393 + /// \brief Returns true when the edge is in the matching.
1.2394 + ///
1.2395 + /// Returns true when the edge is in the matching.
1.2396 + bool matching(const UEdge& edge) const {
1.2397 + return (*_matching)[_ugraph.source(edge)] == _ugraph.direct(edge, true);
1.2398 + }
1.2399 +
1.2400 + /// \brief Returns the incident matching edge.
1.2401 + ///
1.2402 + /// Returns the incident matching edge from given node.
1.2403 + Edge matching(const Node& node) const {
1.2404 + return (*_matching)[node];
1.2405 + }
1.2406 +
1.2407 + /// \brief Returns the mate of the node.
1.2408 + ///
1.2409 + /// Returns the adjancent node in a mathcing edge.
1.2410 + Node mate(const Node& node) const {
1.2411 + return _ugraph.target((*_matching)[node]);
1.2412 + }
1.2413 +
1.2414 + /// @}
1.2415 +
1.2416 + /// \name Dual solution
1.2417 + /// Functions for get the dual solution.
1.2418 +
1.2419 + /// @{
1.2420 +
1.2421 + /// \brief Returns the value of the dual solution.
1.2422 + ///
1.2423 + /// Returns the value of the dual solution. It should be equal to
1.2424 + /// the primal value scaled by \ref dualScale "dual scale".
1.2425 + Value dualValue() const {
1.2426 + Value sum = 0;
1.2427 + for (NodeIt n(_ugraph); n != INVALID; ++n) {
1.2428 + sum += nodeValue(n);
1.2429 + }
1.2430 + for (int i = 0; i < blossomNum(); ++i) {
1.2431 + sum += blossomValue(i) * (blossomSize(i) / 2);
1.2432 + }
1.2433 + return sum;
1.2434 + }
1.2435 +
1.2436 + /// \brief Returns the value of the node.
1.2437 + ///
1.2438 + /// Returns the the value of the node.
1.2439 + Value nodeValue(const Node& n) const {
1.2440 + return (*_node_potential)[n];
1.2441 + }
1.2442 +
1.2443 + /// \brief Returns the number of the blossoms in the basis.
1.2444 + ///
1.2445 + /// Returns the number of the blossoms in the basis.
1.2446 + /// \see BlossomIt
1.2447 + int blossomNum() const {
1.2448 + return _blossom_potential.size();
1.2449 + }
1.2450 +
1.2451 +
1.2452 + /// \brief Returns the number of the nodes in the blossom.
1.2453 + ///
1.2454 + /// Returns the number of the nodes in the blossom.
1.2455 + int blossomSize(int k) const {
1.2456 + return _blossom_potential[k].end - _blossom_potential[k].begin;
1.2457 + }
1.2458 +
1.2459 + /// \brief Returns the value of the blossom.
1.2460 + ///
1.2461 + /// Returns the the value of the blossom.
1.2462 + /// \see BlossomIt
1.2463 + Value blossomValue(int k) const {
1.2464 + return _blossom_potential[k].value;
1.2465 + }
1.2466 +
1.2467 + /// \brief Lemon iterator for get the items of the blossom.
1.2468 + ///
1.2469 + /// Lemon iterator for get the nodes of the blossom. This class
1.2470 + /// provides a common style lemon iterator which gives back a
1.2471 + /// subset of the nodes.
1.2472 + class BlossomIt {
1.2473 + public:
1.2474 +
1.2475 + /// \brief Constructor.
1.2476 + ///
1.2477 + /// Constructor for get the nodes of the variable.
1.2478 + BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
1.2479 + : _algorithm(&algorithm)
1.2480 + {
1.2481 + _index = _algorithm->_blossom_potential[variable].begin;
1.2482 + _last = _algorithm->_blossom_potential[variable].end;
1.2483 + }
1.2484 +
1.2485 + /// \brief Invalid constructor.
1.2486 + ///
1.2487 + /// Invalid constructor.
1.2488 + BlossomIt(Invalid) : _index(-1) {}
1.2489 +
1.2490 + /// \brief Conversion to node.
1.2491 + ///
1.2492 + /// Conversion to node.
1.2493 + operator Node() const {
1.2494 + return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
1.2495 + }
1.2496 +
1.2497 + /// \brief Increment operator.
1.2498 + ///
1.2499 + /// Increment operator.
1.2500 + BlossomIt& operator++() {
1.2501 + ++_index;
1.2502 + if (_index == _last) {
1.2503 + _index = -1;
1.2504 + }
1.2505 + return *this;
1.2506 + }
1.2507 +
1.2508 + bool operator==(const BlossomIt& it) const {
1.2509 + return _index == it._index;
1.2510 + }
1.2511 + bool operator!=(const BlossomIt& it) const {
1.2512 + return _index != it._index;
1.2513 + }
1.2514 +
1.2515 + private:
1.2516 + const MaxWeightedPerfectMatching* _algorithm;
1.2517 + int _last;
1.2518 + int _index;
1.2519 + };
1.2520 +
1.2521 + /// @}
1.2522 +
1.2523 + };
1.2524 +
1.2525
1.2526 } //END OF NAMESPACE LEMON
1.2527