1.1 --- a/lemon/matching.h Tue Dec 20 17:44:38 2011 +0100
1.2 +++ b/lemon/matching.h Tue Dec 20 18:15:14 2011 +0100
1.3 @@ -2,7 +2,7 @@
1.4 *
1.5 * This file is a part of LEMON, a generic C++ optimization library.
1.6 *
1.7 - * Copyright (C) 2003-2009
1.8 + * Copyright (C) 2003-2010
1.9 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 *
1.12 @@ -16,8 +16,8 @@
1.13 *
1.14 */
1.15
1.16 -#ifndef LEMON_MAX_MATCHING_H
1.17 -#define LEMON_MAX_MATCHING_H
1.18 +#ifndef LEMON_MATCHING_H
1.19 +#define LEMON_MATCHING_H
1.20
1.21 #include <vector>
1.22 #include <queue>
1.23 @@ -28,6 +28,7 @@
1.24 #include <lemon/unionfind.h>
1.25 #include <lemon/bin_heap.h>
1.26 #include <lemon/maps.h>
1.27 +#include <lemon/fractional_matching.h>
1.28
1.29 ///\ingroup matching
1.30 ///\file
1.31 @@ -41,7 +42,7 @@
1.32 ///
1.33 /// This class implements Edmonds' alternating forest matching algorithm
1.34 /// for finding a maximum cardinality matching in a general undirected graph.
1.35 - /// It can be started from an arbitrary initial matching
1.36 + /// It can be started from an arbitrary initial matching
1.37 /// (the default is the empty one).
1.38 ///
1.39 /// The dual solution of the problem is a map of the nodes to
1.40 @@ -69,11 +70,11 @@
1.41
1.42 ///\brief Status constants for Gallai-Edmonds decomposition.
1.43 ///
1.44 - ///These constants are used for indicating the Gallai-Edmonds
1.45 + ///These constants are used for indicating the Gallai-Edmonds
1.46 ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
1.47 ///induce a subgraph with factor-critical components, the nodes with
1.48 ///status \c ODD (or \c A) form the canonical barrier, and the nodes
1.49 - ///with status \c MATCHED (or \c C) induce a subgraph having a
1.50 + ///with status \c MATCHED (or \c C) induce a subgraph having a
1.51 ///perfect matching.
1.52 enum Status {
1.53 EVEN = 1, ///< = 1. (\c D is an alias for \c EVEN.)
1.54 @@ -512,7 +513,7 @@
1.55 }
1.56 }
1.57
1.58 - /// \brief Start Edmonds' algorithm with a heuristic improvement
1.59 + /// \brief Start Edmonds' algorithm with a heuristic improvement
1.60 /// for dense graphs
1.61 ///
1.62 /// This function runs Edmonds' algorithm with a heuristic of postponing
1.63 @@ -534,8 +535,8 @@
1.64
1.65 /// \brief Run Edmonds' algorithm
1.66 ///
1.67 - /// This function runs Edmonds' algorithm. An additional heuristic of
1.68 - /// postponing shrinks is used for relatively dense graphs
1.69 + /// This function runs Edmonds' algorithm. An additional heuristic of
1.70 + /// postponing shrinks is used for relatively dense graphs
1.71 /// (for which <tt>m>=2*n</tt> holds).
1.72 void run() {
1.73 if (countEdges(_graph) < 2 * countNodes(_graph)) {
1.74 @@ -556,7 +557,7 @@
1.75
1.76 /// \brief Return the size (cardinality) of the matching.
1.77 ///
1.78 - /// This function returns the size (cardinality) of the current matching.
1.79 + /// This function returns the size (cardinality) of the current matching.
1.80 /// After run() it returns the size of the maximum matching in the graph.
1.81 int matchingSize() const {
1.82 int size = 0;
1.83 @@ -570,7 +571,7 @@
1.84
1.85 /// \brief Return \c true if the given edge is in the matching.
1.86 ///
1.87 - /// This function returns \c true if the given edge is in the current
1.88 + /// This function returns \c true if the given edge is in the current
1.89 /// matching.
1.90 bool matching(const Edge& edge) const {
1.91 return edge == (*_matching)[_graph.u(edge)];
1.92 @@ -579,7 +580,7 @@
1.93 /// \brief Return the matching arc (or edge) incident to the given node.
1.94 ///
1.95 /// This function returns the matching arc (or edge) incident to the
1.96 - /// given node in the current matching or \c INVALID if the node is
1.97 + /// given node in the current matching or \c INVALID if the node is
1.98 /// not covered by the matching.
1.99 Arc matching(const Node& n) const {
1.100 return (*_matching)[n];
1.101 @@ -595,7 +596,7 @@
1.102
1.103 /// \brief Return the mate of the given node.
1.104 ///
1.105 - /// This function returns the mate of the given node in the current
1.106 + /// This function returns the mate of the given node in the current
1.107 /// matching or \c INVALID if the node is not covered by the matching.
1.108 Node mate(const Node& n) const {
1.109 return (*_matching)[n] != INVALID ?
1.110 @@ -605,7 +606,7 @@
1.111 /// @}
1.112
1.113 /// \name Dual Solution
1.114 - /// Functions to get the dual solution, i.e. the Gallai-Edmonds
1.115 + /// Functions to get the dual solution, i.e. the Gallai-Edmonds
1.116 /// decomposition.
1.117
1.118 /// @{
1.119 @@ -648,8 +649,8 @@
1.120 /// on extensive use of priority queues and provides
1.121 /// \f$O(nm\log n)\f$ time complexity.
1.122 ///
1.123 - /// The maximum weighted matching problem is to find a subset of the
1.124 - /// edges in an undirected graph with maximum overall weight for which
1.125 + /// The maximum weighted matching problem is to find a subset of the
1.126 + /// edges in an undirected graph with maximum overall weight for which
1.127 /// each node has at most one incident edge.
1.128 /// It can be formulated with the following linear program.
1.129 /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
1.130 @@ -673,16 +674,16 @@
1.131 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
1.132 \frac{\vert B \vert - 1}{2}z_B\f] */
1.133 ///
1.134 - /// The algorithm can be executed with the run() function.
1.135 + /// The algorithm can be executed with the run() function.
1.136 /// After it the matching (the primal solution) and the dual solution
1.137 - /// can be obtained using the query functions and the
1.138 - /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
1.139 - /// which is able to iterate on the nodes of a blossom.
1.140 + /// can be obtained using the query functions and the
1.141 + /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
1.142 + /// which is able to iterate on the nodes of a blossom.
1.143 /// If the value type is integer, then the dual solution is multiplied
1.144 /// by \ref MaxWeightedMatching::dualScale "4".
1.145 ///
1.146 /// \tparam GR The undirected graph type the algorithm runs on.
1.147 - /// \tparam WM The type edge weight map. The default type is
1.148 + /// \tparam WM The type edge weight map. The default type is
1.149 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
1.150 #ifdef DOXYGEN
1.151 template <typename GR, typename WM>
1.152 @@ -745,7 +746,7 @@
1.153 typedef RangeMap<int> IntIntMap;
1.154
1.155 enum Status {
1.156 - EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
1.157 + EVEN = -1, MATCHED = 0, ODD = 1
1.158 };
1.159
1.160 typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
1.161 @@ -797,6 +798,10 @@
1.162 BinHeap<Value, IntIntMap> *_delta4;
1.163
1.164 Value _delta_sum;
1.165 + int _unmatched;
1.166 +
1.167 + typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
1.168 + FractionalMatching *_fractional;
1.169
1.170 void createStructures() {
1.171 _node_num = countNodes(_graph);
1.172 @@ -863,9 +868,6 @@
1.173 }
1.174
1.175 void destroyStructures() {
1.176 - _node_num = countNodes(_graph);
1.177 - _blossom_num = _node_num * 3 / 2;
1.178 -
1.179 if (_matching) {
1.180 delete _matching;
1.181 }
1.182 @@ -941,10 +943,6 @@
1.183 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1.184 _delta3->push(e, rw / 2);
1.185 }
1.186 - } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1.187 - if (_delta3->state(e) != _delta3->IN_HEAP) {
1.188 - _delta3->push(e, rw);
1.189 - }
1.190 } else {
1.191 typename std::map<int, Arc>::iterator it =
1.192 (*_node_data)[vi].heap_index.find(tree);
1.193 @@ -968,202 +966,6 @@
1.194 _delta2->push(vb, _blossom_set->classPrio(vb) -
1.195 (*_blossom_data)[vb].offset);
1.196 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1.197 - (*_blossom_data)[vb].offset){
1.198 - _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1.199 - (*_blossom_data)[vb].offset);
1.200 - }
1.201 - }
1.202 - }
1.203 - }
1.204 - }
1.205 - }
1.206 - (*_blossom_data)[blossom].offset = 0;
1.207 - }
1.208 -
1.209 - void matchedToOdd(int blossom) {
1.210 - if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.211 - _delta2->erase(blossom);
1.212 - }
1.213 - (*_blossom_data)[blossom].offset += _delta_sum;
1.214 - if (!_blossom_set->trivial(blossom)) {
1.215 - _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
1.216 - (*_blossom_data)[blossom].offset);
1.217 - }
1.218 - }
1.219 -
1.220 - void evenToMatched(int blossom, int tree) {
1.221 - if (!_blossom_set->trivial(blossom)) {
1.222 - (*_blossom_data)[blossom].pot += 2 * _delta_sum;
1.223 - }
1.224 -
1.225 - for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.226 - n != INVALID; ++n) {
1.227 - int ni = (*_node_index)[n];
1.228 - (*_node_data)[ni].pot -= _delta_sum;
1.229 -
1.230 - _delta1->erase(n);
1.231 -
1.232 - for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.233 - Node v = _graph.source(e);
1.234 - int vb = _blossom_set->find(v);
1.235 - int vi = (*_node_index)[v];
1.236 -
1.237 - Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.238 - dualScale * _weight[e];
1.239 -
1.240 - if (vb == blossom) {
1.241 - if (_delta3->state(e) == _delta3->IN_HEAP) {
1.242 - _delta3->erase(e);
1.243 - }
1.244 - } else if ((*_blossom_data)[vb].status == EVEN) {
1.245 -
1.246 - if (_delta3->state(e) == _delta3->IN_HEAP) {
1.247 - _delta3->erase(e);
1.248 - }
1.249 -
1.250 - int vt = _tree_set->find(vb);
1.251 -
1.252 - if (vt != tree) {
1.253 -
1.254 - Arc r = _graph.oppositeArc(e);
1.255 -
1.256 - typename std::map<int, Arc>::iterator it =
1.257 - (*_node_data)[ni].heap_index.find(vt);
1.258 -
1.259 - if (it != (*_node_data)[ni].heap_index.end()) {
1.260 - if ((*_node_data)[ni].heap[it->second] > rw) {
1.261 - (*_node_data)[ni].heap.replace(it->second, r);
1.262 - (*_node_data)[ni].heap.decrease(r, rw);
1.263 - it->second = r;
1.264 - }
1.265 - } else {
1.266 - (*_node_data)[ni].heap.push(r, rw);
1.267 - (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1.268 - }
1.269 -
1.270 - if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1.271 - _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1.272 -
1.273 - if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1.274 - _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.275 - (*_blossom_data)[blossom].offset);
1.276 - } else if ((*_delta2)[blossom] >
1.277 - _blossom_set->classPrio(blossom) -
1.278 - (*_blossom_data)[blossom].offset){
1.279 - _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1.280 - (*_blossom_data)[blossom].offset);
1.281 - }
1.282 - }
1.283 - }
1.284 -
1.285 - } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1.286 - if (_delta3->state(e) == _delta3->IN_HEAP) {
1.287 - _delta3->erase(e);
1.288 - }
1.289 - } else {
1.290 -
1.291 - typename std::map<int, Arc>::iterator it =
1.292 - (*_node_data)[vi].heap_index.find(tree);
1.293 -
1.294 - if (it != (*_node_data)[vi].heap_index.end()) {
1.295 - (*_node_data)[vi].heap.erase(it->second);
1.296 - (*_node_data)[vi].heap_index.erase(it);
1.297 - if ((*_node_data)[vi].heap.empty()) {
1.298 - _blossom_set->increase(v, std::numeric_limits<Value>::max());
1.299 - } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1.300 - _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1.301 - }
1.302 -
1.303 - if ((*_blossom_data)[vb].status == MATCHED) {
1.304 - if (_blossom_set->classPrio(vb) ==
1.305 - std::numeric_limits<Value>::max()) {
1.306 - _delta2->erase(vb);
1.307 - } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1.308 - (*_blossom_data)[vb].offset) {
1.309 - _delta2->increase(vb, _blossom_set->classPrio(vb) -
1.310 - (*_blossom_data)[vb].offset);
1.311 - }
1.312 - }
1.313 - }
1.314 - }
1.315 - }
1.316 - }
1.317 - }
1.318 -
1.319 - void oddToMatched(int blossom) {
1.320 - (*_blossom_data)[blossom].offset -= _delta_sum;
1.321 -
1.322 - if (_blossom_set->classPrio(blossom) !=
1.323 - std::numeric_limits<Value>::max()) {
1.324 - _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.325 - (*_blossom_data)[blossom].offset);
1.326 - }
1.327 -
1.328 - if (!_blossom_set->trivial(blossom)) {
1.329 - _delta4->erase(blossom);
1.330 - }
1.331 - }
1.332 -
1.333 - void oddToEven(int blossom, int tree) {
1.334 - if (!_blossom_set->trivial(blossom)) {
1.335 - _delta4->erase(blossom);
1.336 - (*_blossom_data)[blossom].pot -=
1.337 - 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1.338 - }
1.339 -
1.340 - for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.341 - n != INVALID; ++n) {
1.342 - int ni = (*_node_index)[n];
1.343 -
1.344 - _blossom_set->increase(n, std::numeric_limits<Value>::max());
1.345 -
1.346 - (*_node_data)[ni].heap.clear();
1.347 - (*_node_data)[ni].heap_index.clear();
1.348 - (*_node_data)[ni].pot +=
1.349 - 2 * _delta_sum - (*_blossom_data)[blossom].offset;
1.350 -
1.351 - _delta1->push(n, (*_node_data)[ni].pot);
1.352 -
1.353 - for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.354 - Node v = _graph.source(e);
1.355 - int vb = _blossom_set->find(v);
1.356 - int vi = (*_node_index)[v];
1.357 -
1.358 - Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.359 - dualScale * _weight[e];
1.360 -
1.361 - if ((*_blossom_data)[vb].status == EVEN) {
1.362 - if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1.363 - _delta3->push(e, rw / 2);
1.364 - }
1.365 - } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1.366 - if (_delta3->state(e) != _delta3->IN_HEAP) {
1.367 - _delta3->push(e, rw);
1.368 - }
1.369 - } else {
1.370 -
1.371 - typename std::map<int, Arc>::iterator it =
1.372 - (*_node_data)[vi].heap_index.find(tree);
1.373 -
1.374 - if (it != (*_node_data)[vi].heap_index.end()) {
1.375 - if ((*_node_data)[vi].heap[it->second] > rw) {
1.376 - (*_node_data)[vi].heap.replace(it->second, e);
1.377 - (*_node_data)[vi].heap.decrease(e, rw);
1.378 - it->second = e;
1.379 - }
1.380 - } else {
1.381 - (*_node_data)[vi].heap.push(e, rw);
1.382 - (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1.383 - }
1.384 -
1.385 - if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1.386 - _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1.387 -
1.388 - if ((*_blossom_data)[vb].status == MATCHED) {
1.389 - if (_delta2->state(vb) != _delta2->IN_HEAP) {
1.390 - _delta2->push(vb, _blossom_set->classPrio(vb) -
1.391 - (*_blossom_data)[vb].offset);
1.392 - } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1.393 (*_blossom_data)[vb].offset) {
1.394 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1.395 (*_blossom_data)[vb].offset);
1.396 @@ -1176,43 +978,145 @@
1.397 (*_blossom_data)[blossom].offset = 0;
1.398 }
1.399
1.400 -
1.401 - void matchedToUnmatched(int blossom) {
1.402 + void matchedToOdd(int blossom) {
1.403 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.404 _delta2->erase(blossom);
1.405 }
1.406 + (*_blossom_data)[blossom].offset += _delta_sum;
1.407 + if (!_blossom_set->trivial(blossom)) {
1.408 + _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
1.409 + (*_blossom_data)[blossom].offset);
1.410 + }
1.411 + }
1.412 +
1.413 + void evenToMatched(int blossom, int tree) {
1.414 + if (!_blossom_set->trivial(blossom)) {
1.415 + (*_blossom_data)[blossom].pot += 2 * _delta_sum;
1.416 + }
1.417
1.418 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.419 n != INVALID; ++n) {
1.420 int ni = (*_node_index)[n];
1.421 -
1.422 - _blossom_set->increase(n, std::numeric_limits<Value>::max());
1.423 -
1.424 - (*_node_data)[ni].heap.clear();
1.425 - (*_node_data)[ni].heap_index.clear();
1.426 -
1.427 - for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.428 - Node v = _graph.target(e);
1.429 + (*_node_data)[ni].pot -= _delta_sum;
1.430 +
1.431 + _delta1->erase(n);
1.432 +
1.433 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.434 + Node v = _graph.source(e);
1.435 int vb = _blossom_set->find(v);
1.436 int vi = (*_node_index)[v];
1.437
1.438 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.439 dualScale * _weight[e];
1.440
1.441 - if ((*_blossom_data)[vb].status == EVEN) {
1.442 - if (_delta3->state(e) != _delta3->IN_HEAP) {
1.443 - _delta3->push(e, rw);
1.444 + if (vb == blossom) {
1.445 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.446 + _delta3->erase(e);
1.447 + }
1.448 + } else if ((*_blossom_data)[vb].status == EVEN) {
1.449 +
1.450 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.451 + _delta3->erase(e);
1.452 + }
1.453 +
1.454 + int vt = _tree_set->find(vb);
1.455 +
1.456 + if (vt != tree) {
1.457 +
1.458 + Arc r = _graph.oppositeArc(e);
1.459 +
1.460 + typename std::map<int, Arc>::iterator it =
1.461 + (*_node_data)[ni].heap_index.find(vt);
1.462 +
1.463 + if (it != (*_node_data)[ni].heap_index.end()) {
1.464 + if ((*_node_data)[ni].heap[it->second] > rw) {
1.465 + (*_node_data)[ni].heap.replace(it->second, r);
1.466 + (*_node_data)[ni].heap.decrease(r, rw);
1.467 + it->second = r;
1.468 + }
1.469 + } else {
1.470 + (*_node_data)[ni].heap.push(r, rw);
1.471 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1.472 + }
1.473 +
1.474 + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1.475 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1.476 +
1.477 + if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1.478 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.479 + (*_blossom_data)[blossom].offset);
1.480 + } else if ((*_delta2)[blossom] >
1.481 + _blossom_set->classPrio(blossom) -
1.482 + (*_blossom_data)[blossom].offset){
1.483 + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1.484 + (*_blossom_data)[blossom].offset);
1.485 + }
1.486 + }
1.487 + }
1.488 + } else {
1.489 +
1.490 + typename std::map<int, Arc>::iterator it =
1.491 + (*_node_data)[vi].heap_index.find(tree);
1.492 +
1.493 + if (it != (*_node_data)[vi].heap_index.end()) {
1.494 + (*_node_data)[vi].heap.erase(it->second);
1.495 + (*_node_data)[vi].heap_index.erase(it);
1.496 + if ((*_node_data)[vi].heap.empty()) {
1.497 + _blossom_set->increase(v, std::numeric_limits<Value>::max());
1.498 + } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1.499 + _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1.500 + }
1.501 +
1.502 + if ((*_blossom_data)[vb].status == MATCHED) {
1.503 + if (_blossom_set->classPrio(vb) ==
1.504 + std::numeric_limits<Value>::max()) {
1.505 + _delta2->erase(vb);
1.506 + } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1.507 + (*_blossom_data)[vb].offset) {
1.508 + _delta2->increase(vb, _blossom_set->classPrio(vb) -
1.509 + (*_blossom_data)[vb].offset);
1.510 + }
1.511 + }
1.512 }
1.513 }
1.514 }
1.515 }
1.516 }
1.517
1.518 - void unmatchedToMatched(int blossom) {
1.519 + void oddToMatched(int blossom) {
1.520 + (*_blossom_data)[blossom].offset -= _delta_sum;
1.521 +
1.522 + if (_blossom_set->classPrio(blossom) !=
1.523 + std::numeric_limits<Value>::max()) {
1.524 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.525 + (*_blossom_data)[blossom].offset);
1.526 + }
1.527 +
1.528 + if (!_blossom_set->trivial(blossom)) {
1.529 + _delta4->erase(blossom);
1.530 + }
1.531 + }
1.532 +
1.533 + void oddToEven(int blossom, int tree) {
1.534 + if (!_blossom_set->trivial(blossom)) {
1.535 + _delta4->erase(blossom);
1.536 + (*_blossom_data)[blossom].pot -=
1.537 + 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1.538 + }
1.539 +
1.540 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.541 n != INVALID; ++n) {
1.542 int ni = (*_node_index)[n];
1.543
1.544 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
1.545 +
1.546 + (*_node_data)[ni].heap.clear();
1.547 + (*_node_data)[ni].heap_index.clear();
1.548 + (*_node_data)[ni].pot +=
1.549 + 2 * _delta_sum - (*_blossom_data)[blossom].offset;
1.550 +
1.551 + _delta1->push(n, (*_node_data)[ni].pot);
1.552 +
1.553 for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.554 Node v = _graph.source(e);
1.555 int vb = _blossom_set->find(v);
1.556 @@ -1221,54 +1125,44 @@
1.557 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.558 dualScale * _weight[e];
1.559
1.560 - if (vb == blossom) {
1.561 - if (_delta3->state(e) == _delta3->IN_HEAP) {
1.562 - _delta3->erase(e);
1.563 + if ((*_blossom_data)[vb].status == EVEN) {
1.564 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1.565 + _delta3->push(e, rw / 2);
1.566 }
1.567 - } else if ((*_blossom_data)[vb].status == EVEN) {
1.568 -
1.569 - if (_delta3->state(e) == _delta3->IN_HEAP) {
1.570 - _delta3->erase(e);
1.571 - }
1.572 -
1.573 - int vt = _tree_set->find(vb);
1.574 -
1.575 - Arc r = _graph.oppositeArc(e);
1.576 + } else {
1.577
1.578 typename std::map<int, Arc>::iterator it =
1.579 - (*_node_data)[ni].heap_index.find(vt);
1.580 -
1.581 - if (it != (*_node_data)[ni].heap_index.end()) {
1.582 - if ((*_node_data)[ni].heap[it->second] > rw) {
1.583 - (*_node_data)[ni].heap.replace(it->second, r);
1.584 - (*_node_data)[ni].heap.decrease(r, rw);
1.585 - it->second = r;
1.586 + (*_node_data)[vi].heap_index.find(tree);
1.587 +
1.588 + if (it != (*_node_data)[vi].heap_index.end()) {
1.589 + if ((*_node_data)[vi].heap[it->second] > rw) {
1.590 + (*_node_data)[vi].heap.replace(it->second, e);
1.591 + (*_node_data)[vi].heap.decrease(e, rw);
1.592 + it->second = e;
1.593 }
1.594 } else {
1.595 - (*_node_data)[ni].heap.push(r, rw);
1.596 - (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1.597 + (*_node_data)[vi].heap.push(e, rw);
1.598 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1.599 }
1.600
1.601 - if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1.602 - _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1.603 -
1.604 - if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1.605 - _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.606 - (*_blossom_data)[blossom].offset);
1.607 - } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
1.608 - (*_blossom_data)[blossom].offset){
1.609 - _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1.610 - (*_blossom_data)[blossom].offset);
1.611 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1.612 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1.613 +
1.614 + if ((*_blossom_data)[vb].status == MATCHED) {
1.615 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
1.616 + _delta2->push(vb, _blossom_set->classPrio(vb) -
1.617 + (*_blossom_data)[vb].offset);
1.618 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1.619 + (*_blossom_data)[vb].offset) {
1.620 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1.621 + (*_blossom_data)[vb].offset);
1.622 + }
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 + (*_blossom_data)[blossom].offset = 0;
1.634 }
1.635
1.636 void alternatePath(int even, int tree) {
1.637 @@ -1313,39 +1207,42 @@
1.638 alternatePath(blossom, tree);
1.639 destroyTree(tree);
1.640
1.641 - (*_blossom_data)[blossom].status = UNMATCHED;
1.642 (*_blossom_data)[blossom].base = node;
1.643 - matchedToUnmatched(blossom);
1.644 + (*_blossom_data)[blossom].next = INVALID;
1.645 }
1.646
1.647 -
1.648 void augmentOnEdge(const Edge& edge) {
1.649
1.650 int left = _blossom_set->find(_graph.u(edge));
1.651 int right = _blossom_set->find(_graph.v(edge));
1.652
1.653 - if ((*_blossom_data)[left].status == EVEN) {
1.654 - int left_tree = _tree_set->find(left);
1.655 - alternatePath(left, left_tree);
1.656 - destroyTree(left_tree);
1.657 - } else {
1.658 - (*_blossom_data)[left].status = MATCHED;
1.659 - unmatchedToMatched(left);
1.660 - }
1.661 -
1.662 - if ((*_blossom_data)[right].status == EVEN) {
1.663 - int right_tree = _tree_set->find(right);
1.664 - alternatePath(right, right_tree);
1.665 - destroyTree(right_tree);
1.666 - } else {
1.667 - (*_blossom_data)[right].status = MATCHED;
1.668 - unmatchedToMatched(right);
1.669 - }
1.670 + int left_tree = _tree_set->find(left);
1.671 + alternatePath(left, left_tree);
1.672 + destroyTree(left_tree);
1.673 +
1.674 + int right_tree = _tree_set->find(right);
1.675 + alternatePath(right, right_tree);
1.676 + destroyTree(right_tree);
1.677
1.678 (*_blossom_data)[left].next = _graph.direct(edge, true);
1.679 (*_blossom_data)[right].next = _graph.direct(edge, false);
1.680 }
1.681
1.682 + void augmentOnArc(const Arc& arc) {
1.683 +
1.684 + int left = _blossom_set->find(_graph.source(arc));
1.685 + int right = _blossom_set->find(_graph.target(arc));
1.686 +
1.687 + (*_blossom_data)[left].status = MATCHED;
1.688 +
1.689 + int right_tree = _tree_set->find(right);
1.690 + alternatePath(right, right_tree);
1.691 + destroyTree(right_tree);
1.692 +
1.693 + (*_blossom_data)[left].next = arc;
1.694 + (*_blossom_data)[right].next = _graph.oppositeArc(arc);
1.695 + }
1.696 +
1.697 void extendOnArc(const Arc& arc) {
1.698 int base = _blossom_set->find(_graph.target(arc));
1.699 int tree = _tree_set->find(base);
1.700 @@ -1548,7 +1445,7 @@
1.701 _tree_set->insert(sb, tree);
1.702 (*_blossom_data)[sb].pred = pred;
1.703 (*_blossom_data)[sb].next =
1.704 - _graph.oppositeArc((*_blossom_data)[tb].next);
1.705 + _graph.oppositeArc((*_blossom_data)[tb].next);
1.706
1.707 pred = (*_blossom_data)[ub].next;
1.708
1.709 @@ -1648,7 +1545,7 @@
1.710 }
1.711
1.712 for (int i = 0; i < int(blossoms.size()); ++i) {
1.713 - if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
1.714 + if ((*_blossom_data)[blossoms[i]].next != INVALID) {
1.715
1.716 Value offset = (*_blossom_data)[blossoms[i]].offset;
1.717 (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1.718 @@ -1686,10 +1583,16 @@
1.719 _delta3_index(0), _delta3(0),
1.720 _delta4_index(0), _delta4(0),
1.721
1.722 - _delta_sum() {}
1.723 + _delta_sum(), _unmatched(0),
1.724 +
1.725 + _fractional(0)
1.726 + {}
1.727
1.728 ~MaxWeightedMatching() {
1.729 destroyStructures();
1.730 + if (_fractional) {
1.731 + delete _fractional;
1.732 + }
1.733 }
1.734
1.735 /// \name Execution Control
1.736 @@ -1720,7 +1623,9 @@
1.737 (*_delta2_index)[i] = _delta2->PRE_HEAP;
1.738 (*_delta4_index)[i] = _delta4->PRE_HEAP;
1.739 }
1.740 -
1.741 +
1.742 + _unmatched = _node_num;
1.743 +
1.744 _delta1->clear();
1.745 _delta2->clear();
1.746 _delta3->clear();
1.747 @@ -1764,18 +1669,167 @@
1.748 }
1.749 }
1.750
1.751 + /// \brief Initialize the algorithm with fractional matching
1.752 + ///
1.753 + /// This function initializes the algorithm with a fractional
1.754 + /// matching. This initialization is also called jumpstart heuristic.
1.755 + void fractionalInit() {
1.756 + createStructures();
1.757 +
1.758 + _blossom_node_list.clear();
1.759 + _blossom_potential.clear();
1.760 +
1.761 + if (_fractional == 0) {
1.762 + _fractional = new FractionalMatching(_graph, _weight, false);
1.763 + }
1.764 + _fractional->run();
1.765 +
1.766 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.767 + (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1.768 + }
1.769 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.770 + (*_delta1_index)[n] = _delta1->PRE_HEAP;
1.771 + }
1.772 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.773 + (*_delta3_index)[e] = _delta3->PRE_HEAP;
1.774 + }
1.775 + for (int i = 0; i < _blossom_num; ++i) {
1.776 + (*_delta2_index)[i] = _delta2->PRE_HEAP;
1.777 + (*_delta4_index)[i] = _delta4->PRE_HEAP;
1.778 + }
1.779 +
1.780 + _unmatched = 0;
1.781 +
1.782 + _delta1->clear();
1.783 + _delta2->clear();
1.784 + _delta3->clear();
1.785 + _delta4->clear();
1.786 + _blossom_set->clear();
1.787 + _tree_set->clear();
1.788 +
1.789 + int index = 0;
1.790 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.791 + Value pot = _fractional->nodeValue(n);
1.792 + (*_node_index)[n] = index;
1.793 + (*_node_data)[index].pot = pot;
1.794 + (*_node_data)[index].heap_index.clear();
1.795 + (*_node_data)[index].heap.clear();
1.796 + int blossom =
1.797 + _blossom_set->insert(n, std::numeric_limits<Value>::max());
1.798 +
1.799 + (*_blossom_data)[blossom].status = MATCHED;
1.800 + (*_blossom_data)[blossom].pred = INVALID;
1.801 + (*_blossom_data)[blossom].next = _fractional->matching(n);
1.802 + if (_fractional->matching(n) == INVALID) {
1.803 + (*_blossom_data)[blossom].base = n;
1.804 + }
1.805 + (*_blossom_data)[blossom].pot = 0;
1.806 + (*_blossom_data)[blossom].offset = 0;
1.807 + ++index;
1.808 + }
1.809 +
1.810 + typename Graph::template NodeMap<bool> processed(_graph, false);
1.811 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.812 + if (processed[n]) continue;
1.813 + processed[n] = true;
1.814 + if (_fractional->matching(n) == INVALID) continue;
1.815 + int num = 1;
1.816 + Node v = _graph.target(_fractional->matching(n));
1.817 + while (n != v) {
1.818 + processed[v] = true;
1.819 + v = _graph.target(_fractional->matching(v));
1.820 + ++num;
1.821 + }
1.822 +
1.823 + if (num % 2 == 1) {
1.824 + std::vector<int> subblossoms(num);
1.825 +
1.826 + subblossoms[--num] = _blossom_set->find(n);
1.827 + _delta1->push(n, _fractional->nodeValue(n));
1.828 + v = _graph.target(_fractional->matching(n));
1.829 + while (n != v) {
1.830 + subblossoms[--num] = _blossom_set->find(v);
1.831 + _delta1->push(v, _fractional->nodeValue(v));
1.832 + v = _graph.target(_fractional->matching(v));
1.833 + }
1.834 +
1.835 + int surface =
1.836 + _blossom_set->join(subblossoms.begin(), subblossoms.end());
1.837 + (*_blossom_data)[surface].status = EVEN;
1.838 + (*_blossom_data)[surface].pred = INVALID;
1.839 + (*_blossom_data)[surface].next = INVALID;
1.840 + (*_blossom_data)[surface].pot = 0;
1.841 + (*_blossom_data)[surface].offset = 0;
1.842 +
1.843 + _tree_set->insert(surface);
1.844 + ++_unmatched;
1.845 + }
1.846 + }
1.847 +
1.848 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.849 + int si = (*_node_index)[_graph.u(e)];
1.850 + int sb = _blossom_set->find(_graph.u(e));
1.851 + int ti = (*_node_index)[_graph.v(e)];
1.852 + int tb = _blossom_set->find(_graph.v(e));
1.853 + if ((*_blossom_data)[sb].status == EVEN &&
1.854 + (*_blossom_data)[tb].status == EVEN && sb != tb) {
1.855 + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1.856 + dualScale * _weight[e]) / 2);
1.857 + }
1.858 + }
1.859 +
1.860 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.861 + int nb = _blossom_set->find(n);
1.862 + if ((*_blossom_data)[nb].status != MATCHED) continue;
1.863 + int ni = (*_node_index)[n];
1.864 +
1.865 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.866 + Node v = _graph.target(e);
1.867 + int vb = _blossom_set->find(v);
1.868 + int vi = (*_node_index)[v];
1.869 +
1.870 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.871 + dualScale * _weight[e];
1.872 +
1.873 + if ((*_blossom_data)[vb].status == EVEN) {
1.874 +
1.875 + int vt = _tree_set->find(vb);
1.876 +
1.877 + typename std::map<int, Arc>::iterator it =
1.878 + (*_node_data)[ni].heap_index.find(vt);
1.879 +
1.880 + if (it != (*_node_data)[ni].heap_index.end()) {
1.881 + if ((*_node_data)[ni].heap[it->second] > rw) {
1.882 + (*_node_data)[ni].heap.replace(it->second, e);
1.883 + (*_node_data)[ni].heap.decrease(e, rw);
1.884 + it->second = e;
1.885 + }
1.886 + } else {
1.887 + (*_node_data)[ni].heap.push(e, rw);
1.888 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
1.889 + }
1.890 + }
1.891 + }
1.892 +
1.893 + if (!(*_node_data)[ni].heap.empty()) {
1.894 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1.895 + _delta2->push(nb, _blossom_set->classPrio(nb));
1.896 + }
1.897 + }
1.898 + }
1.899 +
1.900 /// \brief Start the algorithm
1.901 ///
1.902 /// This function starts the algorithm.
1.903 ///
1.904 - /// \pre \ref init() must be called before using this function.
1.905 + /// \pre \ref init() or \ref fractionalInit() must be called
1.906 + /// before using this function.
1.907 void start() {
1.908 enum OpType {
1.909 D1, D2, D3, D4
1.910 };
1.911
1.912 - int unmatched = _node_num;
1.913 - while (unmatched > 0) {
1.914 + while (_unmatched > 0) {
1.915 Value d1 = !_delta1->empty() ?
1.916 _delta1->prio() : std::numeric_limits<Value>::max();
1.917
1.918 @@ -1788,26 +1842,30 @@
1.919 Value d4 = !_delta4->empty() ?
1.920 _delta4->prio() : std::numeric_limits<Value>::max();
1.921
1.922 - _delta_sum = d1; OpType ot = D1;
1.923 + _delta_sum = d3; OpType ot = D3;
1.924 + if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
1.925 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1.926 - if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
1.927 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1.928
1.929 -
1.930 switch (ot) {
1.931 case D1:
1.932 {
1.933 Node n = _delta1->top();
1.934 unmatchNode(n);
1.935 - --unmatched;
1.936 + --_unmatched;
1.937 }
1.938 break;
1.939 case D2:
1.940 {
1.941 int blossom = _delta2->top();
1.942 Node n = _blossom_set->classTop(blossom);
1.943 - Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
1.944 - extendOnArc(e);
1.945 + Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
1.946 + if ((*_blossom_data)[blossom].next == INVALID) {
1.947 + augmentOnArc(a);
1.948 + --_unmatched;
1.949 + } else {
1.950 + extendOnArc(a);
1.951 + }
1.952 }
1.953 break;
1.954 case D3:
1.955 @@ -1820,26 +1878,14 @@
1.956 if (left_blossom == right_blossom) {
1.957 _delta3->pop();
1.958 } else {
1.959 - int left_tree;
1.960 - if ((*_blossom_data)[left_blossom].status == EVEN) {
1.961 - left_tree = _tree_set->find(left_blossom);
1.962 - } else {
1.963 - left_tree = -1;
1.964 - ++unmatched;
1.965 - }
1.966 - int right_tree;
1.967 - if ((*_blossom_data)[right_blossom].status == EVEN) {
1.968 - right_tree = _tree_set->find(right_blossom);
1.969 - } else {
1.970 - right_tree = -1;
1.971 - ++unmatched;
1.972 - }
1.973 + int left_tree = _tree_set->find(left_blossom);
1.974 + int right_tree = _tree_set->find(right_blossom);
1.975
1.976 if (left_tree == right_tree) {
1.977 shrinkOnEdge(e, left_tree);
1.978 } else {
1.979 augmentOnEdge(e);
1.980 - unmatched -= 2;
1.981 + _unmatched -= 2;
1.982 }
1.983 }
1.984 } break;
1.985 @@ -1857,18 +1903,18 @@
1.986 ///
1.987 /// \note mwm.run() is just a shortcut of the following code.
1.988 /// \code
1.989 - /// mwm.init();
1.990 + /// mwm.fractionalInit();
1.991 /// mwm.start();
1.992 /// \endcode
1.993 void run() {
1.994 - init();
1.995 + fractionalInit();
1.996 start();
1.997 }
1.998
1.999 /// @}
1.1000
1.1001 /// \name Primal Solution
1.1002 - /// Functions to get the primal solution, i.e. the maximum weighted
1.1003 + /// Functions to get the primal solution, i.e. the maximum weighted
1.1004 /// matching.\n
1.1005 /// Either \ref run() or \ref start() function should be called before
1.1006 /// using them.
1.1007 @@ -1887,7 +1933,7 @@
1.1008 sum += _weight[(*_matching)[n]];
1.1009 }
1.1010 }
1.1011 - return sum /= 2;
1.1012 + return sum / 2;
1.1013 }
1.1014
1.1015 /// \brief Return the size (cardinality) of the matching.
1.1016 @@ -1907,7 +1953,7 @@
1.1017
1.1018 /// \brief Return \c true if the given edge is in the matching.
1.1019 ///
1.1020 - /// This function returns \c true if the given edge is in the found
1.1021 + /// This function returns \c true if the given edge is in the found
1.1022 /// matching.
1.1023 ///
1.1024 /// \pre Either run() or start() must be called before using this function.
1.1025 @@ -1918,7 +1964,7 @@
1.1026 /// \brief Return the matching arc (or edge) incident to the given node.
1.1027 ///
1.1028 /// This function returns the matching arc (or edge) incident to the
1.1029 - /// given node in the found matching or \c INVALID if the node is
1.1030 + /// given node in the found matching or \c INVALID if the node is
1.1031 /// not covered by the matching.
1.1032 ///
1.1033 /// \pre Either run() or start() must be called before using this function.
1.1034 @@ -1936,7 +1982,7 @@
1.1035
1.1036 /// \brief Return the mate of the given node.
1.1037 ///
1.1038 - /// This function returns the mate of the given node in the found
1.1039 + /// This function returns the mate of the given node in the found
1.1040 /// matching or \c INVALID if the node is not covered by the matching.
1.1041 ///
1.1042 /// \pre Either run() or start() must be called before using this function.
1.1043 @@ -1956,8 +2002,8 @@
1.1044
1.1045 /// \brief Return the value of the dual solution.
1.1046 ///
1.1047 - /// This function returns the value of the dual solution.
1.1048 - /// It should be equal to the primal value scaled by \ref dualScale
1.1049 + /// This function returns the value of the dual solution.
1.1050 + /// It should be equal to the primal value scaled by \ref dualScale
1.1051 /// "dual scale".
1.1052 ///
1.1053 /// \pre Either run() or start() must be called before using this function.
1.1054 @@ -2012,9 +2058,9 @@
1.1055
1.1056 /// \brief Iterator for obtaining the nodes of a blossom.
1.1057 ///
1.1058 - /// This class provides an iterator for obtaining the nodes of the
1.1059 + /// This class provides an iterator for obtaining the nodes of the
1.1060 /// given blossom. It lists a subset of the nodes.
1.1061 - /// Before using this iterator, you must allocate a
1.1062 + /// Before using this iterator, you must allocate a
1.1063 /// MaxWeightedMatching class and execute it.
1.1064 class BlossomIt {
1.1065 public:
1.1066 @@ -2023,8 +2069,8 @@
1.1067 ///
1.1068 /// Constructor to get the nodes of the given variable.
1.1069 ///
1.1070 - /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
1.1071 - /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
1.1072 + /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
1.1073 + /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
1.1074 /// called before initializing this iterator.
1.1075 BlossomIt(const MaxWeightedMatching& algorithm, int variable)
1.1076 : _algorithm(&algorithm)
1.1077 @@ -2077,8 +2123,8 @@
1.1078 /// is based on extensive use of priority queues and provides
1.1079 /// \f$O(nm\log n)\f$ time complexity.
1.1080 ///
1.1081 - /// The maximum weighted perfect matching problem is to find a subset of
1.1082 - /// the edges in an undirected graph with maximum overall weight for which
1.1083 + /// The maximum weighted perfect matching problem is to find a subset of
1.1084 + /// the edges in an undirected graph with maximum overall weight for which
1.1085 /// each node has exactly one incident edge.
1.1086 /// It can be formulated with the following linear program.
1.1087 /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1.1088 @@ -2101,16 +2147,16 @@
1.1089 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
1.1090 \frac{\vert B \vert - 1}{2}z_B\f] */
1.1091 ///
1.1092 - /// The algorithm can be executed with the run() function.
1.1093 + /// The algorithm can be executed with the run() function.
1.1094 /// After it the matching (the primal solution) and the dual solution
1.1095 - /// can be obtained using the query functions and the
1.1096 - /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
1.1097 - /// which is able to iterate on the nodes of a blossom.
1.1098 + /// can be obtained using the query functions and the
1.1099 + /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
1.1100 + /// which is able to iterate on the nodes of a blossom.
1.1101 /// If the value type is integer, then the dual solution is multiplied
1.1102 /// by \ref MaxWeightedMatching::dualScale "4".
1.1103 ///
1.1104 /// \tparam GR The undirected graph type the algorithm runs on.
1.1105 - /// \tparam WM The type edge weight map. The default type is
1.1106 + /// \tparam WM The type edge weight map. The default type is
1.1107 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
1.1108 #ifdef DOXYGEN
1.1109 template <typename GR, typename WM>
1.1110 @@ -2221,6 +2267,11 @@
1.1111 BinHeap<Value, IntIntMap> *_delta4;
1.1112
1.1113 Value _delta_sum;
1.1114 + int _unmatched;
1.1115 +
1.1116 + typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap>
1.1117 + FractionalMatching;
1.1118 + FractionalMatching *_fractional;
1.1119
1.1120 void createStructures() {
1.1121 _node_num = countNodes(_graph);
1.1122 @@ -2282,9 +2333,6 @@
1.1123 }
1.1124
1.1125 void destroyStructures() {
1.1126 - _node_num = countNodes(_graph);
1.1127 - _blossom_num = _node_num * 3 / 2;
1.1128 -
1.1129 if (_matching) {
1.1130 delete _matching;
1.1131 }
1.1132 @@ -2957,10 +3005,16 @@
1.1133 _delta3_index(0), _delta3(0),
1.1134 _delta4_index(0), _delta4(0),
1.1135
1.1136 - _delta_sum() {}
1.1137 + _delta_sum(), _unmatched(0),
1.1138 +
1.1139 + _fractional(0)
1.1140 + {}
1.1141
1.1142 ~MaxWeightedPerfectMatching() {
1.1143 destroyStructures();
1.1144 + if (_fractional) {
1.1145 + delete _fractional;
1.1146 + }
1.1147 }
1.1148
1.1149 /// \name Execution Control
1.1150 @@ -2989,6 +3043,8 @@
1.1151 (*_delta4_index)[i] = _delta4->PRE_HEAP;
1.1152 }
1.1153
1.1154 + _unmatched = _node_num;
1.1155 +
1.1156 _delta2->clear();
1.1157 _delta3->clear();
1.1158 _delta4->clear();
1.1159 @@ -3030,18 +3086,163 @@
1.1160 }
1.1161 }
1.1162
1.1163 + /// \brief Initialize the algorithm with fractional matching
1.1164 + ///
1.1165 + /// This function initializes the algorithm with a fractional
1.1166 + /// matching. This initialization is also called jumpstart heuristic.
1.1167 + void fractionalInit() {
1.1168 + createStructures();
1.1169 +
1.1170 + _blossom_node_list.clear();
1.1171 + _blossom_potential.clear();
1.1172 +
1.1173 + if (_fractional == 0) {
1.1174 + _fractional = new FractionalMatching(_graph, _weight, false);
1.1175 + }
1.1176 + if (!_fractional->run()) {
1.1177 + _unmatched = -1;
1.1178 + return;
1.1179 + }
1.1180 +
1.1181 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.1182 + (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
1.1183 + }
1.1184 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.1185 + (*_delta3_index)[e] = _delta3->PRE_HEAP;
1.1186 + }
1.1187 + for (int i = 0; i < _blossom_num; ++i) {
1.1188 + (*_delta2_index)[i] = _delta2->PRE_HEAP;
1.1189 + (*_delta4_index)[i] = _delta4->PRE_HEAP;
1.1190 + }
1.1191 +
1.1192 + _unmatched = 0;
1.1193 +
1.1194 + _delta2->clear();
1.1195 + _delta3->clear();
1.1196 + _delta4->clear();
1.1197 + _blossom_set->clear();
1.1198 + _tree_set->clear();
1.1199 +
1.1200 + int index = 0;
1.1201 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1202 + Value pot = _fractional->nodeValue(n);
1.1203 + (*_node_index)[n] = index;
1.1204 + (*_node_data)[index].pot = pot;
1.1205 + (*_node_data)[index].heap_index.clear();
1.1206 + (*_node_data)[index].heap.clear();
1.1207 + int blossom =
1.1208 + _blossom_set->insert(n, std::numeric_limits<Value>::max());
1.1209 +
1.1210 + (*_blossom_data)[blossom].status = MATCHED;
1.1211 + (*_blossom_data)[blossom].pred = INVALID;
1.1212 + (*_blossom_data)[blossom].next = _fractional->matching(n);
1.1213 + (*_blossom_data)[blossom].pot = 0;
1.1214 + (*_blossom_data)[blossom].offset = 0;
1.1215 + ++index;
1.1216 + }
1.1217 +
1.1218 + typename Graph::template NodeMap<bool> processed(_graph, false);
1.1219 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1220 + if (processed[n]) continue;
1.1221 + processed[n] = true;
1.1222 + if (_fractional->matching(n) == INVALID) continue;
1.1223 + int num = 1;
1.1224 + Node v = _graph.target(_fractional->matching(n));
1.1225 + while (n != v) {
1.1226 + processed[v] = true;
1.1227 + v = _graph.target(_fractional->matching(v));
1.1228 + ++num;
1.1229 + }
1.1230 +
1.1231 + if (num % 2 == 1) {
1.1232 + std::vector<int> subblossoms(num);
1.1233 +
1.1234 + subblossoms[--num] = _blossom_set->find(n);
1.1235 + v = _graph.target(_fractional->matching(n));
1.1236 + while (n != v) {
1.1237 + subblossoms[--num] = _blossom_set->find(v);
1.1238 + v = _graph.target(_fractional->matching(v));
1.1239 + }
1.1240 +
1.1241 + int surface =
1.1242 + _blossom_set->join(subblossoms.begin(), subblossoms.end());
1.1243 + (*_blossom_data)[surface].status = EVEN;
1.1244 + (*_blossom_data)[surface].pred = INVALID;
1.1245 + (*_blossom_data)[surface].next = INVALID;
1.1246 + (*_blossom_data)[surface].pot = 0;
1.1247 + (*_blossom_data)[surface].offset = 0;
1.1248 +
1.1249 + _tree_set->insert(surface);
1.1250 + ++_unmatched;
1.1251 + }
1.1252 + }
1.1253 +
1.1254 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.1255 + int si = (*_node_index)[_graph.u(e)];
1.1256 + int sb = _blossom_set->find(_graph.u(e));
1.1257 + int ti = (*_node_index)[_graph.v(e)];
1.1258 + int tb = _blossom_set->find(_graph.v(e));
1.1259 + if ((*_blossom_data)[sb].status == EVEN &&
1.1260 + (*_blossom_data)[tb].status == EVEN && sb != tb) {
1.1261 + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1.1262 + dualScale * _weight[e]) / 2);
1.1263 + }
1.1264 + }
1.1265 +
1.1266 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1267 + int nb = _blossom_set->find(n);
1.1268 + if ((*_blossom_data)[nb].status != MATCHED) continue;
1.1269 + int ni = (*_node_index)[n];
1.1270 +
1.1271 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.1272 + Node v = _graph.target(e);
1.1273 + int vb = _blossom_set->find(v);
1.1274 + int vi = (*_node_index)[v];
1.1275 +
1.1276 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.1277 + dualScale * _weight[e];
1.1278 +
1.1279 + if ((*_blossom_data)[vb].status == EVEN) {
1.1280 +
1.1281 + int vt = _tree_set->find(vb);
1.1282 +
1.1283 + typename std::map<int, Arc>::iterator it =
1.1284 + (*_node_data)[ni].heap_index.find(vt);
1.1285 +
1.1286 + if (it != (*_node_data)[ni].heap_index.end()) {
1.1287 + if ((*_node_data)[ni].heap[it->second] > rw) {
1.1288 + (*_node_data)[ni].heap.replace(it->second, e);
1.1289 + (*_node_data)[ni].heap.decrease(e, rw);
1.1290 + it->second = e;
1.1291 + }
1.1292 + } else {
1.1293 + (*_node_data)[ni].heap.push(e, rw);
1.1294 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
1.1295 + }
1.1296 + }
1.1297 + }
1.1298 +
1.1299 + if (!(*_node_data)[ni].heap.empty()) {
1.1300 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1.1301 + _delta2->push(nb, _blossom_set->classPrio(nb));
1.1302 + }
1.1303 + }
1.1304 + }
1.1305 +
1.1306 /// \brief Start the algorithm
1.1307 ///
1.1308 /// This function starts the algorithm.
1.1309 ///
1.1310 - /// \pre \ref init() must be called before using this function.
1.1311 + /// \pre \ref init() or \ref fractionalInit() must be called before
1.1312 + /// using this function.
1.1313 bool start() {
1.1314 enum OpType {
1.1315 D2, D3, D4
1.1316 };
1.1317
1.1318 - int unmatched = _node_num;
1.1319 - while (unmatched > 0) {
1.1320 + if (_unmatched == -1) return false;
1.1321 +
1.1322 + while (_unmatched > 0) {
1.1323 Value d2 = !_delta2->empty() ?
1.1324 _delta2->prio() : std::numeric_limits<Value>::max();
1.1325
1.1326 @@ -3051,8 +3252,8 @@
1.1327 Value d4 = !_delta4->empty() ?
1.1328 _delta4->prio() : std::numeric_limits<Value>::max();
1.1329
1.1330 - _delta_sum = d2; OpType ot = D2;
1.1331 - if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
1.1332 + _delta_sum = d3; OpType ot = D3;
1.1333 + if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1.1334 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1.1335
1.1336 if (_delta_sum == std::numeric_limits<Value>::max()) {
1.1337 @@ -3085,7 +3286,7 @@
1.1338 shrinkOnEdge(e, left_tree);
1.1339 } else {
1.1340 augmentOnEdge(e);
1.1341 - unmatched -= 2;
1.1342 + _unmatched -= 2;
1.1343 }
1.1344 }
1.1345 } break;
1.1346 @@ -3104,18 +3305,18 @@
1.1347 ///
1.1348 /// \note mwpm.run() is just a shortcut of the following code.
1.1349 /// \code
1.1350 - /// mwpm.init();
1.1351 + /// mwpm.fractionalInit();
1.1352 /// mwpm.start();
1.1353 /// \endcode
1.1354 bool run() {
1.1355 - init();
1.1356 + fractionalInit();
1.1357 return start();
1.1358 }
1.1359
1.1360 /// @}
1.1361
1.1362 /// \name Primal Solution
1.1363 - /// Functions to get the primal solution, i.e. the maximum weighted
1.1364 + /// Functions to get the primal solution, i.e. the maximum weighted
1.1365 /// perfect matching.\n
1.1366 /// Either \ref run() or \ref start() function should be called before
1.1367 /// using them.
1.1368 @@ -3134,12 +3335,12 @@
1.1369 sum += _weight[(*_matching)[n]];
1.1370 }
1.1371 }
1.1372 - return sum /= 2;
1.1373 + return sum / 2;
1.1374 }
1.1375
1.1376 /// \brief Return \c true if the given edge is in the matching.
1.1377 ///
1.1378 - /// This function returns \c true if the given edge is in the found
1.1379 + /// This function returns \c true if the given edge is in the found
1.1380 /// matching.
1.1381 ///
1.1382 /// \pre Either run() or start() must be called before using this function.
1.1383 @@ -3150,7 +3351,7 @@
1.1384 /// \brief Return the matching arc (or edge) incident to the given node.
1.1385 ///
1.1386 /// This function returns the matching arc (or edge) incident to the
1.1387 - /// given node in the found matching or \c INVALID if the node is
1.1388 + /// given node in the found matching or \c INVALID if the node is
1.1389 /// not covered by the matching.
1.1390 ///
1.1391 /// \pre Either run() or start() must be called before using this function.
1.1392 @@ -3168,7 +3369,7 @@
1.1393
1.1394 /// \brief Return the mate of the given node.
1.1395 ///
1.1396 - /// This function returns the mate of the given node in the found
1.1397 + /// This function returns the mate of the given node in the found
1.1398 /// matching or \c INVALID if the node is not covered by the matching.
1.1399 ///
1.1400 /// \pre Either run() or start() must be called before using this function.
1.1401 @@ -3187,8 +3388,8 @@
1.1402
1.1403 /// \brief Return the value of the dual solution.
1.1404 ///
1.1405 - /// This function returns the value of the dual solution.
1.1406 - /// It should be equal to the primal value scaled by \ref dualScale
1.1407 + /// This function returns the value of the dual solution.
1.1408 + /// It should be equal to the primal value scaled by \ref dualScale
1.1409 /// "dual scale".
1.1410 ///
1.1411 /// \pre Either run() or start() must be called before using this function.
1.1412 @@ -3243,9 +3444,9 @@
1.1413
1.1414 /// \brief Iterator for obtaining the nodes of a blossom.
1.1415 ///
1.1416 - /// This class provides an iterator for obtaining the nodes of the
1.1417 + /// This class provides an iterator for obtaining the nodes of the
1.1418 /// given blossom. It lists a subset of the nodes.
1.1419 - /// Before using this iterator, you must allocate a
1.1420 + /// Before using this iterator, you must allocate a
1.1421 /// MaxWeightedPerfectMatching class and execute it.
1.1422 class BlossomIt {
1.1423 public:
1.1424 @@ -3254,8 +3455,8 @@
1.1425 ///
1.1426 /// Constructor to get the nodes of the given variable.
1.1427 ///
1.1428 - /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
1.1429 - /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
1.1430 + /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
1.1431 + /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
1.1432 /// must be called before initializing this iterator.
1.1433 BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
1.1434 : _algorithm(&algorithm)
1.1435 @@ -3301,4 +3502,4 @@
1.1436
1.1437 } //END OF NAMESPACE LEMON
1.1438
1.1439 -#endif //LEMON_MAX_MATCHING_H
1.1440 +#endif //LEMON_MATCHING_H