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