1.1 --- a/lemon/Makefile.am Wed Oct 22 14:39:04 2008 +0100
1.2 +++ b/lemon/Makefile.am Wed Oct 22 14:41:18 2008 +0100
1.3 @@ -35,6 +35,7 @@
1.4 lemon/list_graph.h \
1.5 lemon/maps.h \
1.6 lemon/math.h \
1.7 + lemon/max_matching.h \
1.8 lemon/path.h \
1.9 lemon/random.h \
1.10 lemon/smart_graph.h \
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/lemon/max_matching.h Wed Oct 22 14:41:18 2008 +0100
2.3 @@ -0,0 +1,3104 @@
2.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
2.5 + *
2.6 + * This file is a part of LEMON, a generic C++ optimization library.
2.7 + *
2.8 + * Copyright (C) 2003-2008
2.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
2.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
2.11 + *
2.12 + * Permission to use, modify and distribute this software is granted
2.13 + * provided that this copyright notice appears in all copies. For
2.14 + * precise terms see the accompanying LICENSE file.
2.15 + *
2.16 + * This software is provided "AS IS" with no warranty of any kind,
2.17 + * express or implied, and with no claim as to its suitability for any
2.18 + * purpose.
2.19 + *
2.20 + */
2.21 +
2.22 +#ifndef LEMON_MAX_MATCHING_H
2.23 +#define LEMON_MAX_MATCHING_H
2.24 +
2.25 +#include <vector>
2.26 +#include <queue>
2.27 +#include <set>
2.28 +#include <limits>
2.29 +
2.30 +#include <lemon/core.h>
2.31 +#include <lemon/unionfind.h>
2.32 +#include <lemon/bin_heap.h>
2.33 +#include <lemon/maps.h>
2.34 +
2.35 +///\ingroup matching
2.36 +///\file
2.37 +///\brief Maximum matching algorithms in general graphs.
2.38 +
2.39 +namespace lemon {
2.40 +
2.41 + /// \ingroup matching
2.42 + ///
2.43 + /// \brief Edmonds' alternating forest maximum matching algorithm.
2.44 + ///
2.45 + /// This class implements Edmonds' alternating forest matching
2.46 + /// algorithm. The algorithm can be started from an arbitrary initial
2.47 + /// matching (the default is the empty one)
2.48 + ///
2.49 + /// The dual solution of the problem is a map of the nodes to
2.50 + /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c
2.51 + /// MATCHED/C showing the Gallai-Edmonds decomposition of the
2.52 + /// graph. The nodes in \c EVEN/D induce a graph with
2.53 + /// factor-critical components, the nodes in \c ODD/A form the
2.54 + /// barrier, and the nodes in \c MATCHED/C induce a graph having a
2.55 + /// perfect matching. The number of the factor-critical components
2.56 + /// minus the number of barrier nodes is a lower bound on the
2.57 + /// unmatched nodes, and the matching is optimal if and only if this bound is
2.58 + /// tight. This decomposition can be attained by calling \c
2.59 + /// decomposition() after running the algorithm.
2.60 + ///
2.61 + /// \param _Graph The graph type the algorithm runs on.
2.62 + template <typename _Graph>
2.63 + class MaxMatching {
2.64 + public:
2.65 +
2.66 + typedef _Graph Graph;
2.67 + typedef typename Graph::template NodeMap<typename Graph::Arc>
2.68 + MatchingMap;
2.69 +
2.70 + ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
2.71 + ///
2.72 + ///Indicates the Gallai-Edmonds decomposition of the graph. The
2.73 + ///nodes with Status \c EVEN/D induce a graph with factor-critical
2.74 + ///components, the nodes in \c ODD/A form the canonical barrier,
2.75 + ///and the nodes in \c MATCHED/C induce a graph having a perfect
2.76 + ///matching.
2.77 + enum Status {
2.78 + EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2
2.79 + };
2.80 +
2.81 + typedef typename Graph::template NodeMap<Status> StatusMap;
2.82 +
2.83 + private:
2.84 +
2.85 + TEMPLATE_GRAPH_TYPEDEFS(Graph);
2.86 +
2.87 + typedef UnionFindEnum<IntNodeMap> BlossomSet;
2.88 + typedef ExtendFindEnum<IntNodeMap> TreeSet;
2.89 + typedef RangeMap<Node> NodeIntMap;
2.90 + typedef MatchingMap EarMap;
2.91 + typedef std::vector<Node> NodeQueue;
2.92 +
2.93 + const Graph& _graph;
2.94 + MatchingMap* _matching;
2.95 + StatusMap* _status;
2.96 +
2.97 + EarMap* _ear;
2.98 +
2.99 + IntNodeMap* _blossom_set_index;
2.100 + BlossomSet* _blossom_set;
2.101 + NodeIntMap* _blossom_rep;
2.102 +
2.103 + IntNodeMap* _tree_set_index;
2.104 + TreeSet* _tree_set;
2.105 +
2.106 + NodeQueue _node_queue;
2.107 + int _process, _postpone, _last;
2.108 +
2.109 + int _node_num;
2.110 +
2.111 + private:
2.112 +
2.113 + void createStructures() {
2.114 + _node_num = countNodes(_graph);
2.115 + if (!_matching) {
2.116 + _matching = new MatchingMap(_graph);
2.117 + }
2.118 + if (!_status) {
2.119 + _status = new StatusMap(_graph);
2.120 + }
2.121 + if (!_ear) {
2.122 + _ear = new EarMap(_graph);
2.123 + }
2.124 + if (!_blossom_set) {
2.125 + _blossom_set_index = new IntNodeMap(_graph);
2.126 + _blossom_set = new BlossomSet(*_blossom_set_index);
2.127 + }
2.128 + if (!_blossom_rep) {
2.129 + _blossom_rep = new NodeIntMap(_node_num);
2.130 + }
2.131 + if (!_tree_set) {
2.132 + _tree_set_index = new IntNodeMap(_graph);
2.133 + _tree_set = new TreeSet(*_tree_set_index);
2.134 + }
2.135 + _node_queue.resize(_node_num);
2.136 + }
2.137 +
2.138 + void destroyStructures() {
2.139 + if (_matching) {
2.140 + delete _matching;
2.141 + }
2.142 + if (_status) {
2.143 + delete _status;
2.144 + }
2.145 + if (_ear) {
2.146 + delete _ear;
2.147 + }
2.148 + if (_blossom_set) {
2.149 + delete _blossom_set;
2.150 + delete _blossom_set_index;
2.151 + }
2.152 + if (_blossom_rep) {
2.153 + delete _blossom_rep;
2.154 + }
2.155 + if (_tree_set) {
2.156 + delete _tree_set_index;
2.157 + delete _tree_set;
2.158 + }
2.159 + }
2.160 +
2.161 + void processDense(const Node& n) {
2.162 + _process = _postpone = _last = 0;
2.163 + _node_queue[_last++] = n;
2.164 +
2.165 + while (_process != _last) {
2.166 + Node u = _node_queue[_process++];
2.167 + for (OutArcIt a(_graph, u); a != INVALID; ++a) {
2.168 + Node v = _graph.target(a);
2.169 + if ((*_status)[v] == MATCHED) {
2.170 + extendOnArc(a);
2.171 + } else if ((*_status)[v] == UNMATCHED) {
2.172 + augmentOnArc(a);
2.173 + return;
2.174 + }
2.175 + }
2.176 + }
2.177 +
2.178 + while (_postpone != _last) {
2.179 + Node u = _node_queue[_postpone++];
2.180 +
2.181 + for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
2.182 + Node v = _graph.target(a);
2.183 +
2.184 + if ((*_status)[v] == EVEN) {
2.185 + if (_blossom_set->find(u) != _blossom_set->find(v)) {
2.186 + shrinkOnEdge(a);
2.187 + }
2.188 + }
2.189 +
2.190 + while (_process != _last) {
2.191 + Node w = _node_queue[_process++];
2.192 + for (OutArcIt b(_graph, w); b != INVALID; ++b) {
2.193 + Node x = _graph.target(b);
2.194 + if ((*_status)[x] == MATCHED) {
2.195 + extendOnArc(b);
2.196 + } else if ((*_status)[x] == UNMATCHED) {
2.197 + augmentOnArc(b);
2.198 + return;
2.199 + }
2.200 + }
2.201 + }
2.202 + }
2.203 + }
2.204 + }
2.205 +
2.206 + void processSparse(const Node& n) {
2.207 + _process = _last = 0;
2.208 + _node_queue[_last++] = n;
2.209 + while (_process != _last) {
2.210 + Node u = _node_queue[_process++];
2.211 + for (OutArcIt a(_graph, u); a != INVALID; ++a) {
2.212 + Node v = _graph.target(a);
2.213 +
2.214 + if ((*_status)[v] == EVEN) {
2.215 + if (_blossom_set->find(u) != _blossom_set->find(v)) {
2.216 + shrinkOnEdge(a);
2.217 + }
2.218 + } else if ((*_status)[v] == MATCHED) {
2.219 + extendOnArc(a);
2.220 + } else if ((*_status)[v] == UNMATCHED) {
2.221 + augmentOnArc(a);
2.222 + return;
2.223 + }
2.224 + }
2.225 + }
2.226 + }
2.227 +
2.228 + void shrinkOnEdge(const Edge& e) {
2.229 + Node nca = INVALID;
2.230 +
2.231 + {
2.232 + std::set<Node> left_set, right_set;
2.233 +
2.234 + Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
2.235 + left_set.insert(left);
2.236 +
2.237 + Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
2.238 + right_set.insert(right);
2.239 +
2.240 + while (true) {
2.241 + if ((*_matching)[left] == INVALID) break;
2.242 + left = _graph.target((*_matching)[left]);
2.243 + left = (*_blossom_rep)[_blossom_set->
2.244 + find(_graph.target((*_ear)[left]))];
2.245 + if (right_set.find(left) != right_set.end()) {
2.246 + nca = left;
2.247 + break;
2.248 + }
2.249 + left_set.insert(left);
2.250 +
2.251 + if ((*_matching)[right] == INVALID) break;
2.252 + right = _graph.target((*_matching)[right]);
2.253 + right = (*_blossom_rep)[_blossom_set->
2.254 + find(_graph.target((*_ear)[right]))];
2.255 + if (left_set.find(right) != left_set.end()) {
2.256 + nca = right;
2.257 + break;
2.258 + }
2.259 + right_set.insert(right);
2.260 + }
2.261 +
2.262 + if (nca == INVALID) {
2.263 + if ((*_matching)[left] == INVALID) {
2.264 + nca = right;
2.265 + while (left_set.find(nca) == left_set.end()) {
2.266 + nca = _graph.target((*_matching)[nca]);
2.267 + nca =(*_blossom_rep)[_blossom_set->
2.268 + find(_graph.target((*_ear)[nca]))];
2.269 + }
2.270 + } else {
2.271 + nca = left;
2.272 + while (right_set.find(nca) == right_set.end()) {
2.273 + nca = _graph.target((*_matching)[nca]);
2.274 + nca = (*_blossom_rep)[_blossom_set->
2.275 + find(_graph.target((*_ear)[nca]))];
2.276 + }
2.277 + }
2.278 + }
2.279 + }
2.280 +
2.281 + {
2.282 +
2.283 + Node node = _graph.u(e);
2.284 + Arc arc = _graph.direct(e, true);
2.285 + Node base = (*_blossom_rep)[_blossom_set->find(node)];
2.286 +
2.287 + while (base != nca) {
2.288 + _ear->set(node, arc);
2.289 +
2.290 + Node n = node;
2.291 + while (n != base) {
2.292 + n = _graph.target((*_matching)[n]);
2.293 + Arc a = (*_ear)[n];
2.294 + n = _graph.target(a);
2.295 + _ear->set(n, _graph.oppositeArc(a));
2.296 + }
2.297 + node = _graph.target((*_matching)[base]);
2.298 + _tree_set->erase(base);
2.299 + _tree_set->erase(node);
2.300 + _blossom_set->insert(node, _blossom_set->find(base));
2.301 + _status->set(node, EVEN);
2.302 + _node_queue[_last++] = node;
2.303 + arc = _graph.oppositeArc((*_ear)[node]);
2.304 + node = _graph.target((*_ear)[node]);
2.305 + base = (*_blossom_rep)[_blossom_set->find(node)];
2.306 + _blossom_set->join(_graph.target(arc), base);
2.307 + }
2.308 + }
2.309 +
2.310 + _blossom_rep->set(_blossom_set->find(nca), nca);
2.311 +
2.312 + {
2.313 +
2.314 + Node node = _graph.v(e);
2.315 + Arc arc = _graph.direct(e, false);
2.316 + Node base = (*_blossom_rep)[_blossom_set->find(node)];
2.317 +
2.318 + while (base != nca) {
2.319 + _ear->set(node, arc);
2.320 +
2.321 + Node n = node;
2.322 + while (n != base) {
2.323 + n = _graph.target((*_matching)[n]);
2.324 + Arc a = (*_ear)[n];
2.325 + n = _graph.target(a);
2.326 + _ear->set(n, _graph.oppositeArc(a));
2.327 + }
2.328 + node = _graph.target((*_matching)[base]);
2.329 + _tree_set->erase(base);
2.330 + _tree_set->erase(node);
2.331 + _blossom_set->insert(node, _blossom_set->find(base));
2.332 + _status->set(node, EVEN);
2.333 + _node_queue[_last++] = node;
2.334 + arc = _graph.oppositeArc((*_ear)[node]);
2.335 + node = _graph.target((*_ear)[node]);
2.336 + base = (*_blossom_rep)[_blossom_set->find(node)];
2.337 + _blossom_set->join(_graph.target(arc), base);
2.338 + }
2.339 + }
2.340 +
2.341 + _blossom_rep->set(_blossom_set->find(nca), nca);
2.342 + }
2.343 +
2.344 +
2.345 +
2.346 + void extendOnArc(const Arc& a) {
2.347 + Node base = _graph.source(a);
2.348 + Node odd = _graph.target(a);
2.349 +
2.350 + _ear->set(odd, _graph.oppositeArc(a));
2.351 + Node even = _graph.target((*_matching)[odd]);
2.352 + _blossom_rep->set(_blossom_set->insert(even), even);
2.353 + _status->set(odd, ODD);
2.354 + _status->set(even, EVEN);
2.355 + int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
2.356 + _tree_set->insert(odd, tree);
2.357 + _tree_set->insert(even, tree);
2.358 + _node_queue[_last++] = even;
2.359 +
2.360 + }
2.361 +
2.362 + void augmentOnArc(const Arc& a) {
2.363 + Node even = _graph.source(a);
2.364 + Node odd = _graph.target(a);
2.365 +
2.366 + int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
2.367 +
2.368 + _matching->set(odd, _graph.oppositeArc(a));
2.369 + _status->set(odd, MATCHED);
2.370 +
2.371 + Arc arc = (*_matching)[even];
2.372 + _matching->set(even, a);
2.373 +
2.374 + while (arc != INVALID) {
2.375 + odd = _graph.target(arc);
2.376 + arc = (*_ear)[odd];
2.377 + even = _graph.target(arc);
2.378 + _matching->set(odd, arc);
2.379 + arc = (*_matching)[even];
2.380 + _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
2.381 + }
2.382 +
2.383 + for (typename TreeSet::ItemIt it(*_tree_set, tree);
2.384 + it != INVALID; ++it) {
2.385 + if ((*_status)[it] == ODD) {
2.386 + _status->set(it, MATCHED);
2.387 + } else {
2.388 + int blossom = _blossom_set->find(it);
2.389 + for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
2.390 + jt != INVALID; ++jt) {
2.391 + _status->set(jt, MATCHED);
2.392 + }
2.393 + _blossom_set->eraseClass(blossom);
2.394 + }
2.395 + }
2.396 + _tree_set->eraseClass(tree);
2.397 +
2.398 + }
2.399 +
2.400 + public:
2.401 +
2.402 + /// \brief Constructor
2.403 + ///
2.404 + /// Constructor.
2.405 + MaxMatching(const Graph& graph)
2.406 + : _graph(graph), _matching(0), _status(0), _ear(0),
2.407 + _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
2.408 + _tree_set_index(0), _tree_set(0) {}
2.409 +
2.410 + ~MaxMatching() {
2.411 + destroyStructures();
2.412 + }
2.413 +
2.414 + /// \name Execution control
2.415 + /// The simplest way to execute the algorithm is to use the
2.416 + /// \c run() member function.
2.417 + /// \n
2.418 +
2.419 + /// If you need better control on the execution, you must call
2.420 + /// \ref init(), \ref greedyInit() or \ref matchingInit()
2.421 + /// functions first, then you can start the algorithm with the \ref
2.422 + /// startParse() or startDense() functions.
2.423 +
2.424 + ///@{
2.425 +
2.426 + /// \brief Sets the actual matching to the empty matching.
2.427 + ///
2.428 + /// Sets the actual matching to the empty matching.
2.429 + ///
2.430 + void init() {
2.431 + createStructures();
2.432 + for(NodeIt n(_graph); n != INVALID; ++n) {
2.433 + _matching->set(n, INVALID);
2.434 + _status->set(n, UNMATCHED);
2.435 + }
2.436 + }
2.437 +
2.438 + ///\brief Finds an initial matching in a greedy way
2.439 + ///
2.440 + ///It finds an initial matching in a greedy way.
2.441 + void greedyInit() {
2.442 + createStructures();
2.443 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.444 + _matching->set(n, INVALID);
2.445 + _status->set(n, UNMATCHED);
2.446 + }
2.447 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.448 + if ((*_matching)[n] == INVALID) {
2.449 + for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
2.450 + Node v = _graph.target(a);
2.451 + if ((*_matching)[v] == INVALID && v != n) {
2.452 + _matching->set(n, a);
2.453 + _status->set(n, MATCHED);
2.454 + _matching->set(v, _graph.oppositeArc(a));
2.455 + _status->set(v, MATCHED);
2.456 + break;
2.457 + }
2.458 + }
2.459 + }
2.460 + }
2.461 + }
2.462 +
2.463 +
2.464 + /// \brief Initialize the matching from a map containing.
2.465 + ///
2.466 + /// Initialize the matching from a \c bool valued \c Edge map. This
2.467 + /// map must have the property that there are no two incident edges
2.468 + /// with true value, ie. it contains a matching.
2.469 + /// \return %True if the map contains a matching.
2.470 + template <typename MatchingMap>
2.471 + bool matchingInit(const MatchingMap& matching) {
2.472 + createStructures();
2.473 +
2.474 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.475 + _matching->set(n, INVALID);
2.476 + _status->set(n, UNMATCHED);
2.477 + }
2.478 + for(EdgeIt e(_graph); e!=INVALID; ++e) {
2.479 + if (matching[e]) {
2.480 +
2.481 + Node u = _graph.u(e);
2.482 + if ((*_matching)[u] != INVALID) return false;
2.483 + _matching->set(u, _graph.direct(e, true));
2.484 + _status->set(u, MATCHED);
2.485 +
2.486 + Node v = _graph.v(e);
2.487 + if ((*_matching)[v] != INVALID) return false;
2.488 + _matching->set(v, _graph.direct(e, false));
2.489 + _status->set(v, MATCHED);
2.490 + }
2.491 + }
2.492 + return true;
2.493 + }
2.494 +
2.495 + /// \brief Starts Edmonds' algorithm
2.496 + ///
2.497 + /// If runs the original Edmonds' algorithm.
2.498 + void startSparse() {
2.499 + for(NodeIt n(_graph); n != INVALID; ++n) {
2.500 + if ((*_status)[n] == UNMATCHED) {
2.501 + (*_blossom_rep)[_blossom_set->insert(n)] = n;
2.502 + _tree_set->insert(n);
2.503 + _status->set(n, EVEN);
2.504 + processSparse(n);
2.505 + }
2.506 + }
2.507 + }
2.508 +
2.509 + /// \brief Starts Edmonds' algorithm.
2.510 + ///
2.511 + /// It runs Edmonds' algorithm with a heuristic of postponing
2.512 + /// shrinks, therefore resulting in a faster algorithm for dense graphs.
2.513 + void startDense() {
2.514 + for(NodeIt n(_graph); n != INVALID; ++n) {
2.515 + if ((*_status)[n] == UNMATCHED) {
2.516 + (*_blossom_rep)[_blossom_set->insert(n)] = n;
2.517 + _tree_set->insert(n);
2.518 + _status->set(n, EVEN);
2.519 + processDense(n);
2.520 + }
2.521 + }
2.522 + }
2.523 +
2.524 +
2.525 + /// \brief Runs Edmonds' algorithm
2.526 + ///
2.527 + /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>)
2.528 + /// or Edmonds' algorithm with a heuristic of
2.529 + /// postponing shrinks for dense graphs.
2.530 + void run() {
2.531 + if (countEdges(_graph) < 2 * countNodes(_graph)) {
2.532 + greedyInit();
2.533 + startSparse();
2.534 + } else {
2.535 + init();
2.536 + startDense();
2.537 + }
2.538 + }
2.539 +
2.540 + /// @}
2.541 +
2.542 + /// \name Primal solution
2.543 + /// Functions to get the primal solution, ie. the matching.
2.544 +
2.545 + /// @{
2.546 +
2.547 + ///\brief Returns the size of the current matching.
2.548 + ///
2.549 + ///Returns the size of the current matching. After \ref
2.550 + ///run() it returns the size of the maximum matching in the graph.
2.551 + int matchingSize() const {
2.552 + int size = 0;
2.553 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.554 + if ((*_matching)[n] != INVALID) {
2.555 + ++size;
2.556 + }
2.557 + }
2.558 + return size / 2;
2.559 + }
2.560 +
2.561 + /// \brief Returns true when the edge is in the matching.
2.562 + ///
2.563 + /// Returns true when the edge is in the matching.
2.564 + bool matching(const Edge& edge) const {
2.565 + return edge == (*_matching)[_graph.u(edge)];
2.566 + }
2.567 +
2.568 + /// \brief Returns the matching edge incident to the given node.
2.569 + ///
2.570 + /// Returns the matching edge of a \c node in the actual matching or
2.571 + /// INVALID if the \c node is not covered by the actual matching.
2.572 + Arc matching(const Node& n) const {
2.573 + return (*_matching)[n];
2.574 + }
2.575 +
2.576 + ///\brief Returns the mate of a node in the actual matching.
2.577 + ///
2.578 + ///Returns the mate of a \c node in the actual matching or
2.579 + ///INVALID if the \c node is not covered by the actual matching.
2.580 + Node mate(const Node& n) const {
2.581 + return (*_matching)[n] != INVALID ?
2.582 + _graph.target((*_matching)[n]) : INVALID;
2.583 + }
2.584 +
2.585 + /// @}
2.586 +
2.587 + /// \name Dual solution
2.588 + /// Functions to get the dual solution, ie. the decomposition.
2.589 +
2.590 + /// @{
2.591 +
2.592 + /// \brief Returns the class of the node in the Edmonds-Gallai
2.593 + /// decomposition.
2.594 + ///
2.595 + /// Returns the class of the node in the Edmonds-Gallai
2.596 + /// decomposition.
2.597 + Status decomposition(const Node& n) const {
2.598 + return (*_status)[n];
2.599 + }
2.600 +
2.601 + /// \brief Returns true when the node is in the barrier.
2.602 + ///
2.603 + /// Returns true when the node is in the barrier.
2.604 + bool barrier(const Node& n) const {
2.605 + return (*_status)[n] == ODD;
2.606 + }
2.607 +
2.608 + /// @}
2.609 +
2.610 + };
2.611 +
2.612 + /// \ingroup matching
2.613 + ///
2.614 + /// \brief Weighted matching in general graphs
2.615 + ///
2.616 + /// This class provides an efficient implementation of Edmond's
2.617 + /// maximum weighted matching algorithm. The implementation is based
2.618 + /// on extensive use of priority queues and provides
2.619 + /// \f$O(nm\log(n))\f$ time complexity.
2.620 + ///
2.621 + /// The maximum weighted matching problem is to find undirected
2.622 + /// edges in the graph with maximum overall weight and no two of
2.623 + /// them shares their ends. The problem can be formulated with the
2.624 + /// following linear program.
2.625 + /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
2.626 + /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
2.627 + \quad \forall B\in\mathcal{O}\f] */
2.628 + /// \f[x_e \ge 0\quad \forall e\in E\f]
2.629 + /// \f[\max \sum_{e\in E}x_ew_e\f]
2.630 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
2.631 + /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
2.632 + /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
2.633 + /// subsets of the nodes.
2.634 + ///
2.635 + /// The algorithm calculates an optimal matching and a proof of the
2.636 + /// optimality. The solution of the dual problem can be used to check
2.637 + /// the result of the algorithm. The dual linear problem is the
2.638 + /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
2.639 + z_B \ge w_{uv} \quad \forall uv\in E\f] */
2.640 + /// \f[y_u \ge 0 \quad \forall u \in V\f]
2.641 + /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
2.642 + /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
2.643 + \frac{\vert B \vert - 1}{2}z_B\f] */
2.644 + ///
2.645 + /// The algorithm can be executed with \c run() or the \c init() and
2.646 + /// then the \c start() member functions. After it the matching can
2.647 + /// be asked with \c matching() or mate() functions. The dual
2.648 + /// solution can be get with \c nodeValue(), \c blossomNum() and \c
2.649 + /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
2.650 + /// "BlossomIt" nested class, which is able to iterate on the nodes
2.651 + /// of a blossom. If the value type is integral then the dual
2.652 + /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
2.653 + template <typename _Graph,
2.654 + typename _WeightMap = typename _Graph::template EdgeMap<int> >
2.655 + class MaxWeightedMatching {
2.656 + public:
2.657 +
2.658 + typedef _Graph Graph;
2.659 + typedef _WeightMap WeightMap;
2.660 + typedef typename WeightMap::Value Value;
2.661 +
2.662 + /// \brief Scaling factor for dual solution
2.663 + ///
2.664 + /// Scaling factor for dual solution, it is equal to 4 or 1
2.665 + /// according to the value type.
2.666 + static const int dualScale =
2.667 + std::numeric_limits<Value>::is_integer ? 4 : 1;
2.668 +
2.669 + typedef typename Graph::template NodeMap<typename Graph::Arc>
2.670 + MatchingMap;
2.671 +
2.672 + private:
2.673 +
2.674 + TEMPLATE_GRAPH_TYPEDEFS(Graph);
2.675 +
2.676 + typedef typename Graph::template NodeMap<Value> NodePotential;
2.677 + typedef std::vector<Node> BlossomNodeList;
2.678 +
2.679 + struct BlossomVariable {
2.680 + int begin, end;
2.681 + Value value;
2.682 +
2.683 + BlossomVariable(int _begin, int _end, Value _value)
2.684 + : begin(_begin), end(_end), value(_value) {}
2.685 +
2.686 + };
2.687 +
2.688 + typedef std::vector<BlossomVariable> BlossomPotential;
2.689 +
2.690 + const Graph& _graph;
2.691 + const WeightMap& _weight;
2.692 +
2.693 + MatchingMap* _matching;
2.694 +
2.695 + NodePotential* _node_potential;
2.696 +
2.697 + BlossomPotential _blossom_potential;
2.698 + BlossomNodeList _blossom_node_list;
2.699 +
2.700 + int _node_num;
2.701 + int _blossom_num;
2.702 +
2.703 + typedef RangeMap<int> IntIntMap;
2.704 +
2.705 + enum Status {
2.706 + EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
2.707 + };
2.708 +
2.709 + typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
2.710 + struct BlossomData {
2.711 + int tree;
2.712 + Status status;
2.713 + Arc pred, next;
2.714 + Value pot, offset;
2.715 + Node base;
2.716 + };
2.717 +
2.718 + IntNodeMap *_blossom_index;
2.719 + BlossomSet *_blossom_set;
2.720 + RangeMap<BlossomData>* _blossom_data;
2.721 +
2.722 + IntNodeMap *_node_index;
2.723 + IntArcMap *_node_heap_index;
2.724 +
2.725 + struct NodeData {
2.726 +
2.727 + NodeData(IntArcMap& node_heap_index)
2.728 + : heap(node_heap_index) {}
2.729 +
2.730 + int blossom;
2.731 + Value pot;
2.732 + BinHeap<Value, IntArcMap> heap;
2.733 + std::map<int, Arc> heap_index;
2.734 +
2.735 + int tree;
2.736 + };
2.737 +
2.738 + RangeMap<NodeData>* _node_data;
2.739 +
2.740 + typedef ExtendFindEnum<IntIntMap> TreeSet;
2.741 +
2.742 + IntIntMap *_tree_set_index;
2.743 + TreeSet *_tree_set;
2.744 +
2.745 + IntNodeMap *_delta1_index;
2.746 + BinHeap<Value, IntNodeMap> *_delta1;
2.747 +
2.748 + IntIntMap *_delta2_index;
2.749 + BinHeap<Value, IntIntMap> *_delta2;
2.750 +
2.751 + IntEdgeMap *_delta3_index;
2.752 + BinHeap<Value, IntEdgeMap> *_delta3;
2.753 +
2.754 + IntIntMap *_delta4_index;
2.755 + BinHeap<Value, IntIntMap> *_delta4;
2.756 +
2.757 + Value _delta_sum;
2.758 +
2.759 + void createStructures() {
2.760 + _node_num = countNodes(_graph);
2.761 + _blossom_num = _node_num * 3 / 2;
2.762 +
2.763 + if (!_matching) {
2.764 + _matching = new MatchingMap(_graph);
2.765 + }
2.766 + if (!_node_potential) {
2.767 + _node_potential = new NodePotential(_graph);
2.768 + }
2.769 + if (!_blossom_set) {
2.770 + _blossom_index = new IntNodeMap(_graph);
2.771 + _blossom_set = new BlossomSet(*_blossom_index);
2.772 + _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2.773 + }
2.774 +
2.775 + if (!_node_index) {
2.776 + _node_index = new IntNodeMap(_graph);
2.777 + _node_heap_index = new IntArcMap(_graph);
2.778 + _node_data = new RangeMap<NodeData>(_node_num,
2.779 + NodeData(*_node_heap_index));
2.780 + }
2.781 +
2.782 + if (!_tree_set) {
2.783 + _tree_set_index = new IntIntMap(_blossom_num);
2.784 + _tree_set = new TreeSet(*_tree_set_index);
2.785 + }
2.786 + if (!_delta1) {
2.787 + _delta1_index = new IntNodeMap(_graph);
2.788 + _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
2.789 + }
2.790 + if (!_delta2) {
2.791 + _delta2_index = new IntIntMap(_blossom_num);
2.792 + _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
2.793 + }
2.794 + if (!_delta3) {
2.795 + _delta3_index = new IntEdgeMap(_graph);
2.796 + _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
2.797 + }
2.798 + if (!_delta4) {
2.799 + _delta4_index = new IntIntMap(_blossom_num);
2.800 + _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
2.801 + }
2.802 + }
2.803 +
2.804 + void destroyStructures() {
2.805 + _node_num = countNodes(_graph);
2.806 + _blossom_num = _node_num * 3 / 2;
2.807 +
2.808 + if (_matching) {
2.809 + delete _matching;
2.810 + }
2.811 + if (_node_potential) {
2.812 + delete _node_potential;
2.813 + }
2.814 + if (_blossom_set) {
2.815 + delete _blossom_index;
2.816 + delete _blossom_set;
2.817 + delete _blossom_data;
2.818 + }
2.819 +
2.820 + if (_node_index) {
2.821 + delete _node_index;
2.822 + delete _node_heap_index;
2.823 + delete _node_data;
2.824 + }
2.825 +
2.826 + if (_tree_set) {
2.827 + delete _tree_set_index;
2.828 + delete _tree_set;
2.829 + }
2.830 + if (_delta1) {
2.831 + delete _delta1_index;
2.832 + delete _delta1;
2.833 + }
2.834 + if (_delta2) {
2.835 + delete _delta2_index;
2.836 + delete _delta2;
2.837 + }
2.838 + if (_delta3) {
2.839 + delete _delta3_index;
2.840 + delete _delta3;
2.841 + }
2.842 + if (_delta4) {
2.843 + delete _delta4_index;
2.844 + delete _delta4;
2.845 + }
2.846 + }
2.847 +
2.848 + void matchedToEven(int blossom, int tree) {
2.849 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2.850 + _delta2->erase(blossom);
2.851 + }
2.852 +
2.853 + if (!_blossom_set->trivial(blossom)) {
2.854 + (*_blossom_data)[blossom].pot -=
2.855 + 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
2.856 + }
2.857 +
2.858 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2.859 + n != INVALID; ++n) {
2.860 +
2.861 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
2.862 + int ni = (*_node_index)[n];
2.863 +
2.864 + (*_node_data)[ni].heap.clear();
2.865 + (*_node_data)[ni].heap_index.clear();
2.866 +
2.867 + (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
2.868 +
2.869 + _delta1->push(n, (*_node_data)[ni].pot);
2.870 +
2.871 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
2.872 + Node v = _graph.source(e);
2.873 + int vb = _blossom_set->find(v);
2.874 + int vi = (*_node_index)[v];
2.875 +
2.876 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2.877 + dualScale * _weight[e];
2.878 +
2.879 + if ((*_blossom_data)[vb].status == EVEN) {
2.880 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2.881 + _delta3->push(e, rw / 2);
2.882 + }
2.883 + } else if ((*_blossom_data)[vb].status == UNMATCHED) {
2.884 + if (_delta3->state(e) != _delta3->IN_HEAP) {
2.885 + _delta3->push(e, rw);
2.886 + }
2.887 + } else {
2.888 + typename std::map<int, Arc>::iterator it =
2.889 + (*_node_data)[vi].heap_index.find(tree);
2.890 +
2.891 + if (it != (*_node_data)[vi].heap_index.end()) {
2.892 + if ((*_node_data)[vi].heap[it->second] > rw) {
2.893 + (*_node_data)[vi].heap.replace(it->second, e);
2.894 + (*_node_data)[vi].heap.decrease(e, rw);
2.895 + it->second = e;
2.896 + }
2.897 + } else {
2.898 + (*_node_data)[vi].heap.push(e, rw);
2.899 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2.900 + }
2.901 +
2.902 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2.903 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2.904 +
2.905 + if ((*_blossom_data)[vb].status == MATCHED) {
2.906 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
2.907 + _delta2->push(vb, _blossom_set->classPrio(vb) -
2.908 + (*_blossom_data)[vb].offset);
2.909 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2.910 + (*_blossom_data)[vb].offset){
2.911 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2.912 + (*_blossom_data)[vb].offset);
2.913 + }
2.914 + }
2.915 + }
2.916 + }
2.917 + }
2.918 + }
2.919 + (*_blossom_data)[blossom].offset = 0;
2.920 + }
2.921 +
2.922 + void matchedToOdd(int blossom) {
2.923 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2.924 + _delta2->erase(blossom);
2.925 + }
2.926 + (*_blossom_data)[blossom].offset += _delta_sum;
2.927 + if (!_blossom_set->trivial(blossom)) {
2.928 + _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
2.929 + (*_blossom_data)[blossom].offset);
2.930 + }
2.931 + }
2.932 +
2.933 + void evenToMatched(int blossom, int tree) {
2.934 + if (!_blossom_set->trivial(blossom)) {
2.935 + (*_blossom_data)[blossom].pot += 2 * _delta_sum;
2.936 + }
2.937 +
2.938 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2.939 + n != INVALID; ++n) {
2.940 + int ni = (*_node_index)[n];
2.941 + (*_node_data)[ni].pot -= _delta_sum;
2.942 +
2.943 + _delta1->erase(n);
2.944 +
2.945 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
2.946 + Node v = _graph.source(e);
2.947 + int vb = _blossom_set->find(v);
2.948 + int vi = (*_node_index)[v];
2.949 +
2.950 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2.951 + dualScale * _weight[e];
2.952 +
2.953 + if (vb == blossom) {
2.954 + if (_delta3->state(e) == _delta3->IN_HEAP) {
2.955 + _delta3->erase(e);
2.956 + }
2.957 + } else if ((*_blossom_data)[vb].status == EVEN) {
2.958 +
2.959 + if (_delta3->state(e) == _delta3->IN_HEAP) {
2.960 + _delta3->erase(e);
2.961 + }
2.962 +
2.963 + int vt = _tree_set->find(vb);
2.964 +
2.965 + if (vt != tree) {
2.966 +
2.967 + Arc r = _graph.oppositeArc(e);
2.968 +
2.969 + typename std::map<int, Arc>::iterator it =
2.970 + (*_node_data)[ni].heap_index.find(vt);
2.971 +
2.972 + if (it != (*_node_data)[ni].heap_index.end()) {
2.973 + if ((*_node_data)[ni].heap[it->second] > rw) {
2.974 + (*_node_data)[ni].heap.replace(it->second, r);
2.975 + (*_node_data)[ni].heap.decrease(r, rw);
2.976 + it->second = r;
2.977 + }
2.978 + } else {
2.979 + (*_node_data)[ni].heap.push(r, rw);
2.980 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
2.981 + }
2.982 +
2.983 + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
2.984 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
2.985 +
2.986 + if (_delta2->state(blossom) != _delta2->IN_HEAP) {
2.987 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2.988 + (*_blossom_data)[blossom].offset);
2.989 + } else if ((*_delta2)[blossom] >
2.990 + _blossom_set->classPrio(blossom) -
2.991 + (*_blossom_data)[blossom].offset){
2.992 + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
2.993 + (*_blossom_data)[blossom].offset);
2.994 + }
2.995 + }
2.996 + }
2.997 +
2.998 + } else if ((*_blossom_data)[vb].status == UNMATCHED) {
2.999 + if (_delta3->state(e) == _delta3->IN_HEAP) {
2.1000 + _delta3->erase(e);
2.1001 + }
2.1002 + } else {
2.1003 +
2.1004 + typename std::map<int, Arc>::iterator it =
2.1005 + (*_node_data)[vi].heap_index.find(tree);
2.1006 +
2.1007 + if (it != (*_node_data)[vi].heap_index.end()) {
2.1008 + (*_node_data)[vi].heap.erase(it->second);
2.1009 + (*_node_data)[vi].heap_index.erase(it);
2.1010 + if ((*_node_data)[vi].heap.empty()) {
2.1011 + _blossom_set->increase(v, std::numeric_limits<Value>::max());
2.1012 + } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
2.1013 + _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
2.1014 + }
2.1015 +
2.1016 + if ((*_blossom_data)[vb].status == MATCHED) {
2.1017 + if (_blossom_set->classPrio(vb) ==
2.1018 + std::numeric_limits<Value>::max()) {
2.1019 + _delta2->erase(vb);
2.1020 + } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
2.1021 + (*_blossom_data)[vb].offset) {
2.1022 + _delta2->increase(vb, _blossom_set->classPrio(vb) -
2.1023 + (*_blossom_data)[vb].offset);
2.1024 + }
2.1025 + }
2.1026 + }
2.1027 + }
2.1028 + }
2.1029 + }
2.1030 + }
2.1031 +
2.1032 + void oddToMatched(int blossom) {
2.1033 + (*_blossom_data)[blossom].offset -= _delta_sum;
2.1034 +
2.1035 + if (_blossom_set->classPrio(blossom) !=
2.1036 + std::numeric_limits<Value>::max()) {
2.1037 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2.1038 + (*_blossom_data)[blossom].offset);
2.1039 + }
2.1040 +
2.1041 + if (!_blossom_set->trivial(blossom)) {
2.1042 + _delta4->erase(blossom);
2.1043 + }
2.1044 + }
2.1045 +
2.1046 + void oddToEven(int blossom, int tree) {
2.1047 + if (!_blossom_set->trivial(blossom)) {
2.1048 + _delta4->erase(blossom);
2.1049 + (*_blossom_data)[blossom].pot -=
2.1050 + 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
2.1051 + }
2.1052 +
2.1053 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2.1054 + n != INVALID; ++n) {
2.1055 + int ni = (*_node_index)[n];
2.1056 +
2.1057 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
2.1058 +
2.1059 + (*_node_data)[ni].heap.clear();
2.1060 + (*_node_data)[ni].heap_index.clear();
2.1061 + (*_node_data)[ni].pot +=
2.1062 + 2 * _delta_sum - (*_blossom_data)[blossom].offset;
2.1063 +
2.1064 + _delta1->push(n, (*_node_data)[ni].pot);
2.1065 +
2.1066 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
2.1067 + Node v = _graph.source(e);
2.1068 + int vb = _blossom_set->find(v);
2.1069 + int vi = (*_node_index)[v];
2.1070 +
2.1071 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2.1072 + dualScale * _weight[e];
2.1073 +
2.1074 + if ((*_blossom_data)[vb].status == EVEN) {
2.1075 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2.1076 + _delta3->push(e, rw / 2);
2.1077 + }
2.1078 + } else if ((*_blossom_data)[vb].status == UNMATCHED) {
2.1079 + if (_delta3->state(e) != _delta3->IN_HEAP) {
2.1080 + _delta3->push(e, rw);
2.1081 + }
2.1082 + } else {
2.1083 +
2.1084 + typename std::map<int, Arc>::iterator it =
2.1085 + (*_node_data)[vi].heap_index.find(tree);
2.1086 +
2.1087 + if (it != (*_node_data)[vi].heap_index.end()) {
2.1088 + if ((*_node_data)[vi].heap[it->second] > rw) {
2.1089 + (*_node_data)[vi].heap.replace(it->second, e);
2.1090 + (*_node_data)[vi].heap.decrease(e, rw);
2.1091 + it->second = e;
2.1092 + }
2.1093 + } else {
2.1094 + (*_node_data)[vi].heap.push(e, rw);
2.1095 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2.1096 + }
2.1097 +
2.1098 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2.1099 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2.1100 +
2.1101 + if ((*_blossom_data)[vb].status == MATCHED) {
2.1102 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
2.1103 + _delta2->push(vb, _blossom_set->classPrio(vb) -
2.1104 + (*_blossom_data)[vb].offset);
2.1105 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2.1106 + (*_blossom_data)[vb].offset) {
2.1107 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2.1108 + (*_blossom_data)[vb].offset);
2.1109 + }
2.1110 + }
2.1111 + }
2.1112 + }
2.1113 + }
2.1114 + }
2.1115 + (*_blossom_data)[blossom].offset = 0;
2.1116 + }
2.1117 +
2.1118 +
2.1119 + void matchedToUnmatched(int blossom) {
2.1120 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2.1121 + _delta2->erase(blossom);
2.1122 + }
2.1123 +
2.1124 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2.1125 + n != INVALID; ++n) {
2.1126 + int ni = (*_node_index)[n];
2.1127 +
2.1128 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
2.1129 +
2.1130 + (*_node_data)[ni].heap.clear();
2.1131 + (*_node_data)[ni].heap_index.clear();
2.1132 +
2.1133 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2.1134 + Node v = _graph.target(e);
2.1135 + int vb = _blossom_set->find(v);
2.1136 + int vi = (*_node_index)[v];
2.1137 +
2.1138 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2.1139 + dualScale * _weight[e];
2.1140 +
2.1141 + if ((*_blossom_data)[vb].status == EVEN) {
2.1142 + if (_delta3->state(e) != _delta3->IN_HEAP) {
2.1143 + _delta3->push(e, rw);
2.1144 + }
2.1145 + }
2.1146 + }
2.1147 + }
2.1148 + }
2.1149 +
2.1150 + void unmatchedToMatched(int blossom) {
2.1151 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2.1152 + n != INVALID; ++n) {
2.1153 + int ni = (*_node_index)[n];
2.1154 +
2.1155 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
2.1156 + Node v = _graph.source(e);
2.1157 + int vb = _blossom_set->find(v);
2.1158 + int vi = (*_node_index)[v];
2.1159 +
2.1160 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2.1161 + dualScale * _weight[e];
2.1162 +
2.1163 + if (vb == blossom) {
2.1164 + if (_delta3->state(e) == _delta3->IN_HEAP) {
2.1165 + _delta3->erase(e);
2.1166 + }
2.1167 + } else if ((*_blossom_data)[vb].status == EVEN) {
2.1168 +
2.1169 + if (_delta3->state(e) == _delta3->IN_HEAP) {
2.1170 + _delta3->erase(e);
2.1171 + }
2.1172 +
2.1173 + int vt = _tree_set->find(vb);
2.1174 +
2.1175 + Arc r = _graph.oppositeArc(e);
2.1176 +
2.1177 + typename std::map<int, Arc>::iterator it =
2.1178 + (*_node_data)[ni].heap_index.find(vt);
2.1179 +
2.1180 + if (it != (*_node_data)[ni].heap_index.end()) {
2.1181 + if ((*_node_data)[ni].heap[it->second] > rw) {
2.1182 + (*_node_data)[ni].heap.replace(it->second, r);
2.1183 + (*_node_data)[ni].heap.decrease(r, rw);
2.1184 + it->second = r;
2.1185 + }
2.1186 + } else {
2.1187 + (*_node_data)[ni].heap.push(r, rw);
2.1188 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
2.1189 + }
2.1190 +
2.1191 + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
2.1192 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
2.1193 +
2.1194 + if (_delta2->state(blossom) != _delta2->IN_HEAP) {
2.1195 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2.1196 + (*_blossom_data)[blossom].offset);
2.1197 + } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
2.1198 + (*_blossom_data)[blossom].offset){
2.1199 + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
2.1200 + (*_blossom_data)[blossom].offset);
2.1201 + }
2.1202 + }
2.1203 +
2.1204 + } else if ((*_blossom_data)[vb].status == UNMATCHED) {
2.1205 + if (_delta3->state(e) == _delta3->IN_HEAP) {
2.1206 + _delta3->erase(e);
2.1207 + }
2.1208 + }
2.1209 + }
2.1210 + }
2.1211 + }
2.1212 +
2.1213 + void alternatePath(int even, int tree) {
2.1214 + int odd;
2.1215 +
2.1216 + evenToMatched(even, tree);
2.1217 + (*_blossom_data)[even].status = MATCHED;
2.1218 +
2.1219 + while ((*_blossom_data)[even].pred != INVALID) {
2.1220 + odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
2.1221 + (*_blossom_data)[odd].status = MATCHED;
2.1222 + oddToMatched(odd);
2.1223 + (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
2.1224 +
2.1225 + even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
2.1226 + (*_blossom_data)[even].status = MATCHED;
2.1227 + evenToMatched(even, tree);
2.1228 + (*_blossom_data)[even].next =
2.1229 + _graph.oppositeArc((*_blossom_data)[odd].pred);
2.1230 + }
2.1231 +
2.1232 + }
2.1233 +
2.1234 + void destroyTree(int tree) {
2.1235 + for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
2.1236 + if ((*_blossom_data)[b].status == EVEN) {
2.1237 + (*_blossom_data)[b].status = MATCHED;
2.1238 + evenToMatched(b, tree);
2.1239 + } else if ((*_blossom_data)[b].status == ODD) {
2.1240 + (*_blossom_data)[b].status = MATCHED;
2.1241 + oddToMatched(b);
2.1242 + }
2.1243 + }
2.1244 + _tree_set->eraseClass(tree);
2.1245 + }
2.1246 +
2.1247 +
2.1248 + void unmatchNode(const Node& node) {
2.1249 + int blossom = _blossom_set->find(node);
2.1250 + int tree = _tree_set->find(blossom);
2.1251 +
2.1252 + alternatePath(blossom, tree);
2.1253 + destroyTree(tree);
2.1254 +
2.1255 + (*_blossom_data)[blossom].status = UNMATCHED;
2.1256 + (*_blossom_data)[blossom].base = node;
2.1257 + matchedToUnmatched(blossom);
2.1258 + }
2.1259 +
2.1260 +
2.1261 + void augmentOnEdge(const Edge& edge) {
2.1262 +
2.1263 + int left = _blossom_set->find(_graph.u(edge));
2.1264 + int right = _blossom_set->find(_graph.v(edge));
2.1265 +
2.1266 + if ((*_blossom_data)[left].status == EVEN) {
2.1267 + int left_tree = _tree_set->find(left);
2.1268 + alternatePath(left, left_tree);
2.1269 + destroyTree(left_tree);
2.1270 + } else {
2.1271 + (*_blossom_data)[left].status = MATCHED;
2.1272 + unmatchedToMatched(left);
2.1273 + }
2.1274 +
2.1275 + if ((*_blossom_data)[right].status == EVEN) {
2.1276 + int right_tree = _tree_set->find(right);
2.1277 + alternatePath(right, right_tree);
2.1278 + destroyTree(right_tree);
2.1279 + } else {
2.1280 + (*_blossom_data)[right].status = MATCHED;
2.1281 + unmatchedToMatched(right);
2.1282 + }
2.1283 +
2.1284 + (*_blossom_data)[left].next = _graph.direct(edge, true);
2.1285 + (*_blossom_data)[right].next = _graph.direct(edge, false);
2.1286 + }
2.1287 +
2.1288 + void extendOnArc(const Arc& arc) {
2.1289 + int base = _blossom_set->find(_graph.target(arc));
2.1290 + int tree = _tree_set->find(base);
2.1291 +
2.1292 + int odd = _blossom_set->find(_graph.source(arc));
2.1293 + _tree_set->insert(odd, tree);
2.1294 + (*_blossom_data)[odd].status = ODD;
2.1295 + matchedToOdd(odd);
2.1296 + (*_blossom_data)[odd].pred = arc;
2.1297 +
2.1298 + int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
2.1299 + (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
2.1300 + _tree_set->insert(even, tree);
2.1301 + (*_blossom_data)[even].status = EVEN;
2.1302 + matchedToEven(even, tree);
2.1303 + }
2.1304 +
2.1305 + void shrinkOnEdge(const Edge& edge, int tree) {
2.1306 + int nca = -1;
2.1307 + std::vector<int> left_path, right_path;
2.1308 +
2.1309 + {
2.1310 + std::set<int> left_set, right_set;
2.1311 + int left = _blossom_set->find(_graph.u(edge));
2.1312 + left_path.push_back(left);
2.1313 + left_set.insert(left);
2.1314 +
2.1315 + int right = _blossom_set->find(_graph.v(edge));
2.1316 + right_path.push_back(right);
2.1317 + right_set.insert(right);
2.1318 +
2.1319 + while (true) {
2.1320 +
2.1321 + if ((*_blossom_data)[left].pred == INVALID) break;
2.1322 +
2.1323 + left =
2.1324 + _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2.1325 + left_path.push_back(left);
2.1326 + left =
2.1327 + _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2.1328 + left_path.push_back(left);
2.1329 +
2.1330 + left_set.insert(left);
2.1331 +
2.1332 + if (right_set.find(left) != right_set.end()) {
2.1333 + nca = left;
2.1334 + break;
2.1335 + }
2.1336 +
2.1337 + if ((*_blossom_data)[right].pred == INVALID) break;
2.1338 +
2.1339 + right =
2.1340 + _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2.1341 + right_path.push_back(right);
2.1342 + right =
2.1343 + _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2.1344 + right_path.push_back(right);
2.1345 +
2.1346 + right_set.insert(right);
2.1347 +
2.1348 + if (left_set.find(right) != left_set.end()) {
2.1349 + nca = right;
2.1350 + break;
2.1351 + }
2.1352 +
2.1353 + }
2.1354 +
2.1355 + if (nca == -1) {
2.1356 + if ((*_blossom_data)[left].pred == INVALID) {
2.1357 + nca = right;
2.1358 + while (left_set.find(nca) == left_set.end()) {
2.1359 + nca =
2.1360 + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2.1361 + right_path.push_back(nca);
2.1362 + nca =
2.1363 + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2.1364 + right_path.push_back(nca);
2.1365 + }
2.1366 + } else {
2.1367 + nca = left;
2.1368 + while (right_set.find(nca) == right_set.end()) {
2.1369 + nca =
2.1370 + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2.1371 + left_path.push_back(nca);
2.1372 + nca =
2.1373 + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2.1374 + left_path.push_back(nca);
2.1375 + }
2.1376 + }
2.1377 + }
2.1378 + }
2.1379 +
2.1380 + std::vector<int> subblossoms;
2.1381 + Arc prev;
2.1382 +
2.1383 + prev = _graph.direct(edge, true);
2.1384 + for (int i = 0; left_path[i] != nca; i += 2) {
2.1385 + subblossoms.push_back(left_path[i]);
2.1386 + (*_blossom_data)[left_path[i]].next = prev;
2.1387 + _tree_set->erase(left_path[i]);
2.1388 +
2.1389 + subblossoms.push_back(left_path[i + 1]);
2.1390 + (*_blossom_data)[left_path[i + 1]].status = EVEN;
2.1391 + oddToEven(left_path[i + 1], tree);
2.1392 + _tree_set->erase(left_path[i + 1]);
2.1393 + prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
2.1394 + }
2.1395 +
2.1396 + int k = 0;
2.1397 + while (right_path[k] != nca) ++k;
2.1398 +
2.1399 + subblossoms.push_back(nca);
2.1400 + (*_blossom_data)[nca].next = prev;
2.1401 +
2.1402 + for (int i = k - 2; i >= 0; i -= 2) {
2.1403 + subblossoms.push_back(right_path[i + 1]);
2.1404 + (*_blossom_data)[right_path[i + 1]].status = EVEN;
2.1405 + oddToEven(right_path[i + 1], tree);
2.1406 + _tree_set->erase(right_path[i + 1]);
2.1407 +
2.1408 + (*_blossom_data)[right_path[i + 1]].next =
2.1409 + (*_blossom_data)[right_path[i + 1]].pred;
2.1410 +
2.1411 + subblossoms.push_back(right_path[i]);
2.1412 + _tree_set->erase(right_path[i]);
2.1413 + }
2.1414 +
2.1415 + int surface =
2.1416 + _blossom_set->join(subblossoms.begin(), subblossoms.end());
2.1417 +
2.1418 + for (int i = 0; i < int(subblossoms.size()); ++i) {
2.1419 + if (!_blossom_set->trivial(subblossoms[i])) {
2.1420 + (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
2.1421 + }
2.1422 + (*_blossom_data)[subblossoms[i]].status = MATCHED;
2.1423 + }
2.1424 +
2.1425 + (*_blossom_data)[surface].pot = -2 * _delta_sum;
2.1426 + (*_blossom_data)[surface].offset = 0;
2.1427 + (*_blossom_data)[surface].status = EVEN;
2.1428 + (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
2.1429 + (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
2.1430 +
2.1431 + _tree_set->insert(surface, tree);
2.1432 + _tree_set->erase(nca);
2.1433 + }
2.1434 +
2.1435 + void splitBlossom(int blossom) {
2.1436 + Arc next = (*_blossom_data)[blossom].next;
2.1437 + Arc pred = (*_blossom_data)[blossom].pred;
2.1438 +
2.1439 + int tree = _tree_set->find(blossom);
2.1440 +
2.1441 + (*_blossom_data)[blossom].status = MATCHED;
2.1442 + oddToMatched(blossom);
2.1443 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2.1444 + _delta2->erase(blossom);
2.1445 + }
2.1446 +
2.1447 + std::vector<int> subblossoms;
2.1448 + _blossom_set->split(blossom, std::back_inserter(subblossoms));
2.1449 +
2.1450 + Value offset = (*_blossom_data)[blossom].offset;
2.1451 + int b = _blossom_set->find(_graph.source(pred));
2.1452 + int d = _blossom_set->find(_graph.source(next));
2.1453 +
2.1454 + int ib = -1, id = -1;
2.1455 + for (int i = 0; i < int(subblossoms.size()); ++i) {
2.1456 + if (subblossoms[i] == b) ib = i;
2.1457 + if (subblossoms[i] == d) id = i;
2.1458 +
2.1459 + (*_blossom_data)[subblossoms[i]].offset = offset;
2.1460 + if (!_blossom_set->trivial(subblossoms[i])) {
2.1461 + (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
2.1462 + }
2.1463 + if (_blossom_set->classPrio(subblossoms[i]) !=
2.1464 + std::numeric_limits<Value>::max()) {
2.1465 + _delta2->push(subblossoms[i],
2.1466 + _blossom_set->classPrio(subblossoms[i]) -
2.1467 + (*_blossom_data)[subblossoms[i]].offset);
2.1468 + }
2.1469 + }
2.1470 +
2.1471 + if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
2.1472 + for (int i = (id + 1) % subblossoms.size();
2.1473 + i != ib; i = (i + 2) % subblossoms.size()) {
2.1474 + int sb = subblossoms[i];
2.1475 + int tb = subblossoms[(i + 1) % subblossoms.size()];
2.1476 + (*_blossom_data)[sb].next =
2.1477 + _graph.oppositeArc((*_blossom_data)[tb].next);
2.1478 + }
2.1479 +
2.1480 + for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
2.1481 + int sb = subblossoms[i];
2.1482 + int tb = subblossoms[(i + 1) % subblossoms.size()];
2.1483 + int ub = subblossoms[(i + 2) % subblossoms.size()];
2.1484 +
2.1485 + (*_blossom_data)[sb].status = ODD;
2.1486 + matchedToOdd(sb);
2.1487 + _tree_set->insert(sb, tree);
2.1488 + (*_blossom_data)[sb].pred = pred;
2.1489 + (*_blossom_data)[sb].next =
2.1490 + _graph.oppositeArc((*_blossom_data)[tb].next);
2.1491 +
2.1492 + pred = (*_blossom_data)[ub].next;
2.1493 +
2.1494 + (*_blossom_data)[tb].status = EVEN;
2.1495 + matchedToEven(tb, tree);
2.1496 + _tree_set->insert(tb, tree);
2.1497 + (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
2.1498 + }
2.1499 +
2.1500 + (*_blossom_data)[subblossoms[id]].status = ODD;
2.1501 + matchedToOdd(subblossoms[id]);
2.1502 + _tree_set->insert(subblossoms[id], tree);
2.1503 + (*_blossom_data)[subblossoms[id]].next = next;
2.1504 + (*_blossom_data)[subblossoms[id]].pred = pred;
2.1505 +
2.1506 + } else {
2.1507 +
2.1508 + for (int i = (ib + 1) % subblossoms.size();
2.1509 + i != id; i = (i + 2) % subblossoms.size()) {
2.1510 + int sb = subblossoms[i];
2.1511 + int tb = subblossoms[(i + 1) % subblossoms.size()];
2.1512 + (*_blossom_data)[sb].next =
2.1513 + _graph.oppositeArc((*_blossom_data)[tb].next);
2.1514 + }
2.1515 +
2.1516 + for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
2.1517 + int sb = subblossoms[i];
2.1518 + int tb = subblossoms[(i + 1) % subblossoms.size()];
2.1519 + int ub = subblossoms[(i + 2) % subblossoms.size()];
2.1520 +
2.1521 + (*_blossom_data)[sb].status = ODD;
2.1522 + matchedToOdd(sb);
2.1523 + _tree_set->insert(sb, tree);
2.1524 + (*_blossom_data)[sb].next = next;
2.1525 + (*_blossom_data)[sb].pred =
2.1526 + _graph.oppositeArc((*_blossom_data)[tb].next);
2.1527 +
2.1528 + (*_blossom_data)[tb].status = EVEN;
2.1529 + matchedToEven(tb, tree);
2.1530 + _tree_set->insert(tb, tree);
2.1531 + (*_blossom_data)[tb].pred =
2.1532 + (*_blossom_data)[tb].next =
2.1533 + _graph.oppositeArc((*_blossom_data)[ub].next);
2.1534 + next = (*_blossom_data)[ub].next;
2.1535 + }
2.1536 +
2.1537 + (*_blossom_data)[subblossoms[ib]].status = ODD;
2.1538 + matchedToOdd(subblossoms[ib]);
2.1539 + _tree_set->insert(subblossoms[ib], tree);
2.1540 + (*_blossom_data)[subblossoms[ib]].next = next;
2.1541 + (*_blossom_data)[subblossoms[ib]].pred = pred;
2.1542 + }
2.1543 + _tree_set->erase(blossom);
2.1544 + }
2.1545 +
2.1546 + void extractBlossom(int blossom, const Node& base, const Arc& matching) {
2.1547 + if (_blossom_set->trivial(blossom)) {
2.1548 + int bi = (*_node_index)[base];
2.1549 + Value pot = (*_node_data)[bi].pot;
2.1550 +
2.1551 + _matching->set(base, matching);
2.1552 + _blossom_node_list.push_back(base);
2.1553 + _node_potential->set(base, pot);
2.1554 + } else {
2.1555 +
2.1556 + Value pot = (*_blossom_data)[blossom].pot;
2.1557 + int bn = _blossom_node_list.size();
2.1558 +
2.1559 + std::vector<int> subblossoms;
2.1560 + _blossom_set->split(blossom, std::back_inserter(subblossoms));
2.1561 + int b = _blossom_set->find(base);
2.1562 + int ib = -1;
2.1563 + for (int i = 0; i < int(subblossoms.size()); ++i) {
2.1564 + if (subblossoms[i] == b) { ib = i; break; }
2.1565 + }
2.1566 +
2.1567 + for (int i = 1; i < int(subblossoms.size()); i += 2) {
2.1568 + int sb = subblossoms[(ib + i) % subblossoms.size()];
2.1569 + int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
2.1570 +
2.1571 + Arc m = (*_blossom_data)[tb].next;
2.1572 + extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
2.1573 + extractBlossom(tb, _graph.source(m), m);
2.1574 + }
2.1575 + extractBlossom(subblossoms[ib], base, matching);
2.1576 +
2.1577 + int en = _blossom_node_list.size();
2.1578 +
2.1579 + _blossom_potential.push_back(BlossomVariable(bn, en, pot));
2.1580 + }
2.1581 + }
2.1582 +
2.1583 + void extractMatching() {
2.1584 + std::vector<int> blossoms;
2.1585 + for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
2.1586 + blossoms.push_back(c);
2.1587 + }
2.1588 +
2.1589 + for (int i = 0; i < int(blossoms.size()); ++i) {
2.1590 + if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
2.1591 +
2.1592 + Value offset = (*_blossom_data)[blossoms[i]].offset;
2.1593 + (*_blossom_data)[blossoms[i]].pot += 2 * offset;
2.1594 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
2.1595 + n != INVALID; ++n) {
2.1596 + (*_node_data)[(*_node_index)[n]].pot -= offset;
2.1597 + }
2.1598 +
2.1599 + Arc matching = (*_blossom_data)[blossoms[i]].next;
2.1600 + Node base = _graph.source(matching);
2.1601 + extractBlossom(blossoms[i], base, matching);
2.1602 + } else {
2.1603 + Node base = (*_blossom_data)[blossoms[i]].base;
2.1604 + extractBlossom(blossoms[i], base, INVALID);
2.1605 + }
2.1606 + }
2.1607 + }
2.1608 +
2.1609 + public:
2.1610 +
2.1611 + /// \brief Constructor
2.1612 + ///
2.1613 + /// Constructor.
2.1614 + MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
2.1615 + : _graph(graph), _weight(weight), _matching(0),
2.1616 + _node_potential(0), _blossom_potential(), _blossom_node_list(),
2.1617 + _node_num(0), _blossom_num(0),
2.1618 +
2.1619 + _blossom_index(0), _blossom_set(0), _blossom_data(0),
2.1620 + _node_index(0), _node_heap_index(0), _node_data(0),
2.1621 + _tree_set_index(0), _tree_set(0),
2.1622 +
2.1623 + _delta1_index(0), _delta1(0),
2.1624 + _delta2_index(0), _delta2(0),
2.1625 + _delta3_index(0), _delta3(0),
2.1626 + _delta4_index(0), _delta4(0),
2.1627 +
2.1628 + _delta_sum() {}
2.1629 +
2.1630 + ~MaxWeightedMatching() {
2.1631 + destroyStructures();
2.1632 + }
2.1633 +
2.1634 + /// \name Execution control
2.1635 + /// The simplest way to execute the algorithm is to use the
2.1636 + /// \c run() member function.
2.1637 +
2.1638 + ///@{
2.1639 +
2.1640 + /// \brief Initialize the algorithm
2.1641 + ///
2.1642 + /// Initialize the algorithm
2.1643 + void init() {
2.1644 + createStructures();
2.1645 +
2.1646 + for (ArcIt e(_graph); e != INVALID; ++e) {
2.1647 + _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
2.1648 + }
2.1649 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.1650 + _delta1_index->set(n, _delta1->PRE_HEAP);
2.1651 + }
2.1652 + for (EdgeIt e(_graph); e != INVALID; ++e) {
2.1653 + _delta3_index->set(e, _delta3->PRE_HEAP);
2.1654 + }
2.1655 + for (int i = 0; i < _blossom_num; ++i) {
2.1656 + _delta2_index->set(i, _delta2->PRE_HEAP);
2.1657 + _delta4_index->set(i, _delta4->PRE_HEAP);
2.1658 + }
2.1659 +
2.1660 + int index = 0;
2.1661 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.1662 + Value max = 0;
2.1663 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2.1664 + if (_graph.target(e) == n) continue;
2.1665 + if ((dualScale * _weight[e]) / 2 > max) {
2.1666 + max = (dualScale * _weight[e]) / 2;
2.1667 + }
2.1668 + }
2.1669 + _node_index->set(n, index);
2.1670 + (*_node_data)[index].pot = max;
2.1671 + _delta1->push(n, max);
2.1672 + int blossom =
2.1673 + _blossom_set->insert(n, std::numeric_limits<Value>::max());
2.1674 +
2.1675 + _tree_set->insert(blossom);
2.1676 +
2.1677 + (*_blossom_data)[blossom].status = EVEN;
2.1678 + (*_blossom_data)[blossom].pred = INVALID;
2.1679 + (*_blossom_data)[blossom].next = INVALID;
2.1680 + (*_blossom_data)[blossom].pot = 0;
2.1681 + (*_blossom_data)[blossom].offset = 0;
2.1682 + ++index;
2.1683 + }
2.1684 + for (EdgeIt e(_graph); e != INVALID; ++e) {
2.1685 + int si = (*_node_index)[_graph.u(e)];
2.1686 + int ti = (*_node_index)[_graph.v(e)];
2.1687 + if (_graph.u(e) != _graph.v(e)) {
2.1688 + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
2.1689 + dualScale * _weight[e]) / 2);
2.1690 + }
2.1691 + }
2.1692 + }
2.1693 +
2.1694 + /// \brief Starts the algorithm
2.1695 + ///
2.1696 + /// Starts the algorithm
2.1697 + void start() {
2.1698 + enum OpType {
2.1699 + D1, D2, D3, D4
2.1700 + };
2.1701 +
2.1702 + int unmatched = _node_num;
2.1703 + while (unmatched > 0) {
2.1704 + Value d1 = !_delta1->empty() ?
2.1705 + _delta1->prio() : std::numeric_limits<Value>::max();
2.1706 +
2.1707 + Value d2 = !_delta2->empty() ?
2.1708 + _delta2->prio() : std::numeric_limits<Value>::max();
2.1709 +
2.1710 + Value d3 = !_delta3->empty() ?
2.1711 + _delta3->prio() : std::numeric_limits<Value>::max();
2.1712 +
2.1713 + Value d4 = !_delta4->empty() ?
2.1714 + _delta4->prio() : std::numeric_limits<Value>::max();
2.1715 +
2.1716 + _delta_sum = d1; OpType ot = D1;
2.1717 + if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
2.1718 + if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
2.1719 + if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
2.1720 +
2.1721 +
2.1722 + switch (ot) {
2.1723 + case D1:
2.1724 + {
2.1725 + Node n = _delta1->top();
2.1726 + unmatchNode(n);
2.1727 + --unmatched;
2.1728 + }
2.1729 + break;
2.1730 + case D2:
2.1731 + {
2.1732 + int blossom = _delta2->top();
2.1733 + Node n = _blossom_set->classTop(blossom);
2.1734 + Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
2.1735 + extendOnArc(e);
2.1736 + }
2.1737 + break;
2.1738 + case D3:
2.1739 + {
2.1740 + Edge e = _delta3->top();
2.1741 +
2.1742 + int left_blossom = _blossom_set->find(_graph.u(e));
2.1743 + int right_blossom = _blossom_set->find(_graph.v(e));
2.1744 +
2.1745 + if (left_blossom == right_blossom) {
2.1746 + _delta3->pop();
2.1747 + } else {
2.1748 + int left_tree;
2.1749 + if ((*_blossom_data)[left_blossom].status == EVEN) {
2.1750 + left_tree = _tree_set->find(left_blossom);
2.1751 + } else {
2.1752 + left_tree = -1;
2.1753 + ++unmatched;
2.1754 + }
2.1755 + int right_tree;
2.1756 + if ((*_blossom_data)[right_blossom].status == EVEN) {
2.1757 + right_tree = _tree_set->find(right_blossom);
2.1758 + } else {
2.1759 + right_tree = -1;
2.1760 + ++unmatched;
2.1761 + }
2.1762 +
2.1763 + if (left_tree == right_tree) {
2.1764 + shrinkOnEdge(e, left_tree);
2.1765 + } else {
2.1766 + augmentOnEdge(e);
2.1767 + unmatched -= 2;
2.1768 + }
2.1769 + }
2.1770 + } break;
2.1771 + case D4:
2.1772 + splitBlossom(_delta4->top());
2.1773 + break;
2.1774 + }
2.1775 + }
2.1776 + extractMatching();
2.1777 + }
2.1778 +
2.1779 + /// \brief Runs %MaxWeightedMatching algorithm.
2.1780 + ///
2.1781 + /// This method runs the %MaxWeightedMatching algorithm.
2.1782 + ///
2.1783 + /// \note mwm.run() is just a shortcut of the following code.
2.1784 + /// \code
2.1785 + /// mwm.init();
2.1786 + /// mwm.start();
2.1787 + /// \endcode
2.1788 + void run() {
2.1789 + init();
2.1790 + start();
2.1791 + }
2.1792 +
2.1793 + /// @}
2.1794 +
2.1795 + /// \name Primal solution
2.1796 + /// Functions to get the primal solution, ie. the matching.
2.1797 +
2.1798 + /// @{
2.1799 +
2.1800 + /// \brief Returns the weight of the matching.
2.1801 + ///
2.1802 + /// Returns the weight of the matching.
2.1803 + Value matchingValue() const {
2.1804 + Value sum = 0;
2.1805 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.1806 + if ((*_matching)[n] != INVALID) {
2.1807 + sum += _weight[(*_matching)[n]];
2.1808 + }
2.1809 + }
2.1810 + return sum /= 2;
2.1811 + }
2.1812 +
2.1813 + /// \brief Returns the cardinality of the matching.
2.1814 + ///
2.1815 + /// Returns the cardinality of the matching.
2.1816 + int matchingSize() const {
2.1817 + int num = 0;
2.1818 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.1819 + if ((*_matching)[n] != INVALID) {
2.1820 + ++num;
2.1821 + }
2.1822 + }
2.1823 + return num /= 2;
2.1824 + }
2.1825 +
2.1826 + /// \brief Returns true when the edge is in the matching.
2.1827 + ///
2.1828 + /// Returns true when the edge is in the matching.
2.1829 + bool matching(const Edge& edge) const {
2.1830 + return edge == (*_matching)[_graph.u(edge)];
2.1831 + }
2.1832 +
2.1833 + /// \brief Returns the incident matching arc.
2.1834 + ///
2.1835 + /// Returns the incident matching arc from given node. If the
2.1836 + /// node is not matched then it gives back \c INVALID.
2.1837 + Arc matching(const Node& node) const {
2.1838 + return (*_matching)[node];
2.1839 + }
2.1840 +
2.1841 + /// \brief Returns the mate of the node.
2.1842 + ///
2.1843 + /// Returns the adjancent node in a mathcing arc. If the node is
2.1844 + /// not matched then it gives back \c INVALID.
2.1845 + Node mate(const Node& node) const {
2.1846 + return (*_matching)[node] != INVALID ?
2.1847 + _graph.target((*_matching)[node]) : INVALID;
2.1848 + }
2.1849 +
2.1850 + /// @}
2.1851 +
2.1852 + /// \name Dual solution
2.1853 + /// Functions to get the dual solution.
2.1854 +
2.1855 + /// @{
2.1856 +
2.1857 + /// \brief Returns the value of the dual solution.
2.1858 + ///
2.1859 + /// Returns the value of the dual solution. It should be equal to
2.1860 + /// the primal value scaled by \ref dualScale "dual scale".
2.1861 + Value dualValue() const {
2.1862 + Value sum = 0;
2.1863 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.1864 + sum += nodeValue(n);
2.1865 + }
2.1866 + for (int i = 0; i < blossomNum(); ++i) {
2.1867 + sum += blossomValue(i) * (blossomSize(i) / 2);
2.1868 + }
2.1869 + return sum;
2.1870 + }
2.1871 +
2.1872 + /// \brief Returns the value of the node.
2.1873 + ///
2.1874 + /// Returns the the value of the node.
2.1875 + Value nodeValue(const Node& n) const {
2.1876 + return (*_node_potential)[n];
2.1877 + }
2.1878 +
2.1879 + /// \brief Returns the number of the blossoms in the basis.
2.1880 + ///
2.1881 + /// Returns the number of the blossoms in the basis.
2.1882 + /// \see BlossomIt
2.1883 + int blossomNum() const {
2.1884 + return _blossom_potential.size();
2.1885 + }
2.1886 +
2.1887 +
2.1888 + /// \brief Returns the number of the nodes in the blossom.
2.1889 + ///
2.1890 + /// Returns the number of the nodes in the blossom.
2.1891 + int blossomSize(int k) const {
2.1892 + return _blossom_potential[k].end - _blossom_potential[k].begin;
2.1893 + }
2.1894 +
2.1895 + /// \brief Returns the value of the blossom.
2.1896 + ///
2.1897 + /// Returns the the value of the blossom.
2.1898 + /// \see BlossomIt
2.1899 + Value blossomValue(int k) const {
2.1900 + return _blossom_potential[k].value;
2.1901 + }
2.1902 +
2.1903 + /// \brief Iterator for obtaining the nodes of the blossom.
2.1904 + ///
2.1905 + /// Iterator for obtaining the nodes of the blossom. This class
2.1906 + /// provides a common lemon style iterator for listing a
2.1907 + /// subset of the nodes.
2.1908 + class BlossomIt {
2.1909 + public:
2.1910 +
2.1911 + /// \brief Constructor.
2.1912 + ///
2.1913 + /// Constructor to get the nodes of the variable.
2.1914 + BlossomIt(const MaxWeightedMatching& algorithm, int variable)
2.1915 + : _algorithm(&algorithm)
2.1916 + {
2.1917 + _index = _algorithm->_blossom_potential[variable].begin;
2.1918 + _last = _algorithm->_blossom_potential[variable].end;
2.1919 + }
2.1920 +
2.1921 + /// \brief Conversion to node.
2.1922 + ///
2.1923 + /// Conversion to node.
2.1924 + operator Node() const {
2.1925 + return _algorithm->_blossom_node_list[_index];
2.1926 + }
2.1927 +
2.1928 + /// \brief Increment operator.
2.1929 + ///
2.1930 + /// Increment operator.
2.1931 + BlossomIt& operator++() {
2.1932 + ++_index;
2.1933 + return *this;
2.1934 + }
2.1935 +
2.1936 + /// \brief Validity checking
2.1937 + ///
2.1938 + /// Checks whether the iterator is invalid.
2.1939 + bool operator==(Invalid) const { return _index == _last; }
2.1940 +
2.1941 + /// \brief Validity checking
2.1942 + ///
2.1943 + /// Checks whether the iterator is valid.
2.1944 + bool operator!=(Invalid) const { return _index != _last; }
2.1945 +
2.1946 + private:
2.1947 + const MaxWeightedMatching* _algorithm;
2.1948 + int _last;
2.1949 + int _index;
2.1950 + };
2.1951 +
2.1952 + /// @}
2.1953 +
2.1954 + };
2.1955 +
2.1956 + /// \ingroup matching
2.1957 + ///
2.1958 + /// \brief Weighted perfect matching in general graphs
2.1959 + ///
2.1960 + /// This class provides an efficient implementation of Edmond's
2.1961 + /// maximum weighted perfect matching algorithm. The implementation
2.1962 + /// is based on extensive use of priority queues and provides
2.1963 + /// \f$O(nm\log(n))\f$ time complexity.
2.1964 + ///
2.1965 + /// The maximum weighted matching problem is to find undirected
2.1966 + /// edges in the graph with maximum overall weight and no two of
2.1967 + /// them shares their ends and covers all nodes. The problem can be
2.1968 + /// formulated with the following linear program.
2.1969 + /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
2.1970 + /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
2.1971 + \quad \forall B\in\mathcal{O}\f] */
2.1972 + /// \f[x_e \ge 0\quad \forall e\in E\f]
2.1973 + /// \f[\max \sum_{e\in E}x_ew_e\f]
2.1974 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
2.1975 + /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
2.1976 + /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
2.1977 + /// subsets of the nodes.
2.1978 + ///
2.1979 + /// The algorithm calculates an optimal matching and a proof of the
2.1980 + /// optimality. The solution of the dual problem can be used to check
2.1981 + /// the result of the algorithm. The dual linear problem is the
2.1982 + /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
2.1983 + w_{uv} \quad \forall uv\in E\f] */
2.1984 + /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
2.1985 + /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
2.1986 + \frac{\vert B \vert - 1}{2}z_B\f] */
2.1987 + ///
2.1988 + /// The algorithm can be executed with \c run() or the \c init() and
2.1989 + /// then the \c start() member functions. After it the matching can
2.1990 + /// be asked with \c matching() or mate() functions. The dual
2.1991 + /// solution can be get with \c nodeValue(), \c blossomNum() and \c
2.1992 + /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
2.1993 + /// "BlossomIt" nested class which is able to iterate on the nodes
2.1994 + /// of a blossom. If the value type is integral then the dual
2.1995 + /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
2.1996 + template <typename _Graph,
2.1997 + typename _WeightMap = typename _Graph::template EdgeMap<int> >
2.1998 + class MaxWeightedPerfectMatching {
2.1999 + public:
2.2000 +
2.2001 + typedef _Graph Graph;
2.2002 + typedef _WeightMap WeightMap;
2.2003 + typedef typename WeightMap::Value Value;
2.2004 +
2.2005 + /// \brief Scaling factor for dual solution
2.2006 + ///
2.2007 + /// Scaling factor for dual solution, it is equal to 4 or 1
2.2008 + /// according to the value type.
2.2009 + static const int dualScale =
2.2010 + std::numeric_limits<Value>::is_integer ? 4 : 1;
2.2011 +
2.2012 + typedef typename Graph::template NodeMap<typename Graph::Arc>
2.2013 + MatchingMap;
2.2014 +
2.2015 + private:
2.2016 +
2.2017 + TEMPLATE_GRAPH_TYPEDEFS(Graph);
2.2018 +
2.2019 + typedef typename Graph::template NodeMap<Value> NodePotential;
2.2020 + typedef std::vector<Node> BlossomNodeList;
2.2021 +
2.2022 + struct BlossomVariable {
2.2023 + int begin, end;
2.2024 + Value value;
2.2025 +
2.2026 + BlossomVariable(int _begin, int _end, Value _value)
2.2027 + : begin(_begin), end(_end), value(_value) {}
2.2028 +
2.2029 + };
2.2030 +
2.2031 + typedef std::vector<BlossomVariable> BlossomPotential;
2.2032 +
2.2033 + const Graph& _graph;
2.2034 + const WeightMap& _weight;
2.2035 +
2.2036 + MatchingMap* _matching;
2.2037 +
2.2038 + NodePotential* _node_potential;
2.2039 +
2.2040 + BlossomPotential _blossom_potential;
2.2041 + BlossomNodeList _blossom_node_list;
2.2042 +
2.2043 + int _node_num;
2.2044 + int _blossom_num;
2.2045 +
2.2046 + typedef RangeMap<int> IntIntMap;
2.2047 +
2.2048 + enum Status {
2.2049 + EVEN = -1, MATCHED = 0, ODD = 1
2.2050 + };
2.2051 +
2.2052 + typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
2.2053 + struct BlossomData {
2.2054 + int tree;
2.2055 + Status status;
2.2056 + Arc pred, next;
2.2057 + Value pot, offset;
2.2058 + };
2.2059 +
2.2060 + IntNodeMap *_blossom_index;
2.2061 + BlossomSet *_blossom_set;
2.2062 + RangeMap<BlossomData>* _blossom_data;
2.2063 +
2.2064 + IntNodeMap *_node_index;
2.2065 + IntArcMap *_node_heap_index;
2.2066 +
2.2067 + struct NodeData {
2.2068 +
2.2069 + NodeData(IntArcMap& node_heap_index)
2.2070 + : heap(node_heap_index) {}
2.2071 +
2.2072 + int blossom;
2.2073 + Value pot;
2.2074 + BinHeap<Value, IntArcMap> heap;
2.2075 + std::map<int, Arc> heap_index;
2.2076 +
2.2077 + int tree;
2.2078 + };
2.2079 +
2.2080 + RangeMap<NodeData>* _node_data;
2.2081 +
2.2082 + typedef ExtendFindEnum<IntIntMap> TreeSet;
2.2083 +
2.2084 + IntIntMap *_tree_set_index;
2.2085 + TreeSet *_tree_set;
2.2086 +
2.2087 + IntIntMap *_delta2_index;
2.2088 + BinHeap<Value, IntIntMap> *_delta2;
2.2089 +
2.2090 + IntEdgeMap *_delta3_index;
2.2091 + BinHeap<Value, IntEdgeMap> *_delta3;
2.2092 +
2.2093 + IntIntMap *_delta4_index;
2.2094 + BinHeap<Value, IntIntMap> *_delta4;
2.2095 +
2.2096 + Value _delta_sum;
2.2097 +
2.2098 + void createStructures() {
2.2099 + _node_num = countNodes(_graph);
2.2100 + _blossom_num = _node_num * 3 / 2;
2.2101 +
2.2102 + if (!_matching) {
2.2103 + _matching = new MatchingMap(_graph);
2.2104 + }
2.2105 + if (!_node_potential) {
2.2106 + _node_potential = new NodePotential(_graph);
2.2107 + }
2.2108 + if (!_blossom_set) {
2.2109 + _blossom_index = new IntNodeMap(_graph);
2.2110 + _blossom_set = new BlossomSet(*_blossom_index);
2.2111 + _blossom_data = new RangeMap<BlossomData>(_blossom_num);
2.2112 + }
2.2113 +
2.2114 + if (!_node_index) {
2.2115 + _node_index = new IntNodeMap(_graph);
2.2116 + _node_heap_index = new IntArcMap(_graph);
2.2117 + _node_data = new RangeMap<NodeData>(_node_num,
2.2118 + NodeData(*_node_heap_index));
2.2119 + }
2.2120 +
2.2121 + if (!_tree_set) {
2.2122 + _tree_set_index = new IntIntMap(_blossom_num);
2.2123 + _tree_set = new TreeSet(*_tree_set_index);
2.2124 + }
2.2125 + if (!_delta2) {
2.2126 + _delta2_index = new IntIntMap(_blossom_num);
2.2127 + _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
2.2128 + }
2.2129 + if (!_delta3) {
2.2130 + _delta3_index = new IntEdgeMap(_graph);
2.2131 + _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
2.2132 + }
2.2133 + if (!_delta4) {
2.2134 + _delta4_index = new IntIntMap(_blossom_num);
2.2135 + _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
2.2136 + }
2.2137 + }
2.2138 +
2.2139 + void destroyStructures() {
2.2140 + _node_num = countNodes(_graph);
2.2141 + _blossom_num = _node_num * 3 / 2;
2.2142 +
2.2143 + if (_matching) {
2.2144 + delete _matching;
2.2145 + }
2.2146 + if (_node_potential) {
2.2147 + delete _node_potential;
2.2148 + }
2.2149 + if (_blossom_set) {
2.2150 + delete _blossom_index;
2.2151 + delete _blossom_set;
2.2152 + delete _blossom_data;
2.2153 + }
2.2154 +
2.2155 + if (_node_index) {
2.2156 + delete _node_index;
2.2157 + delete _node_heap_index;
2.2158 + delete _node_data;
2.2159 + }
2.2160 +
2.2161 + if (_tree_set) {
2.2162 + delete _tree_set_index;
2.2163 + delete _tree_set;
2.2164 + }
2.2165 + if (_delta2) {
2.2166 + delete _delta2_index;
2.2167 + delete _delta2;
2.2168 + }
2.2169 + if (_delta3) {
2.2170 + delete _delta3_index;
2.2171 + delete _delta3;
2.2172 + }
2.2173 + if (_delta4) {
2.2174 + delete _delta4_index;
2.2175 + delete _delta4;
2.2176 + }
2.2177 + }
2.2178 +
2.2179 + void matchedToEven(int blossom, int tree) {
2.2180 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2.2181 + _delta2->erase(blossom);
2.2182 + }
2.2183 +
2.2184 + if (!_blossom_set->trivial(blossom)) {
2.2185 + (*_blossom_data)[blossom].pot -=
2.2186 + 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
2.2187 + }
2.2188 +
2.2189 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2.2190 + n != INVALID; ++n) {
2.2191 +
2.2192 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
2.2193 + int ni = (*_node_index)[n];
2.2194 +
2.2195 + (*_node_data)[ni].heap.clear();
2.2196 + (*_node_data)[ni].heap_index.clear();
2.2197 +
2.2198 + (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
2.2199 +
2.2200 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
2.2201 + Node v = _graph.source(e);
2.2202 + int vb = _blossom_set->find(v);
2.2203 + int vi = (*_node_index)[v];
2.2204 +
2.2205 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2.2206 + dualScale * _weight[e];
2.2207 +
2.2208 + if ((*_blossom_data)[vb].status == EVEN) {
2.2209 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2.2210 + _delta3->push(e, rw / 2);
2.2211 + }
2.2212 + } else {
2.2213 + typename std::map<int, Arc>::iterator it =
2.2214 + (*_node_data)[vi].heap_index.find(tree);
2.2215 +
2.2216 + if (it != (*_node_data)[vi].heap_index.end()) {
2.2217 + if ((*_node_data)[vi].heap[it->second] > rw) {
2.2218 + (*_node_data)[vi].heap.replace(it->second, e);
2.2219 + (*_node_data)[vi].heap.decrease(e, rw);
2.2220 + it->second = e;
2.2221 + }
2.2222 + } else {
2.2223 + (*_node_data)[vi].heap.push(e, rw);
2.2224 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2.2225 + }
2.2226 +
2.2227 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2.2228 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2.2229 +
2.2230 + if ((*_blossom_data)[vb].status == MATCHED) {
2.2231 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
2.2232 + _delta2->push(vb, _blossom_set->classPrio(vb) -
2.2233 + (*_blossom_data)[vb].offset);
2.2234 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2.2235 + (*_blossom_data)[vb].offset){
2.2236 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2.2237 + (*_blossom_data)[vb].offset);
2.2238 + }
2.2239 + }
2.2240 + }
2.2241 + }
2.2242 + }
2.2243 + }
2.2244 + (*_blossom_data)[blossom].offset = 0;
2.2245 + }
2.2246 +
2.2247 + void matchedToOdd(int blossom) {
2.2248 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2.2249 + _delta2->erase(blossom);
2.2250 + }
2.2251 + (*_blossom_data)[blossom].offset += _delta_sum;
2.2252 + if (!_blossom_set->trivial(blossom)) {
2.2253 + _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
2.2254 + (*_blossom_data)[blossom].offset);
2.2255 + }
2.2256 + }
2.2257 +
2.2258 + void evenToMatched(int blossom, int tree) {
2.2259 + if (!_blossom_set->trivial(blossom)) {
2.2260 + (*_blossom_data)[blossom].pot += 2 * _delta_sum;
2.2261 + }
2.2262 +
2.2263 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2.2264 + n != INVALID; ++n) {
2.2265 + int ni = (*_node_index)[n];
2.2266 + (*_node_data)[ni].pot -= _delta_sum;
2.2267 +
2.2268 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
2.2269 + Node v = _graph.source(e);
2.2270 + int vb = _blossom_set->find(v);
2.2271 + int vi = (*_node_index)[v];
2.2272 +
2.2273 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2.2274 + dualScale * _weight[e];
2.2275 +
2.2276 + if (vb == blossom) {
2.2277 + if (_delta3->state(e) == _delta3->IN_HEAP) {
2.2278 + _delta3->erase(e);
2.2279 + }
2.2280 + } else if ((*_blossom_data)[vb].status == EVEN) {
2.2281 +
2.2282 + if (_delta3->state(e) == _delta3->IN_HEAP) {
2.2283 + _delta3->erase(e);
2.2284 + }
2.2285 +
2.2286 + int vt = _tree_set->find(vb);
2.2287 +
2.2288 + if (vt != tree) {
2.2289 +
2.2290 + Arc r = _graph.oppositeArc(e);
2.2291 +
2.2292 + typename std::map<int, Arc>::iterator it =
2.2293 + (*_node_data)[ni].heap_index.find(vt);
2.2294 +
2.2295 + if (it != (*_node_data)[ni].heap_index.end()) {
2.2296 + if ((*_node_data)[ni].heap[it->second] > rw) {
2.2297 + (*_node_data)[ni].heap.replace(it->second, r);
2.2298 + (*_node_data)[ni].heap.decrease(r, rw);
2.2299 + it->second = r;
2.2300 + }
2.2301 + } else {
2.2302 + (*_node_data)[ni].heap.push(r, rw);
2.2303 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
2.2304 + }
2.2305 +
2.2306 + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
2.2307 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
2.2308 +
2.2309 + if (_delta2->state(blossom) != _delta2->IN_HEAP) {
2.2310 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2.2311 + (*_blossom_data)[blossom].offset);
2.2312 + } else if ((*_delta2)[blossom] >
2.2313 + _blossom_set->classPrio(blossom) -
2.2314 + (*_blossom_data)[blossom].offset){
2.2315 + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
2.2316 + (*_blossom_data)[blossom].offset);
2.2317 + }
2.2318 + }
2.2319 + }
2.2320 + } else {
2.2321 +
2.2322 + typename std::map<int, Arc>::iterator it =
2.2323 + (*_node_data)[vi].heap_index.find(tree);
2.2324 +
2.2325 + if (it != (*_node_data)[vi].heap_index.end()) {
2.2326 + (*_node_data)[vi].heap.erase(it->second);
2.2327 + (*_node_data)[vi].heap_index.erase(it);
2.2328 + if ((*_node_data)[vi].heap.empty()) {
2.2329 + _blossom_set->increase(v, std::numeric_limits<Value>::max());
2.2330 + } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
2.2331 + _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
2.2332 + }
2.2333 +
2.2334 + if ((*_blossom_data)[vb].status == MATCHED) {
2.2335 + if (_blossom_set->classPrio(vb) ==
2.2336 + std::numeric_limits<Value>::max()) {
2.2337 + _delta2->erase(vb);
2.2338 + } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
2.2339 + (*_blossom_data)[vb].offset) {
2.2340 + _delta2->increase(vb, _blossom_set->classPrio(vb) -
2.2341 + (*_blossom_data)[vb].offset);
2.2342 + }
2.2343 + }
2.2344 + }
2.2345 + }
2.2346 + }
2.2347 + }
2.2348 + }
2.2349 +
2.2350 + void oddToMatched(int blossom) {
2.2351 + (*_blossom_data)[blossom].offset -= _delta_sum;
2.2352 +
2.2353 + if (_blossom_set->classPrio(blossom) !=
2.2354 + std::numeric_limits<Value>::max()) {
2.2355 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
2.2356 + (*_blossom_data)[blossom].offset);
2.2357 + }
2.2358 +
2.2359 + if (!_blossom_set->trivial(blossom)) {
2.2360 + _delta4->erase(blossom);
2.2361 + }
2.2362 + }
2.2363 +
2.2364 + void oddToEven(int blossom, int tree) {
2.2365 + if (!_blossom_set->trivial(blossom)) {
2.2366 + _delta4->erase(blossom);
2.2367 + (*_blossom_data)[blossom].pot -=
2.2368 + 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
2.2369 + }
2.2370 +
2.2371 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
2.2372 + n != INVALID; ++n) {
2.2373 + int ni = (*_node_index)[n];
2.2374 +
2.2375 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
2.2376 +
2.2377 + (*_node_data)[ni].heap.clear();
2.2378 + (*_node_data)[ni].heap_index.clear();
2.2379 + (*_node_data)[ni].pot +=
2.2380 + 2 * _delta_sum - (*_blossom_data)[blossom].offset;
2.2381 +
2.2382 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
2.2383 + Node v = _graph.source(e);
2.2384 + int vb = _blossom_set->find(v);
2.2385 + int vi = (*_node_index)[v];
2.2386 +
2.2387 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
2.2388 + dualScale * _weight[e];
2.2389 +
2.2390 + if ((*_blossom_data)[vb].status == EVEN) {
2.2391 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
2.2392 + _delta3->push(e, rw / 2);
2.2393 + }
2.2394 + } else {
2.2395 +
2.2396 + typename std::map<int, Arc>::iterator it =
2.2397 + (*_node_data)[vi].heap_index.find(tree);
2.2398 +
2.2399 + if (it != (*_node_data)[vi].heap_index.end()) {
2.2400 + if ((*_node_data)[vi].heap[it->second] > rw) {
2.2401 + (*_node_data)[vi].heap.replace(it->second, e);
2.2402 + (*_node_data)[vi].heap.decrease(e, rw);
2.2403 + it->second = e;
2.2404 + }
2.2405 + } else {
2.2406 + (*_node_data)[vi].heap.push(e, rw);
2.2407 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
2.2408 + }
2.2409 +
2.2410 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
2.2411 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
2.2412 +
2.2413 + if ((*_blossom_data)[vb].status == MATCHED) {
2.2414 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
2.2415 + _delta2->push(vb, _blossom_set->classPrio(vb) -
2.2416 + (*_blossom_data)[vb].offset);
2.2417 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
2.2418 + (*_blossom_data)[vb].offset) {
2.2419 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
2.2420 + (*_blossom_data)[vb].offset);
2.2421 + }
2.2422 + }
2.2423 + }
2.2424 + }
2.2425 + }
2.2426 + }
2.2427 + (*_blossom_data)[blossom].offset = 0;
2.2428 + }
2.2429 +
2.2430 + void alternatePath(int even, int tree) {
2.2431 + int odd;
2.2432 +
2.2433 + evenToMatched(even, tree);
2.2434 + (*_blossom_data)[even].status = MATCHED;
2.2435 +
2.2436 + while ((*_blossom_data)[even].pred != INVALID) {
2.2437 + odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
2.2438 + (*_blossom_data)[odd].status = MATCHED;
2.2439 + oddToMatched(odd);
2.2440 + (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
2.2441 +
2.2442 + even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
2.2443 + (*_blossom_data)[even].status = MATCHED;
2.2444 + evenToMatched(even, tree);
2.2445 + (*_blossom_data)[even].next =
2.2446 + _graph.oppositeArc((*_blossom_data)[odd].pred);
2.2447 + }
2.2448 +
2.2449 + }
2.2450 +
2.2451 + void destroyTree(int tree) {
2.2452 + for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
2.2453 + if ((*_blossom_data)[b].status == EVEN) {
2.2454 + (*_blossom_data)[b].status = MATCHED;
2.2455 + evenToMatched(b, tree);
2.2456 + } else if ((*_blossom_data)[b].status == ODD) {
2.2457 + (*_blossom_data)[b].status = MATCHED;
2.2458 + oddToMatched(b);
2.2459 + }
2.2460 + }
2.2461 + _tree_set->eraseClass(tree);
2.2462 + }
2.2463 +
2.2464 + void augmentOnEdge(const Edge& edge) {
2.2465 +
2.2466 + int left = _blossom_set->find(_graph.u(edge));
2.2467 + int right = _blossom_set->find(_graph.v(edge));
2.2468 +
2.2469 + int left_tree = _tree_set->find(left);
2.2470 + alternatePath(left, left_tree);
2.2471 + destroyTree(left_tree);
2.2472 +
2.2473 + int right_tree = _tree_set->find(right);
2.2474 + alternatePath(right, right_tree);
2.2475 + destroyTree(right_tree);
2.2476 +
2.2477 + (*_blossom_data)[left].next = _graph.direct(edge, true);
2.2478 + (*_blossom_data)[right].next = _graph.direct(edge, false);
2.2479 + }
2.2480 +
2.2481 + void extendOnArc(const Arc& arc) {
2.2482 + int base = _blossom_set->find(_graph.target(arc));
2.2483 + int tree = _tree_set->find(base);
2.2484 +
2.2485 + int odd = _blossom_set->find(_graph.source(arc));
2.2486 + _tree_set->insert(odd, tree);
2.2487 + (*_blossom_data)[odd].status = ODD;
2.2488 + matchedToOdd(odd);
2.2489 + (*_blossom_data)[odd].pred = arc;
2.2490 +
2.2491 + int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
2.2492 + (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
2.2493 + _tree_set->insert(even, tree);
2.2494 + (*_blossom_data)[even].status = EVEN;
2.2495 + matchedToEven(even, tree);
2.2496 + }
2.2497 +
2.2498 + void shrinkOnEdge(const Edge& edge, int tree) {
2.2499 + int nca = -1;
2.2500 + std::vector<int> left_path, right_path;
2.2501 +
2.2502 + {
2.2503 + std::set<int> left_set, right_set;
2.2504 + int left = _blossom_set->find(_graph.u(edge));
2.2505 + left_path.push_back(left);
2.2506 + left_set.insert(left);
2.2507 +
2.2508 + int right = _blossom_set->find(_graph.v(edge));
2.2509 + right_path.push_back(right);
2.2510 + right_set.insert(right);
2.2511 +
2.2512 + while (true) {
2.2513 +
2.2514 + if ((*_blossom_data)[left].pred == INVALID) break;
2.2515 +
2.2516 + left =
2.2517 + _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2.2518 + left_path.push_back(left);
2.2519 + left =
2.2520 + _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
2.2521 + left_path.push_back(left);
2.2522 +
2.2523 + left_set.insert(left);
2.2524 +
2.2525 + if (right_set.find(left) != right_set.end()) {
2.2526 + nca = left;
2.2527 + break;
2.2528 + }
2.2529 +
2.2530 + if ((*_blossom_data)[right].pred == INVALID) break;
2.2531 +
2.2532 + right =
2.2533 + _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2.2534 + right_path.push_back(right);
2.2535 + right =
2.2536 + _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
2.2537 + right_path.push_back(right);
2.2538 +
2.2539 + right_set.insert(right);
2.2540 +
2.2541 + if (left_set.find(right) != left_set.end()) {
2.2542 + nca = right;
2.2543 + break;
2.2544 + }
2.2545 +
2.2546 + }
2.2547 +
2.2548 + if (nca == -1) {
2.2549 + if ((*_blossom_data)[left].pred == INVALID) {
2.2550 + nca = right;
2.2551 + while (left_set.find(nca) == left_set.end()) {
2.2552 + nca =
2.2553 + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2.2554 + right_path.push_back(nca);
2.2555 + nca =
2.2556 + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2.2557 + right_path.push_back(nca);
2.2558 + }
2.2559 + } else {
2.2560 + nca = left;
2.2561 + while (right_set.find(nca) == right_set.end()) {
2.2562 + nca =
2.2563 + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2.2564 + left_path.push_back(nca);
2.2565 + nca =
2.2566 + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
2.2567 + left_path.push_back(nca);
2.2568 + }
2.2569 + }
2.2570 + }
2.2571 + }
2.2572 +
2.2573 + std::vector<int> subblossoms;
2.2574 + Arc prev;
2.2575 +
2.2576 + prev = _graph.direct(edge, true);
2.2577 + for (int i = 0; left_path[i] != nca; i += 2) {
2.2578 + subblossoms.push_back(left_path[i]);
2.2579 + (*_blossom_data)[left_path[i]].next = prev;
2.2580 + _tree_set->erase(left_path[i]);
2.2581 +
2.2582 + subblossoms.push_back(left_path[i + 1]);
2.2583 + (*_blossom_data)[left_path[i + 1]].status = EVEN;
2.2584 + oddToEven(left_path[i + 1], tree);
2.2585 + _tree_set->erase(left_path[i + 1]);
2.2586 + prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
2.2587 + }
2.2588 +
2.2589 + int k = 0;
2.2590 + while (right_path[k] != nca) ++k;
2.2591 +
2.2592 + subblossoms.push_back(nca);
2.2593 + (*_blossom_data)[nca].next = prev;
2.2594 +
2.2595 + for (int i = k - 2; i >= 0; i -= 2) {
2.2596 + subblossoms.push_back(right_path[i + 1]);
2.2597 + (*_blossom_data)[right_path[i + 1]].status = EVEN;
2.2598 + oddToEven(right_path[i + 1], tree);
2.2599 + _tree_set->erase(right_path[i + 1]);
2.2600 +
2.2601 + (*_blossom_data)[right_path[i + 1]].next =
2.2602 + (*_blossom_data)[right_path[i + 1]].pred;
2.2603 +
2.2604 + subblossoms.push_back(right_path[i]);
2.2605 + _tree_set->erase(right_path[i]);
2.2606 + }
2.2607 +
2.2608 + int surface =
2.2609 + _blossom_set->join(subblossoms.begin(), subblossoms.end());
2.2610 +
2.2611 + for (int i = 0; i < int(subblossoms.size()); ++i) {
2.2612 + if (!_blossom_set->trivial(subblossoms[i])) {
2.2613 + (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
2.2614 + }
2.2615 + (*_blossom_data)[subblossoms[i]].status = MATCHED;
2.2616 + }
2.2617 +
2.2618 + (*_blossom_data)[surface].pot = -2 * _delta_sum;
2.2619 + (*_blossom_data)[surface].offset = 0;
2.2620 + (*_blossom_data)[surface].status = EVEN;
2.2621 + (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
2.2622 + (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
2.2623 +
2.2624 + _tree_set->insert(surface, tree);
2.2625 + _tree_set->erase(nca);
2.2626 + }
2.2627 +
2.2628 + void splitBlossom(int blossom) {
2.2629 + Arc next = (*_blossom_data)[blossom].next;
2.2630 + Arc pred = (*_blossom_data)[blossom].pred;
2.2631 +
2.2632 + int tree = _tree_set->find(blossom);
2.2633 +
2.2634 + (*_blossom_data)[blossom].status = MATCHED;
2.2635 + oddToMatched(blossom);
2.2636 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
2.2637 + _delta2->erase(blossom);
2.2638 + }
2.2639 +
2.2640 + std::vector<int> subblossoms;
2.2641 + _blossom_set->split(blossom, std::back_inserter(subblossoms));
2.2642 +
2.2643 + Value offset = (*_blossom_data)[blossom].offset;
2.2644 + int b = _blossom_set->find(_graph.source(pred));
2.2645 + int d = _blossom_set->find(_graph.source(next));
2.2646 +
2.2647 + int ib = -1, id = -1;
2.2648 + for (int i = 0; i < int(subblossoms.size()); ++i) {
2.2649 + if (subblossoms[i] == b) ib = i;
2.2650 + if (subblossoms[i] == d) id = i;
2.2651 +
2.2652 + (*_blossom_data)[subblossoms[i]].offset = offset;
2.2653 + if (!_blossom_set->trivial(subblossoms[i])) {
2.2654 + (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
2.2655 + }
2.2656 + if (_blossom_set->classPrio(subblossoms[i]) !=
2.2657 + std::numeric_limits<Value>::max()) {
2.2658 + _delta2->push(subblossoms[i],
2.2659 + _blossom_set->classPrio(subblossoms[i]) -
2.2660 + (*_blossom_data)[subblossoms[i]].offset);
2.2661 + }
2.2662 + }
2.2663 +
2.2664 + if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
2.2665 + for (int i = (id + 1) % subblossoms.size();
2.2666 + i != ib; i = (i + 2) % subblossoms.size()) {
2.2667 + int sb = subblossoms[i];
2.2668 + int tb = subblossoms[(i + 1) % subblossoms.size()];
2.2669 + (*_blossom_data)[sb].next =
2.2670 + _graph.oppositeArc((*_blossom_data)[tb].next);
2.2671 + }
2.2672 +
2.2673 + for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
2.2674 + int sb = subblossoms[i];
2.2675 + int tb = subblossoms[(i + 1) % subblossoms.size()];
2.2676 + int ub = subblossoms[(i + 2) % subblossoms.size()];
2.2677 +
2.2678 + (*_blossom_data)[sb].status = ODD;
2.2679 + matchedToOdd(sb);
2.2680 + _tree_set->insert(sb, tree);
2.2681 + (*_blossom_data)[sb].pred = pred;
2.2682 + (*_blossom_data)[sb].next =
2.2683 + _graph.oppositeArc((*_blossom_data)[tb].next);
2.2684 +
2.2685 + pred = (*_blossom_data)[ub].next;
2.2686 +
2.2687 + (*_blossom_data)[tb].status = EVEN;
2.2688 + matchedToEven(tb, tree);
2.2689 + _tree_set->insert(tb, tree);
2.2690 + (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
2.2691 + }
2.2692 +
2.2693 + (*_blossom_data)[subblossoms[id]].status = ODD;
2.2694 + matchedToOdd(subblossoms[id]);
2.2695 + _tree_set->insert(subblossoms[id], tree);
2.2696 + (*_blossom_data)[subblossoms[id]].next = next;
2.2697 + (*_blossom_data)[subblossoms[id]].pred = pred;
2.2698 +
2.2699 + } else {
2.2700 +
2.2701 + for (int i = (ib + 1) % subblossoms.size();
2.2702 + i != id; i = (i + 2) % subblossoms.size()) {
2.2703 + int sb = subblossoms[i];
2.2704 + int tb = subblossoms[(i + 1) % subblossoms.size()];
2.2705 + (*_blossom_data)[sb].next =
2.2706 + _graph.oppositeArc((*_blossom_data)[tb].next);
2.2707 + }
2.2708 +
2.2709 + for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
2.2710 + int sb = subblossoms[i];
2.2711 + int tb = subblossoms[(i + 1) % subblossoms.size()];
2.2712 + int ub = subblossoms[(i + 2) % subblossoms.size()];
2.2713 +
2.2714 + (*_blossom_data)[sb].status = ODD;
2.2715 + matchedToOdd(sb);
2.2716 + _tree_set->insert(sb, tree);
2.2717 + (*_blossom_data)[sb].next = next;
2.2718 + (*_blossom_data)[sb].pred =
2.2719 + _graph.oppositeArc((*_blossom_data)[tb].next);
2.2720 +
2.2721 + (*_blossom_data)[tb].status = EVEN;
2.2722 + matchedToEven(tb, tree);
2.2723 + _tree_set->insert(tb, tree);
2.2724 + (*_blossom_data)[tb].pred =
2.2725 + (*_blossom_data)[tb].next =
2.2726 + _graph.oppositeArc((*_blossom_data)[ub].next);
2.2727 + next = (*_blossom_data)[ub].next;
2.2728 + }
2.2729 +
2.2730 + (*_blossom_data)[subblossoms[ib]].status = ODD;
2.2731 + matchedToOdd(subblossoms[ib]);
2.2732 + _tree_set->insert(subblossoms[ib], tree);
2.2733 + (*_blossom_data)[subblossoms[ib]].next = next;
2.2734 + (*_blossom_data)[subblossoms[ib]].pred = pred;
2.2735 + }
2.2736 + _tree_set->erase(blossom);
2.2737 + }
2.2738 +
2.2739 + void extractBlossom(int blossom, const Node& base, const Arc& matching) {
2.2740 + if (_blossom_set->trivial(blossom)) {
2.2741 + int bi = (*_node_index)[base];
2.2742 + Value pot = (*_node_data)[bi].pot;
2.2743 +
2.2744 + _matching->set(base, matching);
2.2745 + _blossom_node_list.push_back(base);
2.2746 + _node_potential->set(base, pot);
2.2747 + } else {
2.2748 +
2.2749 + Value pot = (*_blossom_data)[blossom].pot;
2.2750 + int bn = _blossom_node_list.size();
2.2751 +
2.2752 + std::vector<int> subblossoms;
2.2753 + _blossom_set->split(blossom, std::back_inserter(subblossoms));
2.2754 + int b = _blossom_set->find(base);
2.2755 + int ib = -1;
2.2756 + for (int i = 0; i < int(subblossoms.size()); ++i) {
2.2757 + if (subblossoms[i] == b) { ib = i; break; }
2.2758 + }
2.2759 +
2.2760 + for (int i = 1; i < int(subblossoms.size()); i += 2) {
2.2761 + int sb = subblossoms[(ib + i) % subblossoms.size()];
2.2762 + int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
2.2763 +
2.2764 + Arc m = (*_blossom_data)[tb].next;
2.2765 + extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
2.2766 + extractBlossom(tb, _graph.source(m), m);
2.2767 + }
2.2768 + extractBlossom(subblossoms[ib], base, matching);
2.2769 +
2.2770 + int en = _blossom_node_list.size();
2.2771 +
2.2772 + _blossom_potential.push_back(BlossomVariable(bn, en, pot));
2.2773 + }
2.2774 + }
2.2775 +
2.2776 + void extractMatching() {
2.2777 + std::vector<int> blossoms;
2.2778 + for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
2.2779 + blossoms.push_back(c);
2.2780 + }
2.2781 +
2.2782 + for (int i = 0; i < int(blossoms.size()); ++i) {
2.2783 +
2.2784 + Value offset = (*_blossom_data)[blossoms[i]].offset;
2.2785 + (*_blossom_data)[blossoms[i]].pot += 2 * offset;
2.2786 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
2.2787 + n != INVALID; ++n) {
2.2788 + (*_node_data)[(*_node_index)[n]].pot -= offset;
2.2789 + }
2.2790 +
2.2791 + Arc matching = (*_blossom_data)[blossoms[i]].next;
2.2792 + Node base = _graph.source(matching);
2.2793 + extractBlossom(blossoms[i], base, matching);
2.2794 + }
2.2795 + }
2.2796 +
2.2797 + public:
2.2798 +
2.2799 + /// \brief Constructor
2.2800 + ///
2.2801 + /// Constructor.
2.2802 + MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
2.2803 + : _graph(graph), _weight(weight), _matching(0),
2.2804 + _node_potential(0), _blossom_potential(), _blossom_node_list(),
2.2805 + _node_num(0), _blossom_num(0),
2.2806 +
2.2807 + _blossom_index(0), _blossom_set(0), _blossom_data(0),
2.2808 + _node_index(0), _node_heap_index(0), _node_data(0),
2.2809 + _tree_set_index(0), _tree_set(0),
2.2810 +
2.2811 + _delta2_index(0), _delta2(0),
2.2812 + _delta3_index(0), _delta3(0),
2.2813 + _delta4_index(0), _delta4(0),
2.2814 +
2.2815 + _delta_sum() {}
2.2816 +
2.2817 + ~MaxWeightedPerfectMatching() {
2.2818 + destroyStructures();
2.2819 + }
2.2820 +
2.2821 + /// \name Execution control
2.2822 + /// The simplest way to execute the algorithm is to use the
2.2823 + /// \c run() member function.
2.2824 +
2.2825 + ///@{
2.2826 +
2.2827 + /// \brief Initialize the algorithm
2.2828 + ///
2.2829 + /// Initialize the algorithm
2.2830 + void init() {
2.2831 + createStructures();
2.2832 +
2.2833 + for (ArcIt e(_graph); e != INVALID; ++e) {
2.2834 + _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
2.2835 + }
2.2836 + for (EdgeIt e(_graph); e != INVALID; ++e) {
2.2837 + _delta3_index->set(e, _delta3->PRE_HEAP);
2.2838 + }
2.2839 + for (int i = 0; i < _blossom_num; ++i) {
2.2840 + _delta2_index->set(i, _delta2->PRE_HEAP);
2.2841 + _delta4_index->set(i, _delta4->PRE_HEAP);
2.2842 + }
2.2843 +
2.2844 + int index = 0;
2.2845 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.2846 + Value max = - std::numeric_limits<Value>::max();
2.2847 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2.2848 + if (_graph.target(e) == n) continue;
2.2849 + if ((dualScale * _weight[e]) / 2 > max) {
2.2850 + max = (dualScale * _weight[e]) / 2;
2.2851 + }
2.2852 + }
2.2853 + _node_index->set(n, index);
2.2854 + (*_node_data)[index].pot = max;
2.2855 + int blossom =
2.2856 + _blossom_set->insert(n, std::numeric_limits<Value>::max());
2.2857 +
2.2858 + _tree_set->insert(blossom);
2.2859 +
2.2860 + (*_blossom_data)[blossom].status = EVEN;
2.2861 + (*_blossom_data)[blossom].pred = INVALID;
2.2862 + (*_blossom_data)[blossom].next = INVALID;
2.2863 + (*_blossom_data)[blossom].pot = 0;
2.2864 + (*_blossom_data)[blossom].offset = 0;
2.2865 + ++index;
2.2866 + }
2.2867 + for (EdgeIt e(_graph); e != INVALID; ++e) {
2.2868 + int si = (*_node_index)[_graph.u(e)];
2.2869 + int ti = (*_node_index)[_graph.v(e)];
2.2870 + if (_graph.u(e) != _graph.v(e)) {
2.2871 + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
2.2872 + dualScale * _weight[e]) / 2);
2.2873 + }
2.2874 + }
2.2875 + }
2.2876 +
2.2877 + /// \brief Starts the algorithm
2.2878 + ///
2.2879 + /// Starts the algorithm
2.2880 + bool start() {
2.2881 + enum OpType {
2.2882 + D2, D3, D4
2.2883 + };
2.2884 +
2.2885 + int unmatched = _node_num;
2.2886 + while (unmatched > 0) {
2.2887 + Value d2 = !_delta2->empty() ?
2.2888 + _delta2->prio() : std::numeric_limits<Value>::max();
2.2889 +
2.2890 + Value d3 = !_delta3->empty() ?
2.2891 + _delta3->prio() : std::numeric_limits<Value>::max();
2.2892 +
2.2893 + Value d4 = !_delta4->empty() ?
2.2894 + _delta4->prio() : std::numeric_limits<Value>::max();
2.2895 +
2.2896 + _delta_sum = d2; OpType ot = D2;
2.2897 + if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
2.2898 + if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
2.2899 +
2.2900 + if (_delta_sum == std::numeric_limits<Value>::max()) {
2.2901 + return false;
2.2902 + }
2.2903 +
2.2904 + switch (ot) {
2.2905 + case D2:
2.2906 + {
2.2907 + int blossom = _delta2->top();
2.2908 + Node n = _blossom_set->classTop(blossom);
2.2909 + Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
2.2910 + extendOnArc(e);
2.2911 + }
2.2912 + break;
2.2913 + case D3:
2.2914 + {
2.2915 + Edge e = _delta3->top();
2.2916 +
2.2917 + int left_blossom = _blossom_set->find(_graph.u(e));
2.2918 + int right_blossom = _blossom_set->find(_graph.v(e));
2.2919 +
2.2920 + if (left_blossom == right_blossom) {
2.2921 + _delta3->pop();
2.2922 + } else {
2.2923 + int left_tree = _tree_set->find(left_blossom);
2.2924 + int right_tree = _tree_set->find(right_blossom);
2.2925 +
2.2926 + if (left_tree == right_tree) {
2.2927 + shrinkOnEdge(e, left_tree);
2.2928 + } else {
2.2929 + augmentOnEdge(e);
2.2930 + unmatched -= 2;
2.2931 + }
2.2932 + }
2.2933 + } break;
2.2934 + case D4:
2.2935 + splitBlossom(_delta4->top());
2.2936 + break;
2.2937 + }
2.2938 + }
2.2939 + extractMatching();
2.2940 + return true;
2.2941 + }
2.2942 +
2.2943 + /// \brief Runs %MaxWeightedPerfectMatching algorithm.
2.2944 + ///
2.2945 + /// This method runs the %MaxWeightedPerfectMatching algorithm.
2.2946 + ///
2.2947 + /// \note mwm.run() is just a shortcut of the following code.
2.2948 + /// \code
2.2949 + /// mwm.init();
2.2950 + /// mwm.start();
2.2951 + /// \endcode
2.2952 + bool run() {
2.2953 + init();
2.2954 + return start();
2.2955 + }
2.2956 +
2.2957 + /// @}
2.2958 +
2.2959 + /// \name Primal solution
2.2960 + /// Functions to get the primal solution, ie. the matching.
2.2961 +
2.2962 + /// @{
2.2963 +
2.2964 + /// \brief Returns the matching value.
2.2965 + ///
2.2966 + /// Returns the matching value.
2.2967 + Value matchingValue() const {
2.2968 + Value sum = 0;
2.2969 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.2970 + if ((*_matching)[n] != INVALID) {
2.2971 + sum += _weight[(*_matching)[n]];
2.2972 + }
2.2973 + }
2.2974 + return sum /= 2;
2.2975 + }
2.2976 +
2.2977 + /// \brief Returns true when the edge is in the matching.
2.2978 + ///
2.2979 + /// Returns true when the edge is in the matching.
2.2980 + bool matching(const Edge& edge) const {
2.2981 + return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
2.2982 + }
2.2983 +
2.2984 + /// \brief Returns the incident matching edge.
2.2985 + ///
2.2986 + /// Returns the incident matching arc from given edge.
2.2987 + Arc matching(const Node& node) const {
2.2988 + return (*_matching)[node];
2.2989 + }
2.2990 +
2.2991 + /// \brief Returns the mate of the node.
2.2992 + ///
2.2993 + /// Returns the adjancent node in a mathcing arc.
2.2994 + Node mate(const Node& node) const {
2.2995 + return _graph.target((*_matching)[node]);
2.2996 + }
2.2997 +
2.2998 + /// @}
2.2999 +
2.3000 + /// \name Dual solution
2.3001 + /// Functions to get the dual solution.
2.3002 +
2.3003 + /// @{
2.3004 +
2.3005 + /// \brief Returns the value of the dual solution.
2.3006 + ///
2.3007 + /// Returns the value of the dual solution. It should be equal to
2.3008 + /// the primal value scaled by \ref dualScale "dual scale".
2.3009 + Value dualValue() const {
2.3010 + Value sum = 0;
2.3011 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.3012 + sum += nodeValue(n);
2.3013 + }
2.3014 + for (int i = 0; i < blossomNum(); ++i) {
2.3015 + sum += blossomValue(i) * (blossomSize(i) / 2);
2.3016 + }
2.3017 + return sum;
2.3018 + }
2.3019 +
2.3020 + /// \brief Returns the value of the node.
2.3021 + ///
2.3022 + /// Returns the the value of the node.
2.3023 + Value nodeValue(const Node& n) const {
2.3024 + return (*_node_potential)[n];
2.3025 + }
2.3026 +
2.3027 + /// \brief Returns the number of the blossoms in the basis.
2.3028 + ///
2.3029 + /// Returns the number of the blossoms in the basis.
2.3030 + /// \see BlossomIt
2.3031 + int blossomNum() const {
2.3032 + return _blossom_potential.size();
2.3033 + }
2.3034 +
2.3035 +
2.3036 + /// \brief Returns the number of the nodes in the blossom.
2.3037 + ///
2.3038 + /// Returns the number of the nodes in the blossom.
2.3039 + int blossomSize(int k) const {
2.3040 + return _blossom_potential[k].end - _blossom_potential[k].begin;
2.3041 + }
2.3042 +
2.3043 + /// \brief Returns the value of the blossom.
2.3044 + ///
2.3045 + /// Returns the the value of the blossom.
2.3046 + /// \see BlossomIt
2.3047 + Value blossomValue(int k) const {
2.3048 + return _blossom_potential[k].value;
2.3049 + }
2.3050 +
2.3051 + /// \brief Iterator for obtaining the nodes of the blossom.
2.3052 + ///
2.3053 + /// Iterator for obtaining the nodes of the blossom. This class
2.3054 + /// provides a common lemon style iterator for listing a
2.3055 + /// subset of the nodes.
2.3056 + class BlossomIt {
2.3057 + public:
2.3058 +
2.3059 + /// \brief Constructor.
2.3060 + ///
2.3061 + /// Constructor to get the nodes of the variable.
2.3062 + BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
2.3063 + : _algorithm(&algorithm)
2.3064 + {
2.3065 + _index = _algorithm->_blossom_potential[variable].begin;
2.3066 + _last = _algorithm->_blossom_potential[variable].end;
2.3067 + }
2.3068 +
2.3069 + /// \brief Conversion to node.
2.3070 + ///
2.3071 + /// Conversion to node.
2.3072 + operator Node() const {
2.3073 + return _algorithm->_blossom_node_list[_index];
2.3074 + }
2.3075 +
2.3076 + /// \brief Increment operator.
2.3077 + ///
2.3078 + /// Increment operator.
2.3079 + BlossomIt& operator++() {
2.3080 + ++_index;
2.3081 + return *this;
2.3082 + }
2.3083 +
2.3084 + /// \brief Validity checking
2.3085 + ///
2.3086 + /// Checks whether the iterator is invalid.
2.3087 + bool operator==(Invalid) const { return _index == _last; }
2.3088 +
2.3089 + /// \brief Validity checking
2.3090 + ///
2.3091 + /// Checks whether the iterator is valid.
2.3092 + bool operator!=(Invalid) const { return _index != _last; }
2.3093 +
2.3094 + private:
2.3095 + const MaxWeightedPerfectMatching* _algorithm;
2.3096 + int _last;
2.3097 + int _index;
2.3098 + };
2.3099 +
2.3100 + /// @}
2.3101 +
2.3102 + };
2.3103 +
2.3104 +
2.3105 +} //END OF NAMESPACE LEMON
2.3106 +
2.3107 +#endif //LEMON_MAX_MATCHING_H
3.1 --- a/test/CMakeLists.txt Wed Oct 22 14:39:04 2008 +0100
3.2 +++ b/test/CMakeLists.txt Wed Oct 22 14:41:18 2008 +0100
3.3 @@ -16,6 +16,7 @@
3.4 heap_test
3.5 kruskal_test
3.6 maps_test
3.7 + max_matching_test
3.8 random_test
3.9 path_test
3.10 time_measure_test
4.1 --- a/test/Makefile.am Wed Oct 22 14:39:04 2008 +0100
4.2 +++ b/test/Makefile.am Wed Oct 22 14:41:18 2008 +0100
4.3 @@ -19,6 +19,7 @@
4.4 test/heap_test \
4.5 test/kruskal_test \
4.6 test/maps_test \
4.7 + test/max_matching_test \
4.8 test/random_test \
4.9 test/path_test \
4.10 test/test_tools_fail \
4.11 @@ -42,6 +43,7 @@
4.12 test_heap_test_SOURCES = test/heap_test.cc
4.13 test_kruskal_test_SOURCES = test/kruskal_test.cc
4.14 test_maps_test_SOURCES = test/maps_test.cc
4.15 +test_max_matching_test_SOURCES = test/max_matching_test.cc
4.16 test_path_test_SOURCES = test/path_test.cc
4.17 test_random_test_SOURCES = test/random_test.cc
4.18 test_test_tools_fail_SOURCES = test/test_tools_fail.cc
5.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
5.2 +++ b/test/max_matching_test.cc Wed Oct 22 14:41:18 2008 +0100
5.3 @@ -0,0 +1,310 @@
5.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
5.5 + *
5.6 + * This file is a part of LEMON, a generic C++ optimization library.
5.7 + *
5.8 + * Copyright (C) 2003-2008
5.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
5.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
5.11 + *
5.12 + * Permission to use, modify and distribute this software is granted
5.13 + * provided that this copyright notice appears in all copies. For
5.14 + * precise terms see the accompanying LICENSE file.
5.15 + *
5.16 + * This software is provided "AS IS" with no warranty of any kind,
5.17 + * express or implied, and with no claim as to its suitability for any
5.18 + * purpose.
5.19 + *
5.20 + */
5.21 +
5.22 +#include <iostream>
5.23 +#include <sstream>
5.24 +#include <vector>
5.25 +#include <queue>
5.26 +#include <lemon/math.h>
5.27 +#include <cstdlib>
5.28 +
5.29 +#include <lemon/max_matching.h>
5.30 +#include <lemon/smart_graph.h>
5.31 +#include <lemon/lgf_reader.h>
5.32 +
5.33 +#include "test_tools.h"
5.34 +
5.35 +using namespace std;
5.36 +using namespace lemon;
5.37 +
5.38 +GRAPH_TYPEDEFS(SmartGraph);
5.39 +
5.40 +
5.41 +const int lgfn = 3;
5.42 +const std::string lgf[lgfn] = {
5.43 + "@nodes\n"
5.44 + "label\n"
5.45 + "0\n"
5.46 + "1\n"
5.47 + "2\n"
5.48 + "3\n"
5.49 + "4\n"
5.50 + "5\n"
5.51 + "6\n"
5.52 + "7\n"
5.53 + "@edges\n"
5.54 + " label weight\n"
5.55 + "7 4 0 984\n"
5.56 + "0 7 1 73\n"
5.57 + "7 1 2 204\n"
5.58 + "2 3 3 583\n"
5.59 + "2 7 4 565\n"
5.60 + "2 1 5 582\n"
5.61 + "0 4 6 551\n"
5.62 + "2 5 7 385\n"
5.63 + "1 5 8 561\n"
5.64 + "5 3 9 484\n"
5.65 + "7 5 10 904\n"
5.66 + "3 6 11 47\n"
5.67 + "7 6 12 888\n"
5.68 + "3 0 13 747\n"
5.69 + "6 1 14 310\n",
5.70 +
5.71 + "@nodes\n"
5.72 + "label\n"
5.73 + "0\n"
5.74 + "1\n"
5.75 + "2\n"
5.76 + "3\n"
5.77 + "4\n"
5.78 + "5\n"
5.79 + "6\n"
5.80 + "7\n"
5.81 + "@edges\n"
5.82 + " label weight\n"
5.83 + "2 5 0 710\n"
5.84 + "0 5 1 241\n"
5.85 + "2 4 2 856\n"
5.86 + "2 6 3 762\n"
5.87 + "4 1 4 747\n"
5.88 + "6 1 5 962\n"
5.89 + "4 7 6 723\n"
5.90 + "1 7 7 661\n"
5.91 + "2 3 8 376\n"
5.92 + "1 0 9 416\n"
5.93 + "6 7 10 391\n",
5.94 +
5.95 + "@nodes\n"
5.96 + "label\n"
5.97 + "0\n"
5.98 + "1\n"
5.99 + "2\n"
5.100 + "3\n"
5.101 + "4\n"
5.102 + "5\n"
5.103 + "6\n"
5.104 + "7\n"
5.105 + "@edges\n"
5.106 + " label weight\n"
5.107 + "6 2 0 553\n"
5.108 + "0 7 1 653\n"
5.109 + "6 3 2 22\n"
5.110 + "4 7 3 846\n"
5.111 + "7 2 4 981\n"
5.112 + "7 6 5 250\n"
5.113 + "5 2 6 539\n",
5.114 +};
5.115 +
5.116 +void checkMatching(const SmartGraph& graph,
5.117 + const MaxMatching<SmartGraph>& mm) {
5.118 + int num = 0;
5.119 +
5.120 + IntNodeMap comp_index(graph);
5.121 + UnionFind<IntNodeMap> comp(comp_index);
5.122 +
5.123 + int barrier_num = 0;
5.124 +
5.125 + for (NodeIt n(graph); n != INVALID; ++n) {
5.126 + check(mm.decomposition(n) == MaxMatching<SmartGraph>::EVEN ||
5.127 + mm.matching(n) != INVALID, "Wrong Gallai-Edmonds decomposition");
5.128 + if (mm.decomposition(n) == MaxMatching<SmartGraph>::ODD) {
5.129 + ++barrier_num;
5.130 + } else {
5.131 + comp.insert(n);
5.132 + }
5.133 + }
5.134 +
5.135 + for (EdgeIt e(graph); e != INVALID; ++e) {
5.136 + if (mm.matching(e)) {
5.137 + check(e == mm.matching(graph.u(e)), "Wrong matching");
5.138 + check(e == mm.matching(graph.v(e)), "Wrong matching");
5.139 + ++num;
5.140 + }
5.141 + check(mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::EVEN ||
5.142 + mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::MATCHED,
5.143 + "Wrong Gallai-Edmonds decomposition");
5.144 +
5.145 + check(mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::EVEN ||
5.146 + mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::MATCHED,
5.147 + "Wrong Gallai-Edmonds decomposition");
5.148 +
5.149 + if (mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::ODD &&
5.150 + mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::ODD) {
5.151 + comp.join(graph.u(e), graph.v(e));
5.152 + }
5.153 + }
5.154 +
5.155 + std::set<int> comp_root;
5.156 + int odd_comp_num = 0;
5.157 + for (NodeIt n(graph); n != INVALID; ++n) {
5.158 + if (mm.decomposition(n) != MaxMatching<SmartGraph>::ODD) {
5.159 + int root = comp.find(n);
5.160 + if (comp_root.find(root) == comp_root.end()) {
5.161 + comp_root.insert(root);
5.162 + if (comp.size(n) % 2 == 1) {
5.163 + ++odd_comp_num;
5.164 + }
5.165 + }
5.166 + }
5.167 + }
5.168 +
5.169 + check(mm.matchingSize() == num, "Wrong matching");
5.170 + check(2 * num == countNodes(graph) - (odd_comp_num - barrier_num),
5.171 + "Wrong matching");
5.172 + return;
5.173 +}
5.174 +
5.175 +void checkWeightedMatching(const SmartGraph& graph,
5.176 + const SmartGraph::EdgeMap<int>& weight,
5.177 + const MaxWeightedMatching<SmartGraph>& mwm) {
5.178 + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
5.179 + if (graph.u(e) == graph.v(e)) continue;
5.180 + int rw = mwm.nodeValue(graph.u(e)) + mwm.nodeValue(graph.v(e));
5.181 +
5.182 + for (int i = 0; i < mwm.blossomNum(); ++i) {
5.183 + bool s = false, t = false;
5.184 + for (MaxWeightedMatching<SmartGraph>::BlossomIt n(mwm, i);
5.185 + n != INVALID; ++n) {
5.186 + if (graph.u(e) == n) s = true;
5.187 + if (graph.v(e) == n) t = true;
5.188 + }
5.189 + if (s == true && t == true) {
5.190 + rw += mwm.blossomValue(i);
5.191 + }
5.192 + }
5.193 + rw -= weight[e] * mwm.dualScale;
5.194 +
5.195 + check(rw >= 0, "Negative reduced weight");
5.196 + check(rw == 0 || !mwm.matching(e),
5.197 + "Non-zero reduced weight on matching edge");
5.198 + }
5.199 +
5.200 + int pv = 0;
5.201 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
5.202 + if (mwm.matching(n) != INVALID) {
5.203 + check(mwm.nodeValue(n) >= 0, "Invalid node value");
5.204 + pv += weight[mwm.matching(n)];
5.205 + SmartGraph::Node o = graph.target(mwm.matching(n));
5.206 + check(mwm.mate(n) == o, "Invalid matching");
5.207 + check(mwm.matching(n) == graph.oppositeArc(mwm.matching(o)),
5.208 + "Invalid matching");
5.209 + } else {
5.210 + check(mwm.mate(n) == INVALID, "Invalid matching");
5.211 + check(mwm.nodeValue(n) == 0, "Invalid matching");
5.212 + }
5.213 + }
5.214 +
5.215 + int dv = 0;
5.216 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
5.217 + dv += mwm.nodeValue(n);
5.218 + }
5.219 +
5.220 + for (int i = 0; i < mwm.blossomNum(); ++i) {
5.221 + check(mwm.blossomValue(i) >= 0, "Invalid blossom value");
5.222 + check(mwm.blossomSize(i) % 2 == 1, "Even blossom size");
5.223 + dv += mwm.blossomValue(i) * ((mwm.blossomSize(i) - 1) / 2);
5.224 + }
5.225 +
5.226 + check(pv * mwm.dualScale == dv * 2, "Wrong duality");
5.227 +
5.228 + return;
5.229 +}
5.230 +
5.231 +void checkWeightedPerfectMatching(const SmartGraph& graph,
5.232 + const SmartGraph::EdgeMap<int>& weight,
5.233 + const MaxWeightedPerfectMatching<SmartGraph>& mwpm) {
5.234 + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
5.235 + if (graph.u(e) == graph.v(e)) continue;
5.236 + int rw = mwpm.nodeValue(graph.u(e)) + mwpm.nodeValue(graph.v(e));
5.237 +
5.238 + for (int i = 0; i < mwpm.blossomNum(); ++i) {
5.239 + bool s = false, t = false;
5.240 + for (MaxWeightedPerfectMatching<SmartGraph>::BlossomIt n(mwpm, i);
5.241 + n != INVALID; ++n) {
5.242 + if (graph.u(e) == n) s = true;
5.243 + if (graph.v(e) == n) t = true;
5.244 + }
5.245 + if (s == true && t == true) {
5.246 + rw += mwpm.blossomValue(i);
5.247 + }
5.248 + }
5.249 + rw -= weight[e] * mwpm.dualScale;
5.250 +
5.251 + check(rw >= 0, "Negative reduced weight");
5.252 + check(rw == 0 || !mwpm.matching(e),
5.253 + "Non-zero reduced weight on matching edge");
5.254 + }
5.255 +
5.256 + int pv = 0;
5.257 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
5.258 + check(mwpm.matching(n) != INVALID, "Non perfect");
5.259 + pv += weight[mwpm.matching(n)];
5.260 + SmartGraph::Node o = graph.target(mwpm.matching(n));
5.261 + check(mwpm.mate(n) == o, "Invalid matching");
5.262 + check(mwpm.matching(n) == graph.oppositeArc(mwpm.matching(o)),
5.263 + "Invalid matching");
5.264 + }
5.265 +
5.266 + int dv = 0;
5.267 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
5.268 + dv += mwpm.nodeValue(n);
5.269 + }
5.270 +
5.271 + for (int i = 0; i < mwpm.blossomNum(); ++i) {
5.272 + check(mwpm.blossomValue(i) >= 0, "Invalid blossom value");
5.273 + check(mwpm.blossomSize(i) % 2 == 1, "Even blossom size");
5.274 + dv += mwpm.blossomValue(i) * ((mwpm.blossomSize(i) - 1) / 2);
5.275 + }
5.276 +
5.277 + check(pv * mwpm.dualScale == dv * 2, "Wrong duality");
5.278 +
5.279 + return;
5.280 +}
5.281 +
5.282 +
5.283 +int main() {
5.284 +
5.285 + for (int i = 0; i < lgfn; ++i) {
5.286 + SmartGraph graph;
5.287 + SmartGraph::EdgeMap<int> weight(graph);
5.288 +
5.289 + istringstream lgfs(lgf[i]);
5.290 + graphReader(graph, lgfs).
5.291 + edgeMap("weight", weight).run();
5.292 +
5.293 + MaxMatching<SmartGraph> mm(graph);
5.294 + mm.run();
5.295 + checkMatching(graph, mm);
5.296 +
5.297 + MaxWeightedMatching<SmartGraph> mwm(graph, weight);
5.298 + mwm.run();
5.299 + checkWeightedMatching(graph, weight, mwm);
5.300 +
5.301 + MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
5.302 + bool perfect = mwpm.run();
5.303 +
5.304 + check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
5.305 + "Perfect matching found");
5.306 +
5.307 + if (perfect) {
5.308 + checkWeightedPerfectMatching(graph, weight, mwpm);
5.309 + }
5.310 + }
5.311 +
5.312 + return 0;
5.313 +}