1.1 --- a/lemon/matching.h Mon Aug 31 20:27:38 2009 +0200
1.2 +++ b/lemon/matching.h Sun Sep 20 21:38:24 2009 +0200
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 @@ -41,7 +41,7 @@
1.15 ///
1.16 /// This class implements Edmonds' alternating forest matching algorithm
1.17 /// for finding a maximum cardinality matching in a general undirected graph.
1.18 - /// It can be started from an arbitrary initial matching
1.19 + /// It can be started from an arbitrary initial matching
1.20 /// (the default is the empty one).
1.21 ///
1.22 /// The dual solution of the problem is a map of the nodes to
1.23 @@ -69,11 +69,11 @@
1.24
1.25 ///\brief Status constants for Gallai-Edmonds decomposition.
1.26 ///
1.27 - ///These constants are used for indicating the Gallai-Edmonds
1.28 + ///These constants are used for indicating the Gallai-Edmonds
1.29 ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
1.30 ///induce a subgraph with factor-critical components, the nodes with
1.31 ///status \c ODD (or \c A) form the canonical barrier, and the nodes
1.32 - ///with status \c MATCHED (or \c C) induce a subgraph having a
1.33 + ///with status \c MATCHED (or \c C) induce a subgraph having a
1.34 ///perfect matching.
1.35 enum Status {
1.36 EVEN = 1, ///< = 1. (\c D is an alias for \c EVEN.)
1.37 @@ -512,7 +512,7 @@
1.38 }
1.39 }
1.40
1.41 - /// \brief Start Edmonds' algorithm with a heuristic improvement
1.42 + /// \brief Start Edmonds' algorithm with a heuristic improvement
1.43 /// for dense graphs
1.44 ///
1.45 /// This function runs Edmonds' algorithm with a heuristic of postponing
1.46 @@ -534,8 +534,8 @@
1.47
1.48 /// \brief Run Edmonds' algorithm
1.49 ///
1.50 - /// This function runs Edmonds' algorithm. An additional heuristic of
1.51 - /// postponing shrinks is used for relatively dense graphs
1.52 + /// This function runs Edmonds' algorithm. An additional heuristic of
1.53 + /// postponing shrinks is used for relatively dense graphs
1.54 /// (for which <tt>m>=2*n</tt> holds).
1.55 void run() {
1.56 if (countEdges(_graph) < 2 * countNodes(_graph)) {
1.57 @@ -556,7 +556,7 @@
1.58
1.59 /// \brief Return the size (cardinality) of the matching.
1.60 ///
1.61 - /// This function returns the size (cardinality) of the current matching.
1.62 + /// This function returns the size (cardinality) of the current matching.
1.63 /// After run() it returns the size of the maximum matching in the graph.
1.64 int matchingSize() const {
1.65 int size = 0;
1.66 @@ -570,7 +570,7 @@
1.67
1.68 /// \brief Return \c true if the given edge is in the matching.
1.69 ///
1.70 - /// This function returns \c true if the given edge is in the current
1.71 + /// This function returns \c true if the given edge is in the current
1.72 /// matching.
1.73 bool matching(const Edge& edge) const {
1.74 return edge == (*_matching)[_graph.u(edge)];
1.75 @@ -579,7 +579,7 @@
1.76 /// \brief Return the matching arc (or edge) incident to the given node.
1.77 ///
1.78 /// This function returns the matching arc (or edge) incident to the
1.79 - /// given node in the current matching or \c INVALID if the node is
1.80 + /// given node in the current matching or \c INVALID if the node is
1.81 /// not covered by the matching.
1.82 Arc matching(const Node& n) const {
1.83 return (*_matching)[n];
1.84 @@ -595,7 +595,7 @@
1.85
1.86 /// \brief Return the mate of the given node.
1.87 ///
1.88 - /// This function returns the mate of the given node in the current
1.89 + /// This function returns the mate of the given node in the current
1.90 /// matching or \c INVALID if the node is not covered by the matching.
1.91 Node mate(const Node& n) const {
1.92 return (*_matching)[n] != INVALID ?
1.93 @@ -605,7 +605,7 @@
1.94 /// @}
1.95
1.96 /// \name Dual Solution
1.97 - /// Functions to get the dual solution, i.e. the Gallai-Edmonds
1.98 + /// Functions to get the dual solution, i.e. the Gallai-Edmonds
1.99 /// decomposition.
1.100
1.101 /// @{
1.102 @@ -648,8 +648,8 @@
1.103 /// on extensive use of priority queues and provides
1.104 /// \f$O(nm\log n)\f$ time complexity.
1.105 ///
1.106 - /// The maximum weighted matching problem is to find a subset of the
1.107 - /// edges in an undirected graph with maximum overall weight for which
1.108 + /// The maximum weighted matching problem is to find a subset of the
1.109 + /// edges in an undirected graph with maximum overall weight for which
1.110 /// each node has at most one incident edge.
1.111 /// It can be formulated with the following linear program.
1.112 /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
1.113 @@ -673,16 +673,16 @@
1.114 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
1.115 \frac{\vert B \vert - 1}{2}z_B\f] */
1.116 ///
1.117 - /// The algorithm can be executed with the run() function.
1.118 + /// The algorithm can be executed with the run() function.
1.119 /// After it the matching (the primal solution) and the dual solution
1.120 - /// can be obtained using the query functions and the
1.121 - /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
1.122 - /// which is able to iterate on the nodes of a blossom.
1.123 + /// can be obtained using the query functions and the
1.124 + /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
1.125 + /// which is able to iterate on the nodes of a blossom.
1.126 /// If the value type is integer, then the dual solution is multiplied
1.127 /// by \ref MaxWeightedMatching::dualScale "4".
1.128 ///
1.129 /// \tparam GR The undirected graph type the algorithm runs on.
1.130 - /// \tparam WM The type edge weight map. The default type is
1.131 + /// \tparam WM The type edge weight map. The default type is
1.132 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
1.133 #ifdef DOXYGEN
1.134 template <typename GR, typename WM>
1.135 @@ -745,7 +745,7 @@
1.136 typedef RangeMap<int> IntIntMap;
1.137
1.138 enum Status {
1.139 - EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
1.140 + EVEN = -1, MATCHED = 0, ODD = 1
1.141 };
1.142
1.143 typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
1.144 @@ -844,9 +844,6 @@
1.145 }
1.146
1.147 void destroyStructures() {
1.148 - _node_num = countNodes(_graph);
1.149 - _blossom_num = _node_num * 3 / 2;
1.150 -
1.151 if (_matching) {
1.152 delete _matching;
1.153 }
1.154 @@ -922,10 +919,6 @@
1.155 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1.156 _delta3->push(e, rw / 2);
1.157 }
1.158 - } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1.159 - if (_delta3->state(e) != _delta3->IN_HEAP) {
1.160 - _delta3->push(e, rw);
1.161 - }
1.162 } else {
1.163 typename std::map<int, Arc>::iterator it =
1.164 (*_node_data)[vi].heap_index.find(tree);
1.165 @@ -949,202 +942,6 @@
1.166 _delta2->push(vb, _blossom_set->classPrio(vb) -
1.167 (*_blossom_data)[vb].offset);
1.168 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1.169 - (*_blossom_data)[vb].offset){
1.170 - _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1.171 - (*_blossom_data)[vb].offset);
1.172 - }
1.173 - }
1.174 - }
1.175 - }
1.176 - }
1.177 - }
1.178 - (*_blossom_data)[blossom].offset = 0;
1.179 - }
1.180 -
1.181 - void matchedToOdd(int blossom) {
1.182 - if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.183 - _delta2->erase(blossom);
1.184 - }
1.185 - (*_blossom_data)[blossom].offset += _delta_sum;
1.186 - if (!_blossom_set->trivial(blossom)) {
1.187 - _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
1.188 - (*_blossom_data)[blossom].offset);
1.189 - }
1.190 - }
1.191 -
1.192 - void evenToMatched(int blossom, int tree) {
1.193 - if (!_blossom_set->trivial(blossom)) {
1.194 - (*_blossom_data)[blossom].pot += 2 * _delta_sum;
1.195 - }
1.196 -
1.197 - for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.198 - n != INVALID; ++n) {
1.199 - int ni = (*_node_index)[n];
1.200 - (*_node_data)[ni].pot -= _delta_sum;
1.201 -
1.202 - _delta1->erase(n);
1.203 -
1.204 - for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.205 - Node v = _graph.source(e);
1.206 - int vb = _blossom_set->find(v);
1.207 - int vi = (*_node_index)[v];
1.208 -
1.209 - Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.210 - dualScale * _weight[e];
1.211 -
1.212 - if (vb == blossom) {
1.213 - if (_delta3->state(e) == _delta3->IN_HEAP) {
1.214 - _delta3->erase(e);
1.215 - }
1.216 - } else if ((*_blossom_data)[vb].status == EVEN) {
1.217 -
1.218 - if (_delta3->state(e) == _delta3->IN_HEAP) {
1.219 - _delta3->erase(e);
1.220 - }
1.221 -
1.222 - int vt = _tree_set->find(vb);
1.223 -
1.224 - if (vt != tree) {
1.225 -
1.226 - Arc r = _graph.oppositeArc(e);
1.227 -
1.228 - typename std::map<int, Arc>::iterator it =
1.229 - (*_node_data)[ni].heap_index.find(vt);
1.230 -
1.231 - if (it != (*_node_data)[ni].heap_index.end()) {
1.232 - if ((*_node_data)[ni].heap[it->second] > rw) {
1.233 - (*_node_data)[ni].heap.replace(it->second, r);
1.234 - (*_node_data)[ni].heap.decrease(r, rw);
1.235 - it->second = r;
1.236 - }
1.237 - } else {
1.238 - (*_node_data)[ni].heap.push(r, rw);
1.239 - (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1.240 - }
1.241 -
1.242 - if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1.243 - _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1.244 -
1.245 - if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1.246 - _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.247 - (*_blossom_data)[blossom].offset);
1.248 - } else if ((*_delta2)[blossom] >
1.249 - _blossom_set->classPrio(blossom) -
1.250 - (*_blossom_data)[blossom].offset){
1.251 - _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1.252 - (*_blossom_data)[blossom].offset);
1.253 - }
1.254 - }
1.255 - }
1.256 -
1.257 - } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1.258 - if (_delta3->state(e) == _delta3->IN_HEAP) {
1.259 - _delta3->erase(e);
1.260 - }
1.261 - } else {
1.262 -
1.263 - typename std::map<int, Arc>::iterator it =
1.264 - (*_node_data)[vi].heap_index.find(tree);
1.265 -
1.266 - if (it != (*_node_data)[vi].heap_index.end()) {
1.267 - (*_node_data)[vi].heap.erase(it->second);
1.268 - (*_node_data)[vi].heap_index.erase(it);
1.269 - if ((*_node_data)[vi].heap.empty()) {
1.270 - _blossom_set->increase(v, std::numeric_limits<Value>::max());
1.271 - } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1.272 - _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1.273 - }
1.274 -
1.275 - if ((*_blossom_data)[vb].status == MATCHED) {
1.276 - if (_blossom_set->classPrio(vb) ==
1.277 - std::numeric_limits<Value>::max()) {
1.278 - _delta2->erase(vb);
1.279 - } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1.280 - (*_blossom_data)[vb].offset) {
1.281 - _delta2->increase(vb, _blossom_set->classPrio(vb) -
1.282 - (*_blossom_data)[vb].offset);
1.283 - }
1.284 - }
1.285 - }
1.286 - }
1.287 - }
1.288 - }
1.289 - }
1.290 -
1.291 - void oddToMatched(int blossom) {
1.292 - (*_blossom_data)[blossom].offset -= _delta_sum;
1.293 -
1.294 - if (_blossom_set->classPrio(blossom) !=
1.295 - std::numeric_limits<Value>::max()) {
1.296 - _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.297 - (*_blossom_data)[blossom].offset);
1.298 - }
1.299 -
1.300 - if (!_blossom_set->trivial(blossom)) {
1.301 - _delta4->erase(blossom);
1.302 - }
1.303 - }
1.304 -
1.305 - void oddToEven(int blossom, int tree) {
1.306 - if (!_blossom_set->trivial(blossom)) {
1.307 - _delta4->erase(blossom);
1.308 - (*_blossom_data)[blossom].pot -=
1.309 - 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1.310 - }
1.311 -
1.312 - for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.313 - n != INVALID; ++n) {
1.314 - int ni = (*_node_index)[n];
1.315 -
1.316 - _blossom_set->increase(n, std::numeric_limits<Value>::max());
1.317 -
1.318 - (*_node_data)[ni].heap.clear();
1.319 - (*_node_data)[ni].heap_index.clear();
1.320 - (*_node_data)[ni].pot +=
1.321 - 2 * _delta_sum - (*_blossom_data)[blossom].offset;
1.322 -
1.323 - _delta1->push(n, (*_node_data)[ni].pot);
1.324 -
1.325 - for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.326 - Node v = _graph.source(e);
1.327 - int vb = _blossom_set->find(v);
1.328 - int vi = (*_node_index)[v];
1.329 -
1.330 - Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.331 - dualScale * _weight[e];
1.332 -
1.333 - if ((*_blossom_data)[vb].status == EVEN) {
1.334 - if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1.335 - _delta3->push(e, rw / 2);
1.336 - }
1.337 - } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1.338 - if (_delta3->state(e) != _delta3->IN_HEAP) {
1.339 - _delta3->push(e, rw);
1.340 - }
1.341 - } else {
1.342 -
1.343 - typename std::map<int, Arc>::iterator it =
1.344 - (*_node_data)[vi].heap_index.find(tree);
1.345 -
1.346 - if (it != (*_node_data)[vi].heap_index.end()) {
1.347 - if ((*_node_data)[vi].heap[it->second] > rw) {
1.348 - (*_node_data)[vi].heap.replace(it->second, e);
1.349 - (*_node_data)[vi].heap.decrease(e, rw);
1.350 - it->second = e;
1.351 - }
1.352 - } else {
1.353 - (*_node_data)[vi].heap.push(e, rw);
1.354 - (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1.355 - }
1.356 -
1.357 - if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1.358 - _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1.359 -
1.360 - if ((*_blossom_data)[vb].status == MATCHED) {
1.361 - if (_delta2->state(vb) != _delta2->IN_HEAP) {
1.362 - _delta2->push(vb, _blossom_set->classPrio(vb) -
1.363 - (*_blossom_data)[vb].offset);
1.364 - } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1.365 (*_blossom_data)[vb].offset) {
1.366 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1.367 (*_blossom_data)[vb].offset);
1.368 @@ -1157,43 +954,145 @@
1.369 (*_blossom_data)[blossom].offset = 0;
1.370 }
1.371
1.372 -
1.373 - void matchedToUnmatched(int blossom) {
1.374 + void matchedToOdd(int blossom) {
1.375 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.376 _delta2->erase(blossom);
1.377 }
1.378 + (*_blossom_data)[blossom].offset += _delta_sum;
1.379 + if (!_blossom_set->trivial(blossom)) {
1.380 + _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
1.381 + (*_blossom_data)[blossom].offset);
1.382 + }
1.383 + }
1.384 +
1.385 + void evenToMatched(int blossom, int tree) {
1.386 + if (!_blossom_set->trivial(blossom)) {
1.387 + (*_blossom_data)[blossom].pot += 2 * _delta_sum;
1.388 + }
1.389
1.390 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.391 n != INVALID; ++n) {
1.392 int ni = (*_node_index)[n];
1.393 -
1.394 - _blossom_set->increase(n, std::numeric_limits<Value>::max());
1.395 -
1.396 - (*_node_data)[ni].heap.clear();
1.397 - (*_node_data)[ni].heap_index.clear();
1.398 -
1.399 - for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.400 - Node v = _graph.target(e);
1.401 + (*_node_data)[ni].pot -= _delta_sum;
1.402 +
1.403 + _delta1->erase(n);
1.404 +
1.405 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.406 + Node v = _graph.source(e);
1.407 int vb = _blossom_set->find(v);
1.408 int vi = (*_node_index)[v];
1.409
1.410 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.411 dualScale * _weight[e];
1.412
1.413 - if ((*_blossom_data)[vb].status == EVEN) {
1.414 - if (_delta3->state(e) != _delta3->IN_HEAP) {
1.415 - _delta3->push(e, rw);
1.416 + if (vb == blossom) {
1.417 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.418 + _delta3->erase(e);
1.419 + }
1.420 + } else if ((*_blossom_data)[vb].status == EVEN) {
1.421 +
1.422 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.423 + _delta3->erase(e);
1.424 + }
1.425 +
1.426 + int vt = _tree_set->find(vb);
1.427 +
1.428 + if (vt != tree) {
1.429 +
1.430 + Arc r = _graph.oppositeArc(e);
1.431 +
1.432 + typename std::map<int, Arc>::iterator it =
1.433 + (*_node_data)[ni].heap_index.find(vt);
1.434 +
1.435 + if (it != (*_node_data)[ni].heap_index.end()) {
1.436 + if ((*_node_data)[ni].heap[it->second] > rw) {
1.437 + (*_node_data)[ni].heap.replace(it->second, r);
1.438 + (*_node_data)[ni].heap.decrease(r, rw);
1.439 + it->second = r;
1.440 + }
1.441 + } else {
1.442 + (*_node_data)[ni].heap.push(r, rw);
1.443 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1.444 + }
1.445 +
1.446 + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1.447 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1.448 +
1.449 + if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1.450 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.451 + (*_blossom_data)[blossom].offset);
1.452 + } else if ((*_delta2)[blossom] >
1.453 + _blossom_set->classPrio(blossom) -
1.454 + (*_blossom_data)[blossom].offset){
1.455 + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1.456 + (*_blossom_data)[blossom].offset);
1.457 + }
1.458 + }
1.459 + }
1.460 + } else {
1.461 +
1.462 + typename std::map<int, Arc>::iterator it =
1.463 + (*_node_data)[vi].heap_index.find(tree);
1.464 +
1.465 + if (it != (*_node_data)[vi].heap_index.end()) {
1.466 + (*_node_data)[vi].heap.erase(it->second);
1.467 + (*_node_data)[vi].heap_index.erase(it);
1.468 + if ((*_node_data)[vi].heap.empty()) {
1.469 + _blossom_set->increase(v, std::numeric_limits<Value>::max());
1.470 + } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1.471 + _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1.472 + }
1.473 +
1.474 + if ((*_blossom_data)[vb].status == MATCHED) {
1.475 + if (_blossom_set->classPrio(vb) ==
1.476 + std::numeric_limits<Value>::max()) {
1.477 + _delta2->erase(vb);
1.478 + } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1.479 + (*_blossom_data)[vb].offset) {
1.480 + _delta2->increase(vb, _blossom_set->classPrio(vb) -
1.481 + (*_blossom_data)[vb].offset);
1.482 + }
1.483 + }
1.484 }
1.485 }
1.486 }
1.487 }
1.488 }
1.489
1.490 - void unmatchedToMatched(int blossom) {
1.491 + void oddToMatched(int blossom) {
1.492 + (*_blossom_data)[blossom].offset -= _delta_sum;
1.493 +
1.494 + if (_blossom_set->classPrio(blossom) !=
1.495 + std::numeric_limits<Value>::max()) {
1.496 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.497 + (*_blossom_data)[blossom].offset);
1.498 + }
1.499 +
1.500 + if (!_blossom_set->trivial(blossom)) {
1.501 + _delta4->erase(blossom);
1.502 + }
1.503 + }
1.504 +
1.505 + void oddToEven(int blossom, int tree) {
1.506 + if (!_blossom_set->trivial(blossom)) {
1.507 + _delta4->erase(blossom);
1.508 + (*_blossom_data)[blossom].pot -=
1.509 + 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1.510 + }
1.511 +
1.512 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.513 n != INVALID; ++n) {
1.514 int ni = (*_node_index)[n];
1.515
1.516 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
1.517 +
1.518 + (*_node_data)[ni].heap.clear();
1.519 + (*_node_data)[ni].heap_index.clear();
1.520 + (*_node_data)[ni].pot +=
1.521 + 2 * _delta_sum - (*_blossom_data)[blossom].offset;
1.522 +
1.523 + _delta1->push(n, (*_node_data)[ni].pot);
1.524 +
1.525 for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.526 Node v = _graph.source(e);
1.527 int vb = _blossom_set->find(v);
1.528 @@ -1202,54 +1101,44 @@
1.529 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.530 dualScale * _weight[e];
1.531
1.532 - if (vb == blossom) {
1.533 - if (_delta3->state(e) == _delta3->IN_HEAP) {
1.534 - _delta3->erase(e);
1.535 + if ((*_blossom_data)[vb].status == EVEN) {
1.536 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1.537 + _delta3->push(e, rw / 2);
1.538 }
1.539 - } else if ((*_blossom_data)[vb].status == EVEN) {
1.540 -
1.541 - if (_delta3->state(e) == _delta3->IN_HEAP) {
1.542 - _delta3->erase(e);
1.543 - }
1.544 -
1.545 - int vt = _tree_set->find(vb);
1.546 -
1.547 - Arc r = _graph.oppositeArc(e);
1.548 + } else {
1.549
1.550 typename std::map<int, Arc>::iterator it =
1.551 - (*_node_data)[ni].heap_index.find(vt);
1.552 -
1.553 - if (it != (*_node_data)[ni].heap_index.end()) {
1.554 - if ((*_node_data)[ni].heap[it->second] > rw) {
1.555 - (*_node_data)[ni].heap.replace(it->second, r);
1.556 - (*_node_data)[ni].heap.decrease(r, rw);
1.557 - it->second = r;
1.558 + (*_node_data)[vi].heap_index.find(tree);
1.559 +
1.560 + if (it != (*_node_data)[vi].heap_index.end()) {
1.561 + if ((*_node_data)[vi].heap[it->second] > rw) {
1.562 + (*_node_data)[vi].heap.replace(it->second, e);
1.563 + (*_node_data)[vi].heap.decrease(e, rw);
1.564 + it->second = e;
1.565 }
1.566 } else {
1.567 - (*_node_data)[ni].heap.push(r, rw);
1.568 - (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1.569 + (*_node_data)[vi].heap.push(e, rw);
1.570 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1.571 }
1.572
1.573 - if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1.574 - _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1.575 -
1.576 - if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1.577 - _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.578 - (*_blossom_data)[blossom].offset);
1.579 - } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
1.580 - (*_blossom_data)[blossom].offset){
1.581 - _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1.582 - (*_blossom_data)[blossom].offset);
1.583 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1.584 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1.585 +
1.586 + if ((*_blossom_data)[vb].status == MATCHED) {
1.587 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
1.588 + _delta2->push(vb, _blossom_set->classPrio(vb) -
1.589 + (*_blossom_data)[vb].offset);
1.590 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1.591 + (*_blossom_data)[vb].offset) {
1.592 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1.593 + (*_blossom_data)[vb].offset);
1.594 + }
1.595 }
1.596 }
1.597 -
1.598 - } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1.599 - if (_delta3->state(e) == _delta3->IN_HEAP) {
1.600 - _delta3->erase(e);
1.601 - }
1.602 }
1.603 }
1.604 }
1.605 + (*_blossom_data)[blossom].offset = 0;
1.606 }
1.607
1.608 void alternatePath(int even, int tree) {
1.609 @@ -1294,39 +1183,42 @@
1.610 alternatePath(blossom, tree);
1.611 destroyTree(tree);
1.612
1.613 - (*_blossom_data)[blossom].status = UNMATCHED;
1.614 (*_blossom_data)[blossom].base = node;
1.615 - matchedToUnmatched(blossom);
1.616 + (*_blossom_data)[blossom].next = INVALID;
1.617 }
1.618
1.619 -
1.620 void augmentOnEdge(const Edge& edge) {
1.621
1.622 int left = _blossom_set->find(_graph.u(edge));
1.623 int right = _blossom_set->find(_graph.v(edge));
1.624
1.625 - if ((*_blossom_data)[left].status == EVEN) {
1.626 - int left_tree = _tree_set->find(left);
1.627 - alternatePath(left, left_tree);
1.628 - destroyTree(left_tree);
1.629 - } else {
1.630 - (*_blossom_data)[left].status = MATCHED;
1.631 - unmatchedToMatched(left);
1.632 - }
1.633 -
1.634 - if ((*_blossom_data)[right].status == EVEN) {
1.635 - int right_tree = _tree_set->find(right);
1.636 - alternatePath(right, right_tree);
1.637 - destroyTree(right_tree);
1.638 - } else {
1.639 - (*_blossom_data)[right].status = MATCHED;
1.640 - unmatchedToMatched(right);
1.641 - }
1.642 + int left_tree = _tree_set->find(left);
1.643 + alternatePath(left, left_tree);
1.644 + destroyTree(left_tree);
1.645 +
1.646 + int right_tree = _tree_set->find(right);
1.647 + alternatePath(right, right_tree);
1.648 + destroyTree(right_tree);
1.649
1.650 (*_blossom_data)[left].next = _graph.direct(edge, true);
1.651 (*_blossom_data)[right].next = _graph.direct(edge, false);
1.652 }
1.653
1.654 + void augmentOnArc(const Arc& arc) {
1.655 +
1.656 + int left = _blossom_set->find(_graph.source(arc));
1.657 + int right = _blossom_set->find(_graph.target(arc));
1.658 +
1.659 + (*_blossom_data)[left].status = MATCHED;
1.660 +
1.661 + int right_tree = _tree_set->find(right);
1.662 + alternatePath(right, right_tree);
1.663 + destroyTree(right_tree);
1.664 +
1.665 + (*_blossom_data)[left].next = arc;
1.666 + (*_blossom_data)[right].next = _graph.oppositeArc(arc);
1.667 + }
1.668 +
1.669 void extendOnArc(const Arc& arc) {
1.670 int base = _blossom_set->find(_graph.target(arc));
1.671 int tree = _tree_set->find(base);
1.672 @@ -1529,7 +1421,7 @@
1.673 _tree_set->insert(sb, tree);
1.674 (*_blossom_data)[sb].pred = pred;
1.675 (*_blossom_data)[sb].next =
1.676 - _graph.oppositeArc((*_blossom_data)[tb].next);
1.677 + _graph.oppositeArc((*_blossom_data)[tb].next);
1.678
1.679 pred = (*_blossom_data)[ub].next;
1.680
1.681 @@ -1629,7 +1521,7 @@
1.682 }
1.683
1.684 for (int i = 0; i < int(blossoms.size()); ++i) {
1.685 - if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
1.686 + if ((*_blossom_data)[blossoms[i]].next != INVALID) {
1.687
1.688 Value offset = (*_blossom_data)[blossoms[i]].offset;
1.689 (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1.690 @@ -1757,12 +1649,11 @@
1.691 Value d4 = !_delta4->empty() ?
1.692 _delta4->prio() : std::numeric_limits<Value>::max();
1.693
1.694 - _delta_sum = d1; OpType ot = D1;
1.695 + _delta_sum = d3; OpType ot = D3;
1.696 + if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
1.697 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1.698 - if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
1.699 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1.700
1.701 -
1.702 switch (ot) {
1.703 case D1:
1.704 {
1.705 @@ -1775,8 +1666,13 @@
1.706 {
1.707 int blossom = _delta2->top();
1.708 Node n = _blossom_set->classTop(blossom);
1.709 - Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
1.710 - extendOnArc(e);
1.711 + Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
1.712 + if ((*_blossom_data)[blossom].next == INVALID) {
1.713 + augmentOnArc(a);
1.714 + --unmatched;
1.715 + } else {
1.716 + extendOnArc(a);
1.717 + }
1.718 }
1.719 break;
1.720 case D3:
1.721 @@ -1789,20 +1685,8 @@
1.722 if (left_blossom == right_blossom) {
1.723 _delta3->pop();
1.724 } else {
1.725 - int left_tree;
1.726 - if ((*_blossom_data)[left_blossom].status == EVEN) {
1.727 - left_tree = _tree_set->find(left_blossom);
1.728 - } else {
1.729 - left_tree = -1;
1.730 - ++unmatched;
1.731 - }
1.732 - int right_tree;
1.733 - if ((*_blossom_data)[right_blossom].status == EVEN) {
1.734 - right_tree = _tree_set->find(right_blossom);
1.735 - } else {
1.736 - right_tree = -1;
1.737 - ++unmatched;
1.738 - }
1.739 + int left_tree = _tree_set->find(left_blossom);
1.740 + int right_tree = _tree_set->find(right_blossom);
1.741
1.742 if (left_tree == right_tree) {
1.743 shrinkOnEdge(e, left_tree);
1.744 @@ -1837,7 +1721,7 @@
1.745 /// @}
1.746
1.747 /// \name Primal Solution
1.748 - /// Functions to get the primal solution, i.e. the maximum weighted
1.749 + /// Functions to get the primal solution, i.e. the maximum weighted
1.750 /// matching.\n
1.751 /// Either \ref run() or \ref start() function should be called before
1.752 /// using them.
1.753 @@ -1856,7 +1740,7 @@
1.754 sum += _weight[(*_matching)[n]];
1.755 }
1.756 }
1.757 - return sum /= 2;
1.758 + return sum / 2;
1.759 }
1.760
1.761 /// \brief Return the size (cardinality) of the matching.
1.762 @@ -1876,7 +1760,7 @@
1.763
1.764 /// \brief Return \c true if the given edge is in the matching.
1.765 ///
1.766 - /// This function returns \c true if the given edge is in the found
1.767 + /// This function returns \c true if the given edge is in the found
1.768 /// matching.
1.769 ///
1.770 /// \pre Either run() or start() must be called before using this function.
1.771 @@ -1887,7 +1771,7 @@
1.772 /// \brief Return the matching arc (or edge) incident to the given node.
1.773 ///
1.774 /// This function returns the matching arc (or edge) incident to the
1.775 - /// given node in the found matching or \c INVALID if the node is
1.776 + /// given node in the found matching or \c INVALID if the node is
1.777 /// not covered by the matching.
1.778 ///
1.779 /// \pre Either run() or start() must be called before using this function.
1.780 @@ -1905,7 +1789,7 @@
1.781
1.782 /// \brief Return the mate of the given node.
1.783 ///
1.784 - /// This function returns the mate of the given node in the found
1.785 + /// This function returns the mate of the given node in the found
1.786 /// matching or \c INVALID if the node is not covered by the matching.
1.787 ///
1.788 /// \pre Either run() or start() must be called before using this function.
1.789 @@ -1925,8 +1809,8 @@
1.790
1.791 /// \brief Return the value of the dual solution.
1.792 ///
1.793 - /// This function returns the value of the dual solution.
1.794 - /// It should be equal to the primal value scaled by \ref dualScale
1.795 + /// This function returns the value of the dual solution.
1.796 + /// It should be equal to the primal value scaled by \ref dualScale
1.797 /// "dual scale".
1.798 ///
1.799 /// \pre Either run() or start() must be called before using this function.
1.800 @@ -1981,9 +1865,9 @@
1.801
1.802 /// \brief Iterator for obtaining the nodes of a blossom.
1.803 ///
1.804 - /// This class provides an iterator for obtaining the nodes of the
1.805 + /// This class provides an iterator for obtaining the nodes of the
1.806 /// given blossom. It lists a subset of the nodes.
1.807 - /// Before using this iterator, you must allocate a
1.808 + /// Before using this iterator, you must allocate a
1.809 /// MaxWeightedMatching class and execute it.
1.810 class BlossomIt {
1.811 public:
1.812 @@ -1992,8 +1876,8 @@
1.813 ///
1.814 /// Constructor to get the nodes of the given variable.
1.815 ///
1.816 - /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
1.817 - /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
1.818 + /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
1.819 + /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
1.820 /// called before initializing this iterator.
1.821 BlossomIt(const MaxWeightedMatching& algorithm, int variable)
1.822 : _algorithm(&algorithm)
1.823 @@ -2046,8 +1930,8 @@
1.824 /// is based on extensive use of priority queues and provides
1.825 /// \f$O(nm\log n)\f$ time complexity.
1.826 ///
1.827 - /// The maximum weighted perfect matching problem is to find a subset of
1.828 - /// the edges in an undirected graph with maximum overall weight for which
1.829 + /// The maximum weighted perfect matching problem is to find a subset of
1.830 + /// the edges in an undirected graph with maximum overall weight for which
1.831 /// each node has exactly one incident edge.
1.832 /// It can be formulated with the following linear program.
1.833 /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1.834 @@ -2070,16 +1954,16 @@
1.835 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
1.836 \frac{\vert B \vert - 1}{2}z_B\f] */
1.837 ///
1.838 - /// The algorithm can be executed with the run() function.
1.839 + /// The algorithm can be executed with the run() function.
1.840 /// After it the matching (the primal solution) and the dual solution
1.841 - /// can be obtained using the query functions and the
1.842 - /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
1.843 - /// which is able to iterate on the nodes of a blossom.
1.844 + /// can be obtained using the query functions and the
1.845 + /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
1.846 + /// which is able to iterate on the nodes of a blossom.
1.847 /// If the value type is integer, then the dual solution is multiplied
1.848 /// by \ref MaxWeightedMatching::dualScale "4".
1.849 ///
1.850 /// \tparam GR The undirected graph type the algorithm runs on.
1.851 - /// \tparam WM The type edge weight map. The default type is
1.852 + /// \tparam WM The type edge weight map. The default type is
1.853 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
1.854 #ifdef DOXYGEN
1.855 template <typename GR, typename WM>
1.856 @@ -2233,9 +2117,6 @@
1.857 }
1.858
1.859 void destroyStructures() {
1.860 - _node_num = countNodes(_graph);
1.861 - _blossom_num = _node_num * 3 / 2;
1.862 -
1.863 if (_matching) {
1.864 delete _matching;
1.865 }
1.866 @@ -2991,8 +2872,8 @@
1.867 Value d4 = !_delta4->empty() ?
1.868 _delta4->prio() : std::numeric_limits<Value>::max();
1.869
1.870 - _delta_sum = d2; OpType ot = D2;
1.871 - if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
1.872 + _delta_sum = d3; OpType ot = D3;
1.873 + if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1.874 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1.875
1.876 if (_delta_sum == std::numeric_limits<Value>::max()) {
1.877 @@ -3055,7 +2936,7 @@
1.878 /// @}
1.879
1.880 /// \name Primal Solution
1.881 - /// Functions to get the primal solution, i.e. the maximum weighted
1.882 + /// Functions to get the primal solution, i.e. the maximum weighted
1.883 /// perfect matching.\n
1.884 /// Either \ref run() or \ref start() function should be called before
1.885 /// using them.
1.886 @@ -3074,12 +2955,12 @@
1.887 sum += _weight[(*_matching)[n]];
1.888 }
1.889 }
1.890 - return sum /= 2;
1.891 + return sum / 2;
1.892 }
1.893
1.894 /// \brief Return \c true if the given edge is in the matching.
1.895 ///
1.896 - /// This function returns \c true if the given edge is in the found
1.897 + /// This function returns \c true if the given edge is in the found
1.898 /// matching.
1.899 ///
1.900 /// \pre Either run() or start() must be called before using this function.
1.901 @@ -3090,7 +2971,7 @@
1.902 /// \brief Return the matching arc (or edge) incident to the given node.
1.903 ///
1.904 /// This function returns the matching arc (or edge) incident to the
1.905 - /// given node in the found matching or \c INVALID if the node is
1.906 + /// given node in the found matching or \c INVALID if the node is
1.907 /// not covered by the matching.
1.908 ///
1.909 /// \pre Either run() or start() must be called before using this function.
1.910 @@ -3108,7 +2989,7 @@
1.911
1.912 /// \brief Return the mate of the given node.
1.913 ///
1.914 - /// This function returns the mate of the given node in the found
1.915 + /// This function returns the mate of the given node in the found
1.916 /// matching or \c INVALID if the node is not covered by the matching.
1.917 ///
1.918 /// \pre Either run() or start() must be called before using this function.
1.919 @@ -3127,8 +3008,8 @@
1.920
1.921 /// \brief Return the value of the dual solution.
1.922 ///
1.923 - /// This function returns the value of the dual solution.
1.924 - /// It should be equal to the primal value scaled by \ref dualScale
1.925 + /// This function returns the value of the dual solution.
1.926 + /// It should be equal to the primal value scaled by \ref dualScale
1.927 /// "dual scale".
1.928 ///
1.929 /// \pre Either run() or start() must be called before using this function.
1.930 @@ -3183,9 +3064,9 @@
1.931
1.932 /// \brief Iterator for obtaining the nodes of a blossom.
1.933 ///
1.934 - /// This class provides an iterator for obtaining the nodes of the
1.935 + /// This class provides an iterator for obtaining the nodes of the
1.936 /// given blossom. It lists a subset of the nodes.
1.937 - /// Before using this iterator, you must allocate a
1.938 + /// Before using this iterator, you must allocate a
1.939 /// MaxWeightedPerfectMatching class and execute it.
1.940 class BlossomIt {
1.941 public:
1.942 @@ -3194,8 +3075,8 @@
1.943 ///
1.944 /// Constructor to get the nodes of the given variable.
1.945 ///
1.946 - /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
1.947 - /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
1.948 + /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
1.949 + /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
1.950 /// must be called before initializing this iterator.
1.951 BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
1.952 : _algorithm(&algorithm)
1.953 @@ -3241,4 +3122,4 @@
1.954
1.955 } //END OF NAMESPACE LEMON
1.956
1.957 -#endif //LEMON_MAX_MATCHING_H
1.958 +#endif //LEMON_MATCHING_H