1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/max_matching.h Wed Oct 22 14:41:18 2008 +0100
1.3 @@ -0,0 +1,3104 @@
1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library.
1.7 + *
1.8 + * Copyright (C) 2003-2008
1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 + *
1.12 + * Permission to use, modify and distribute this software is granted
1.13 + * provided that this copyright notice appears in all copies. For
1.14 + * precise terms see the accompanying LICENSE file.
1.15 + *
1.16 + * This software is provided "AS IS" with no warranty of any kind,
1.17 + * express or implied, and with no claim as to its suitability for any
1.18 + * purpose.
1.19 + *
1.20 + */
1.21 +
1.22 +#ifndef LEMON_MAX_MATCHING_H
1.23 +#define LEMON_MAX_MATCHING_H
1.24 +
1.25 +#include <vector>
1.26 +#include <queue>
1.27 +#include <set>
1.28 +#include <limits>
1.29 +
1.30 +#include <lemon/core.h>
1.31 +#include <lemon/unionfind.h>
1.32 +#include <lemon/bin_heap.h>
1.33 +#include <lemon/maps.h>
1.34 +
1.35 +///\ingroup matching
1.36 +///\file
1.37 +///\brief Maximum matching algorithms in general graphs.
1.38 +
1.39 +namespace lemon {
1.40 +
1.41 + /// \ingroup matching
1.42 + ///
1.43 + /// \brief Edmonds' alternating forest maximum matching algorithm.
1.44 + ///
1.45 + /// This class implements Edmonds' alternating forest matching
1.46 + /// algorithm. The algorithm can be started from an arbitrary initial
1.47 + /// matching (the default is the empty one)
1.48 + ///
1.49 + /// The dual solution of the problem is a map of the nodes to
1.50 + /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c
1.51 + /// MATCHED/C showing the Gallai-Edmonds decomposition of the
1.52 + /// graph. The nodes in \c EVEN/D induce a graph with
1.53 + /// factor-critical components, the nodes in \c ODD/A form the
1.54 + /// barrier, and the nodes in \c MATCHED/C induce a graph having a
1.55 + /// perfect matching. The number of the factor-critical components
1.56 + /// minus the number of barrier nodes is a lower bound on the
1.57 + /// unmatched nodes, and the matching is optimal if and only if this bound is
1.58 + /// tight. This decomposition can be attained by calling \c
1.59 + /// decomposition() after running the algorithm.
1.60 + ///
1.61 + /// \param _Graph The graph type the algorithm runs on.
1.62 + template <typename _Graph>
1.63 + class MaxMatching {
1.64 + public:
1.65 +
1.66 + typedef _Graph Graph;
1.67 + typedef typename Graph::template NodeMap<typename Graph::Arc>
1.68 + MatchingMap;
1.69 +
1.70 + ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
1.71 + ///
1.72 + ///Indicates the Gallai-Edmonds decomposition of the graph. The
1.73 + ///nodes with Status \c EVEN/D induce a graph with factor-critical
1.74 + ///components, the nodes in \c ODD/A form the canonical barrier,
1.75 + ///and the nodes in \c MATCHED/C induce a graph having a perfect
1.76 + ///matching.
1.77 + enum Status {
1.78 + EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2
1.79 + };
1.80 +
1.81 + typedef typename Graph::template NodeMap<Status> StatusMap;
1.82 +
1.83 + private:
1.84 +
1.85 + TEMPLATE_GRAPH_TYPEDEFS(Graph);
1.86 +
1.87 + typedef UnionFindEnum<IntNodeMap> BlossomSet;
1.88 + typedef ExtendFindEnum<IntNodeMap> TreeSet;
1.89 + typedef RangeMap<Node> NodeIntMap;
1.90 + typedef MatchingMap EarMap;
1.91 + typedef std::vector<Node> NodeQueue;
1.92 +
1.93 + const Graph& _graph;
1.94 + MatchingMap* _matching;
1.95 + StatusMap* _status;
1.96 +
1.97 + EarMap* _ear;
1.98 +
1.99 + IntNodeMap* _blossom_set_index;
1.100 + BlossomSet* _blossom_set;
1.101 + NodeIntMap* _blossom_rep;
1.102 +
1.103 + IntNodeMap* _tree_set_index;
1.104 + TreeSet* _tree_set;
1.105 +
1.106 + NodeQueue _node_queue;
1.107 + int _process, _postpone, _last;
1.108 +
1.109 + int _node_num;
1.110 +
1.111 + private:
1.112 +
1.113 + void createStructures() {
1.114 + _node_num = countNodes(_graph);
1.115 + if (!_matching) {
1.116 + _matching = new MatchingMap(_graph);
1.117 + }
1.118 + if (!_status) {
1.119 + _status = new StatusMap(_graph);
1.120 + }
1.121 + if (!_ear) {
1.122 + _ear = new EarMap(_graph);
1.123 + }
1.124 + if (!_blossom_set) {
1.125 + _blossom_set_index = new IntNodeMap(_graph);
1.126 + _blossom_set = new BlossomSet(*_blossom_set_index);
1.127 + }
1.128 + if (!_blossom_rep) {
1.129 + _blossom_rep = new NodeIntMap(_node_num);
1.130 + }
1.131 + if (!_tree_set) {
1.132 + _tree_set_index = new IntNodeMap(_graph);
1.133 + _tree_set = new TreeSet(*_tree_set_index);
1.134 + }
1.135 + _node_queue.resize(_node_num);
1.136 + }
1.137 +
1.138 + void destroyStructures() {
1.139 + if (_matching) {
1.140 + delete _matching;
1.141 + }
1.142 + if (_status) {
1.143 + delete _status;
1.144 + }
1.145 + if (_ear) {
1.146 + delete _ear;
1.147 + }
1.148 + if (_blossom_set) {
1.149 + delete _blossom_set;
1.150 + delete _blossom_set_index;
1.151 + }
1.152 + if (_blossom_rep) {
1.153 + delete _blossom_rep;
1.154 + }
1.155 + if (_tree_set) {
1.156 + delete _tree_set_index;
1.157 + delete _tree_set;
1.158 + }
1.159 + }
1.160 +
1.161 + void processDense(const Node& n) {
1.162 + _process = _postpone = _last = 0;
1.163 + _node_queue[_last++] = n;
1.164 +
1.165 + while (_process != _last) {
1.166 + Node u = _node_queue[_process++];
1.167 + for (OutArcIt a(_graph, u); a != INVALID; ++a) {
1.168 + Node v = _graph.target(a);
1.169 + if ((*_status)[v] == MATCHED) {
1.170 + extendOnArc(a);
1.171 + } else if ((*_status)[v] == UNMATCHED) {
1.172 + augmentOnArc(a);
1.173 + return;
1.174 + }
1.175 + }
1.176 + }
1.177 +
1.178 + while (_postpone != _last) {
1.179 + Node u = _node_queue[_postpone++];
1.180 +
1.181 + for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
1.182 + Node v = _graph.target(a);
1.183 +
1.184 + if ((*_status)[v] == EVEN) {
1.185 + if (_blossom_set->find(u) != _blossom_set->find(v)) {
1.186 + shrinkOnEdge(a);
1.187 + }
1.188 + }
1.189 +
1.190 + while (_process != _last) {
1.191 + Node w = _node_queue[_process++];
1.192 + for (OutArcIt b(_graph, w); b != INVALID; ++b) {
1.193 + Node x = _graph.target(b);
1.194 + if ((*_status)[x] == MATCHED) {
1.195 + extendOnArc(b);
1.196 + } else if ((*_status)[x] == UNMATCHED) {
1.197 + augmentOnArc(b);
1.198 + return;
1.199 + }
1.200 + }
1.201 + }
1.202 + }
1.203 + }
1.204 + }
1.205 +
1.206 + void processSparse(const Node& n) {
1.207 + _process = _last = 0;
1.208 + _node_queue[_last++] = n;
1.209 + while (_process != _last) {
1.210 + Node u = _node_queue[_process++];
1.211 + for (OutArcIt a(_graph, u); a != INVALID; ++a) {
1.212 + Node v = _graph.target(a);
1.213 +
1.214 + if ((*_status)[v] == EVEN) {
1.215 + if (_blossom_set->find(u) != _blossom_set->find(v)) {
1.216 + shrinkOnEdge(a);
1.217 + }
1.218 + } else if ((*_status)[v] == MATCHED) {
1.219 + extendOnArc(a);
1.220 + } else if ((*_status)[v] == UNMATCHED) {
1.221 + augmentOnArc(a);
1.222 + return;
1.223 + }
1.224 + }
1.225 + }
1.226 + }
1.227 +
1.228 + void shrinkOnEdge(const Edge& e) {
1.229 + Node nca = INVALID;
1.230 +
1.231 + {
1.232 + std::set<Node> left_set, right_set;
1.233 +
1.234 + Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
1.235 + left_set.insert(left);
1.236 +
1.237 + Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
1.238 + right_set.insert(right);
1.239 +
1.240 + while (true) {
1.241 + if ((*_matching)[left] == INVALID) break;
1.242 + left = _graph.target((*_matching)[left]);
1.243 + left = (*_blossom_rep)[_blossom_set->
1.244 + find(_graph.target((*_ear)[left]))];
1.245 + if (right_set.find(left) != right_set.end()) {
1.246 + nca = left;
1.247 + break;
1.248 + }
1.249 + left_set.insert(left);
1.250 +
1.251 + if ((*_matching)[right] == INVALID) break;
1.252 + right = _graph.target((*_matching)[right]);
1.253 + right = (*_blossom_rep)[_blossom_set->
1.254 + find(_graph.target((*_ear)[right]))];
1.255 + if (left_set.find(right) != left_set.end()) {
1.256 + nca = right;
1.257 + break;
1.258 + }
1.259 + right_set.insert(right);
1.260 + }
1.261 +
1.262 + if (nca == INVALID) {
1.263 + if ((*_matching)[left] == INVALID) {
1.264 + nca = right;
1.265 + while (left_set.find(nca) == left_set.end()) {
1.266 + nca = _graph.target((*_matching)[nca]);
1.267 + nca =(*_blossom_rep)[_blossom_set->
1.268 + find(_graph.target((*_ear)[nca]))];
1.269 + }
1.270 + } else {
1.271 + nca = left;
1.272 + while (right_set.find(nca) == right_set.end()) {
1.273 + nca = _graph.target((*_matching)[nca]);
1.274 + nca = (*_blossom_rep)[_blossom_set->
1.275 + find(_graph.target((*_ear)[nca]))];
1.276 + }
1.277 + }
1.278 + }
1.279 + }
1.280 +
1.281 + {
1.282 +
1.283 + Node node = _graph.u(e);
1.284 + Arc arc = _graph.direct(e, true);
1.285 + Node base = (*_blossom_rep)[_blossom_set->find(node)];
1.286 +
1.287 + while (base != nca) {
1.288 + _ear->set(node, arc);
1.289 +
1.290 + Node n = node;
1.291 + while (n != base) {
1.292 + n = _graph.target((*_matching)[n]);
1.293 + Arc a = (*_ear)[n];
1.294 + n = _graph.target(a);
1.295 + _ear->set(n, _graph.oppositeArc(a));
1.296 + }
1.297 + node = _graph.target((*_matching)[base]);
1.298 + _tree_set->erase(base);
1.299 + _tree_set->erase(node);
1.300 + _blossom_set->insert(node, _blossom_set->find(base));
1.301 + _status->set(node, EVEN);
1.302 + _node_queue[_last++] = node;
1.303 + arc = _graph.oppositeArc((*_ear)[node]);
1.304 + node = _graph.target((*_ear)[node]);
1.305 + base = (*_blossom_rep)[_blossom_set->find(node)];
1.306 + _blossom_set->join(_graph.target(arc), base);
1.307 + }
1.308 + }
1.309 +
1.310 + _blossom_rep->set(_blossom_set->find(nca), nca);
1.311 +
1.312 + {
1.313 +
1.314 + Node node = _graph.v(e);
1.315 + Arc arc = _graph.direct(e, false);
1.316 + Node base = (*_blossom_rep)[_blossom_set->find(node)];
1.317 +
1.318 + while (base != nca) {
1.319 + _ear->set(node, arc);
1.320 +
1.321 + Node n = node;
1.322 + while (n != base) {
1.323 + n = _graph.target((*_matching)[n]);
1.324 + Arc a = (*_ear)[n];
1.325 + n = _graph.target(a);
1.326 + _ear->set(n, _graph.oppositeArc(a));
1.327 + }
1.328 + node = _graph.target((*_matching)[base]);
1.329 + _tree_set->erase(base);
1.330 + _tree_set->erase(node);
1.331 + _blossom_set->insert(node, _blossom_set->find(base));
1.332 + _status->set(node, EVEN);
1.333 + _node_queue[_last++] = node;
1.334 + arc = _graph.oppositeArc((*_ear)[node]);
1.335 + node = _graph.target((*_ear)[node]);
1.336 + base = (*_blossom_rep)[_blossom_set->find(node)];
1.337 + _blossom_set->join(_graph.target(arc), base);
1.338 + }
1.339 + }
1.340 +
1.341 + _blossom_rep->set(_blossom_set->find(nca), nca);
1.342 + }
1.343 +
1.344 +
1.345 +
1.346 + void extendOnArc(const Arc& a) {
1.347 + Node base = _graph.source(a);
1.348 + Node odd = _graph.target(a);
1.349 +
1.350 + _ear->set(odd, _graph.oppositeArc(a));
1.351 + Node even = _graph.target((*_matching)[odd]);
1.352 + _blossom_rep->set(_blossom_set->insert(even), even);
1.353 + _status->set(odd, ODD);
1.354 + _status->set(even, EVEN);
1.355 + int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
1.356 + _tree_set->insert(odd, tree);
1.357 + _tree_set->insert(even, tree);
1.358 + _node_queue[_last++] = even;
1.359 +
1.360 + }
1.361 +
1.362 + void augmentOnArc(const Arc& a) {
1.363 + Node even = _graph.source(a);
1.364 + Node odd = _graph.target(a);
1.365 +
1.366 + int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
1.367 +
1.368 + _matching->set(odd, _graph.oppositeArc(a));
1.369 + _status->set(odd, MATCHED);
1.370 +
1.371 + Arc arc = (*_matching)[even];
1.372 + _matching->set(even, a);
1.373 +
1.374 + while (arc != INVALID) {
1.375 + odd = _graph.target(arc);
1.376 + arc = (*_ear)[odd];
1.377 + even = _graph.target(arc);
1.378 + _matching->set(odd, arc);
1.379 + arc = (*_matching)[even];
1.380 + _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
1.381 + }
1.382 +
1.383 + for (typename TreeSet::ItemIt it(*_tree_set, tree);
1.384 + it != INVALID; ++it) {
1.385 + if ((*_status)[it] == ODD) {
1.386 + _status->set(it, MATCHED);
1.387 + } else {
1.388 + int blossom = _blossom_set->find(it);
1.389 + for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
1.390 + jt != INVALID; ++jt) {
1.391 + _status->set(jt, MATCHED);
1.392 + }
1.393 + _blossom_set->eraseClass(blossom);
1.394 + }
1.395 + }
1.396 + _tree_set->eraseClass(tree);
1.397 +
1.398 + }
1.399 +
1.400 + public:
1.401 +
1.402 + /// \brief Constructor
1.403 + ///
1.404 + /// Constructor.
1.405 + MaxMatching(const Graph& graph)
1.406 + : _graph(graph), _matching(0), _status(0), _ear(0),
1.407 + _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
1.408 + _tree_set_index(0), _tree_set(0) {}
1.409 +
1.410 + ~MaxMatching() {
1.411 + destroyStructures();
1.412 + }
1.413 +
1.414 + /// \name Execution control
1.415 + /// The simplest way to execute the algorithm is to use the
1.416 + /// \c run() member function.
1.417 + /// \n
1.418 +
1.419 + /// If you need better control on the execution, you must call
1.420 + /// \ref init(), \ref greedyInit() or \ref matchingInit()
1.421 + /// functions first, then you can start the algorithm with the \ref
1.422 + /// startParse() or startDense() functions.
1.423 +
1.424 + ///@{
1.425 +
1.426 + /// \brief Sets the actual matching to the empty matching.
1.427 + ///
1.428 + /// Sets the actual matching to the empty matching.
1.429 + ///
1.430 + void init() {
1.431 + createStructures();
1.432 + for(NodeIt n(_graph); n != INVALID; ++n) {
1.433 + _matching->set(n, INVALID);
1.434 + _status->set(n, UNMATCHED);
1.435 + }
1.436 + }
1.437 +
1.438 + ///\brief Finds an initial matching in a greedy way
1.439 + ///
1.440 + ///It finds an initial matching in a greedy way.
1.441 + void greedyInit() {
1.442 + createStructures();
1.443 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.444 + _matching->set(n, INVALID);
1.445 + _status->set(n, UNMATCHED);
1.446 + }
1.447 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.448 + if ((*_matching)[n] == INVALID) {
1.449 + for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
1.450 + Node v = _graph.target(a);
1.451 + if ((*_matching)[v] == INVALID && v != n) {
1.452 + _matching->set(n, a);
1.453 + _status->set(n, MATCHED);
1.454 + _matching->set(v, _graph.oppositeArc(a));
1.455 + _status->set(v, MATCHED);
1.456 + break;
1.457 + }
1.458 + }
1.459 + }
1.460 + }
1.461 + }
1.462 +
1.463 +
1.464 + /// \brief Initialize the matching from a map containing.
1.465 + ///
1.466 + /// Initialize the matching from a \c bool valued \c Edge map. This
1.467 + /// map must have the property that there are no two incident edges
1.468 + /// with true value, ie. it contains a matching.
1.469 + /// \return %True if the map contains a matching.
1.470 + template <typename MatchingMap>
1.471 + bool matchingInit(const MatchingMap& matching) {
1.472 + createStructures();
1.473 +
1.474 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.475 + _matching->set(n, INVALID);
1.476 + _status->set(n, UNMATCHED);
1.477 + }
1.478 + for(EdgeIt e(_graph); e!=INVALID; ++e) {
1.479 + if (matching[e]) {
1.480 +
1.481 + Node u = _graph.u(e);
1.482 + if ((*_matching)[u] != INVALID) return false;
1.483 + _matching->set(u, _graph.direct(e, true));
1.484 + _status->set(u, MATCHED);
1.485 +
1.486 + Node v = _graph.v(e);
1.487 + if ((*_matching)[v] != INVALID) return false;
1.488 + _matching->set(v, _graph.direct(e, false));
1.489 + _status->set(v, MATCHED);
1.490 + }
1.491 + }
1.492 + return true;
1.493 + }
1.494 +
1.495 + /// \brief Starts Edmonds' algorithm
1.496 + ///
1.497 + /// If runs the original Edmonds' algorithm.
1.498 + void startSparse() {
1.499 + for(NodeIt n(_graph); n != INVALID; ++n) {
1.500 + if ((*_status)[n] == UNMATCHED) {
1.501 + (*_blossom_rep)[_blossom_set->insert(n)] = n;
1.502 + _tree_set->insert(n);
1.503 + _status->set(n, EVEN);
1.504 + processSparse(n);
1.505 + }
1.506 + }
1.507 + }
1.508 +
1.509 + /// \brief Starts Edmonds' algorithm.
1.510 + ///
1.511 + /// It runs Edmonds' algorithm with a heuristic of postponing
1.512 + /// shrinks, therefore resulting in a faster algorithm for dense graphs.
1.513 + void startDense() {
1.514 + for(NodeIt n(_graph); n != INVALID; ++n) {
1.515 + if ((*_status)[n] == UNMATCHED) {
1.516 + (*_blossom_rep)[_blossom_set->insert(n)] = n;
1.517 + _tree_set->insert(n);
1.518 + _status->set(n, EVEN);
1.519 + processDense(n);
1.520 + }
1.521 + }
1.522 + }
1.523 +
1.524 +
1.525 + /// \brief Runs Edmonds' algorithm
1.526 + ///
1.527 + /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>)
1.528 + /// or Edmonds' algorithm with a heuristic of
1.529 + /// postponing shrinks for dense graphs.
1.530 + void run() {
1.531 + if (countEdges(_graph) < 2 * countNodes(_graph)) {
1.532 + greedyInit();
1.533 + startSparse();
1.534 + } else {
1.535 + init();
1.536 + startDense();
1.537 + }
1.538 + }
1.539 +
1.540 + /// @}
1.541 +
1.542 + /// \name Primal solution
1.543 + /// Functions to get the primal solution, ie. the matching.
1.544 +
1.545 + /// @{
1.546 +
1.547 + ///\brief Returns the size of the current matching.
1.548 + ///
1.549 + ///Returns the size of the current matching. After \ref
1.550 + ///run() it returns the size of the maximum matching in the graph.
1.551 + int matchingSize() const {
1.552 + int size = 0;
1.553 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.554 + if ((*_matching)[n] != INVALID) {
1.555 + ++size;
1.556 + }
1.557 + }
1.558 + return size / 2;
1.559 + }
1.560 +
1.561 + /// \brief Returns true when the edge is in the matching.
1.562 + ///
1.563 + /// Returns true when the edge is in the matching.
1.564 + bool matching(const Edge& edge) const {
1.565 + return edge == (*_matching)[_graph.u(edge)];
1.566 + }
1.567 +
1.568 + /// \brief Returns the matching edge incident to the given node.
1.569 + ///
1.570 + /// Returns the matching edge of a \c node in the actual matching or
1.571 + /// INVALID if the \c node is not covered by the actual matching.
1.572 + Arc matching(const Node& n) const {
1.573 + return (*_matching)[n];
1.574 + }
1.575 +
1.576 + ///\brief Returns the mate of a node in the actual matching.
1.577 + ///
1.578 + ///Returns the mate of a \c node in the actual matching or
1.579 + ///INVALID if the \c node is not covered by the actual matching.
1.580 + Node mate(const Node& n) const {
1.581 + return (*_matching)[n] != INVALID ?
1.582 + _graph.target((*_matching)[n]) : INVALID;
1.583 + }
1.584 +
1.585 + /// @}
1.586 +
1.587 + /// \name Dual solution
1.588 + /// Functions to get the dual solution, ie. the decomposition.
1.589 +
1.590 + /// @{
1.591 +
1.592 + /// \brief Returns the class of the node in the Edmonds-Gallai
1.593 + /// decomposition.
1.594 + ///
1.595 + /// Returns the class of the node in the Edmonds-Gallai
1.596 + /// decomposition.
1.597 + Status decomposition(const Node& n) const {
1.598 + return (*_status)[n];
1.599 + }
1.600 +
1.601 + /// \brief Returns true when the node is in the barrier.
1.602 + ///
1.603 + /// Returns true when the node is in the barrier.
1.604 + bool barrier(const Node& n) const {
1.605 + return (*_status)[n] == ODD;
1.606 + }
1.607 +
1.608 + /// @}
1.609 +
1.610 + };
1.611 +
1.612 + /// \ingroup matching
1.613 + ///
1.614 + /// \brief Weighted matching in general graphs
1.615 + ///
1.616 + /// This class provides an efficient implementation of Edmond's
1.617 + /// maximum weighted matching algorithm. The implementation is based
1.618 + /// on extensive use of priority queues and provides
1.619 + /// \f$O(nm\log(n))\f$ time complexity.
1.620 + ///
1.621 + /// The maximum weighted matching problem is to find undirected
1.622 + /// edges in the graph with maximum overall weight and no two of
1.623 + /// them shares their ends. The problem can be formulated with the
1.624 + /// following linear program.
1.625 + /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
1.626 + /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
1.627 + \quad \forall B\in\mathcal{O}\f] */
1.628 + /// \f[x_e \ge 0\quad \forall e\in E\f]
1.629 + /// \f[\max \sum_{e\in E}x_ew_e\f]
1.630 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1.631 + /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
1.632 + /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
1.633 + /// subsets of the nodes.
1.634 + ///
1.635 + /// The algorithm calculates an optimal matching and a proof of the
1.636 + /// optimality. The solution of the dual problem can be used to check
1.637 + /// the result of the algorithm. The dual linear problem is the
1.638 + /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
1.639 + z_B \ge w_{uv} \quad \forall uv\in E\f] */
1.640 + /// \f[y_u \ge 0 \quad \forall u \in V\f]
1.641 + /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
1.642 + /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
1.643 + \frac{\vert B \vert - 1}{2}z_B\f] */
1.644 + ///
1.645 + /// The algorithm can be executed with \c run() or the \c init() and
1.646 + /// then the \c start() member functions. After it the matching can
1.647 + /// be asked with \c matching() or mate() functions. The dual
1.648 + /// solution can be get with \c nodeValue(), \c blossomNum() and \c
1.649 + /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
1.650 + /// "BlossomIt" nested class, which is able to iterate on the nodes
1.651 + /// of a blossom. If the value type is integral then the dual
1.652 + /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
1.653 + template <typename _Graph,
1.654 + typename _WeightMap = typename _Graph::template EdgeMap<int> >
1.655 + class MaxWeightedMatching {
1.656 + public:
1.657 +
1.658 + typedef _Graph Graph;
1.659 + typedef _WeightMap WeightMap;
1.660 + typedef typename WeightMap::Value Value;
1.661 +
1.662 + /// \brief Scaling factor for dual solution
1.663 + ///
1.664 + /// Scaling factor for dual solution, it is equal to 4 or 1
1.665 + /// according to the value type.
1.666 + static const int dualScale =
1.667 + std::numeric_limits<Value>::is_integer ? 4 : 1;
1.668 +
1.669 + typedef typename Graph::template NodeMap<typename Graph::Arc>
1.670 + MatchingMap;
1.671 +
1.672 + private:
1.673 +
1.674 + TEMPLATE_GRAPH_TYPEDEFS(Graph);
1.675 +
1.676 + typedef typename Graph::template NodeMap<Value> NodePotential;
1.677 + typedef std::vector<Node> BlossomNodeList;
1.678 +
1.679 + struct BlossomVariable {
1.680 + int begin, end;
1.681 + Value value;
1.682 +
1.683 + BlossomVariable(int _begin, int _end, Value _value)
1.684 + : begin(_begin), end(_end), value(_value) {}
1.685 +
1.686 + };
1.687 +
1.688 + typedef std::vector<BlossomVariable> BlossomPotential;
1.689 +
1.690 + const Graph& _graph;
1.691 + const WeightMap& _weight;
1.692 +
1.693 + MatchingMap* _matching;
1.694 +
1.695 + NodePotential* _node_potential;
1.696 +
1.697 + BlossomPotential _blossom_potential;
1.698 + BlossomNodeList _blossom_node_list;
1.699 +
1.700 + int _node_num;
1.701 + int _blossom_num;
1.702 +
1.703 + typedef RangeMap<int> IntIntMap;
1.704 +
1.705 + enum Status {
1.706 + EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
1.707 + };
1.708 +
1.709 + typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
1.710 + struct BlossomData {
1.711 + int tree;
1.712 + Status status;
1.713 + Arc pred, next;
1.714 + Value pot, offset;
1.715 + Node base;
1.716 + };
1.717 +
1.718 + IntNodeMap *_blossom_index;
1.719 + BlossomSet *_blossom_set;
1.720 + RangeMap<BlossomData>* _blossom_data;
1.721 +
1.722 + IntNodeMap *_node_index;
1.723 + IntArcMap *_node_heap_index;
1.724 +
1.725 + struct NodeData {
1.726 +
1.727 + NodeData(IntArcMap& node_heap_index)
1.728 + : heap(node_heap_index) {}
1.729 +
1.730 + int blossom;
1.731 + Value pot;
1.732 + BinHeap<Value, IntArcMap> heap;
1.733 + std::map<int, Arc> heap_index;
1.734 +
1.735 + int tree;
1.736 + };
1.737 +
1.738 + RangeMap<NodeData>* _node_data;
1.739 +
1.740 + typedef ExtendFindEnum<IntIntMap> TreeSet;
1.741 +
1.742 + IntIntMap *_tree_set_index;
1.743 + TreeSet *_tree_set;
1.744 +
1.745 + IntNodeMap *_delta1_index;
1.746 + BinHeap<Value, IntNodeMap> *_delta1;
1.747 +
1.748 + IntIntMap *_delta2_index;
1.749 + BinHeap<Value, IntIntMap> *_delta2;
1.750 +
1.751 + IntEdgeMap *_delta3_index;
1.752 + BinHeap<Value, IntEdgeMap> *_delta3;
1.753 +
1.754 + IntIntMap *_delta4_index;
1.755 + BinHeap<Value, IntIntMap> *_delta4;
1.756 +
1.757 + Value _delta_sum;
1.758 +
1.759 + void createStructures() {
1.760 + _node_num = countNodes(_graph);
1.761 + _blossom_num = _node_num * 3 / 2;
1.762 +
1.763 + if (!_matching) {
1.764 + _matching = new MatchingMap(_graph);
1.765 + }
1.766 + if (!_node_potential) {
1.767 + _node_potential = new NodePotential(_graph);
1.768 + }
1.769 + if (!_blossom_set) {
1.770 + _blossom_index = new IntNodeMap(_graph);
1.771 + _blossom_set = new BlossomSet(*_blossom_index);
1.772 + _blossom_data = new RangeMap<BlossomData>(_blossom_num);
1.773 + }
1.774 +
1.775 + if (!_node_index) {
1.776 + _node_index = new IntNodeMap(_graph);
1.777 + _node_heap_index = new IntArcMap(_graph);
1.778 + _node_data = new RangeMap<NodeData>(_node_num,
1.779 + NodeData(*_node_heap_index));
1.780 + }
1.781 +
1.782 + if (!_tree_set) {
1.783 + _tree_set_index = new IntIntMap(_blossom_num);
1.784 + _tree_set = new TreeSet(*_tree_set_index);
1.785 + }
1.786 + if (!_delta1) {
1.787 + _delta1_index = new IntNodeMap(_graph);
1.788 + _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
1.789 + }
1.790 + if (!_delta2) {
1.791 + _delta2_index = new IntIntMap(_blossom_num);
1.792 + _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
1.793 + }
1.794 + if (!_delta3) {
1.795 + _delta3_index = new IntEdgeMap(_graph);
1.796 + _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
1.797 + }
1.798 + if (!_delta4) {
1.799 + _delta4_index = new IntIntMap(_blossom_num);
1.800 + _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
1.801 + }
1.802 + }
1.803 +
1.804 + void destroyStructures() {
1.805 + _node_num = countNodes(_graph);
1.806 + _blossom_num = _node_num * 3 / 2;
1.807 +
1.808 + if (_matching) {
1.809 + delete _matching;
1.810 + }
1.811 + if (_node_potential) {
1.812 + delete _node_potential;
1.813 + }
1.814 + if (_blossom_set) {
1.815 + delete _blossom_index;
1.816 + delete _blossom_set;
1.817 + delete _blossom_data;
1.818 + }
1.819 +
1.820 + if (_node_index) {
1.821 + delete _node_index;
1.822 + delete _node_heap_index;
1.823 + delete _node_data;
1.824 + }
1.825 +
1.826 + if (_tree_set) {
1.827 + delete _tree_set_index;
1.828 + delete _tree_set;
1.829 + }
1.830 + if (_delta1) {
1.831 + delete _delta1_index;
1.832 + delete _delta1;
1.833 + }
1.834 + if (_delta2) {
1.835 + delete _delta2_index;
1.836 + delete _delta2;
1.837 + }
1.838 + if (_delta3) {
1.839 + delete _delta3_index;
1.840 + delete _delta3;
1.841 + }
1.842 + if (_delta4) {
1.843 + delete _delta4_index;
1.844 + delete _delta4;
1.845 + }
1.846 + }
1.847 +
1.848 + void matchedToEven(int blossom, int tree) {
1.849 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.850 + _delta2->erase(blossom);
1.851 + }
1.852 +
1.853 + if (!_blossom_set->trivial(blossom)) {
1.854 + (*_blossom_data)[blossom].pot -=
1.855 + 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
1.856 + }
1.857 +
1.858 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.859 + n != INVALID; ++n) {
1.860 +
1.861 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
1.862 + int ni = (*_node_index)[n];
1.863 +
1.864 + (*_node_data)[ni].heap.clear();
1.865 + (*_node_data)[ni].heap_index.clear();
1.866 +
1.867 + (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
1.868 +
1.869 + _delta1->push(n, (*_node_data)[ni].pot);
1.870 +
1.871 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.872 + Node v = _graph.source(e);
1.873 + int vb = _blossom_set->find(v);
1.874 + int vi = (*_node_index)[v];
1.875 +
1.876 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.877 + dualScale * _weight[e];
1.878 +
1.879 + if ((*_blossom_data)[vb].status == EVEN) {
1.880 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1.881 + _delta3->push(e, rw / 2);
1.882 + }
1.883 + } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1.884 + if (_delta3->state(e) != _delta3->IN_HEAP) {
1.885 + _delta3->push(e, rw);
1.886 + }
1.887 + } else {
1.888 + typename std::map<int, Arc>::iterator it =
1.889 + (*_node_data)[vi].heap_index.find(tree);
1.890 +
1.891 + if (it != (*_node_data)[vi].heap_index.end()) {
1.892 + if ((*_node_data)[vi].heap[it->second] > rw) {
1.893 + (*_node_data)[vi].heap.replace(it->second, e);
1.894 + (*_node_data)[vi].heap.decrease(e, rw);
1.895 + it->second = e;
1.896 + }
1.897 + } else {
1.898 + (*_node_data)[vi].heap.push(e, rw);
1.899 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1.900 + }
1.901 +
1.902 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1.903 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1.904 +
1.905 + if ((*_blossom_data)[vb].status == MATCHED) {
1.906 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
1.907 + _delta2->push(vb, _blossom_set->classPrio(vb) -
1.908 + (*_blossom_data)[vb].offset);
1.909 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1.910 + (*_blossom_data)[vb].offset){
1.911 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1.912 + (*_blossom_data)[vb].offset);
1.913 + }
1.914 + }
1.915 + }
1.916 + }
1.917 + }
1.918 + }
1.919 + (*_blossom_data)[blossom].offset = 0;
1.920 + }
1.921 +
1.922 + void matchedToOdd(int blossom) {
1.923 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.924 + _delta2->erase(blossom);
1.925 + }
1.926 + (*_blossom_data)[blossom].offset += _delta_sum;
1.927 + if (!_blossom_set->trivial(blossom)) {
1.928 + _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
1.929 + (*_blossom_data)[blossom].offset);
1.930 + }
1.931 + }
1.932 +
1.933 + void evenToMatched(int blossom, int tree) {
1.934 + if (!_blossom_set->trivial(blossom)) {
1.935 + (*_blossom_data)[blossom].pot += 2 * _delta_sum;
1.936 + }
1.937 +
1.938 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.939 + n != INVALID; ++n) {
1.940 + int ni = (*_node_index)[n];
1.941 + (*_node_data)[ni].pot -= _delta_sum;
1.942 +
1.943 + _delta1->erase(n);
1.944 +
1.945 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.946 + Node v = _graph.source(e);
1.947 + int vb = _blossom_set->find(v);
1.948 + int vi = (*_node_index)[v];
1.949 +
1.950 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.951 + dualScale * _weight[e];
1.952 +
1.953 + if (vb == blossom) {
1.954 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.955 + _delta3->erase(e);
1.956 + }
1.957 + } else if ((*_blossom_data)[vb].status == EVEN) {
1.958 +
1.959 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.960 + _delta3->erase(e);
1.961 + }
1.962 +
1.963 + int vt = _tree_set->find(vb);
1.964 +
1.965 + if (vt != tree) {
1.966 +
1.967 + Arc r = _graph.oppositeArc(e);
1.968 +
1.969 + typename std::map<int, Arc>::iterator it =
1.970 + (*_node_data)[ni].heap_index.find(vt);
1.971 +
1.972 + if (it != (*_node_data)[ni].heap_index.end()) {
1.973 + if ((*_node_data)[ni].heap[it->second] > rw) {
1.974 + (*_node_data)[ni].heap.replace(it->second, r);
1.975 + (*_node_data)[ni].heap.decrease(r, rw);
1.976 + it->second = r;
1.977 + }
1.978 + } else {
1.979 + (*_node_data)[ni].heap.push(r, rw);
1.980 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1.981 + }
1.982 +
1.983 + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1.984 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1.985 +
1.986 + if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1.987 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.988 + (*_blossom_data)[blossom].offset);
1.989 + } else if ((*_delta2)[blossom] >
1.990 + _blossom_set->classPrio(blossom) -
1.991 + (*_blossom_data)[blossom].offset){
1.992 + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1.993 + (*_blossom_data)[blossom].offset);
1.994 + }
1.995 + }
1.996 + }
1.997 +
1.998 + } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1.999 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.1000 + _delta3->erase(e);
1.1001 + }
1.1002 + } else {
1.1003 +
1.1004 + typename std::map<int, Arc>::iterator it =
1.1005 + (*_node_data)[vi].heap_index.find(tree);
1.1006 +
1.1007 + if (it != (*_node_data)[vi].heap_index.end()) {
1.1008 + (*_node_data)[vi].heap.erase(it->second);
1.1009 + (*_node_data)[vi].heap_index.erase(it);
1.1010 + if ((*_node_data)[vi].heap.empty()) {
1.1011 + _blossom_set->increase(v, std::numeric_limits<Value>::max());
1.1012 + } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1.1013 + _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1.1014 + }
1.1015 +
1.1016 + if ((*_blossom_data)[vb].status == MATCHED) {
1.1017 + if (_blossom_set->classPrio(vb) ==
1.1018 + std::numeric_limits<Value>::max()) {
1.1019 + _delta2->erase(vb);
1.1020 + } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1.1021 + (*_blossom_data)[vb].offset) {
1.1022 + _delta2->increase(vb, _blossom_set->classPrio(vb) -
1.1023 + (*_blossom_data)[vb].offset);
1.1024 + }
1.1025 + }
1.1026 + }
1.1027 + }
1.1028 + }
1.1029 + }
1.1030 + }
1.1031 +
1.1032 + void oddToMatched(int blossom) {
1.1033 + (*_blossom_data)[blossom].offset -= _delta_sum;
1.1034 +
1.1035 + if (_blossom_set->classPrio(blossom) !=
1.1036 + std::numeric_limits<Value>::max()) {
1.1037 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.1038 + (*_blossom_data)[blossom].offset);
1.1039 + }
1.1040 +
1.1041 + if (!_blossom_set->trivial(blossom)) {
1.1042 + _delta4->erase(blossom);
1.1043 + }
1.1044 + }
1.1045 +
1.1046 + void oddToEven(int blossom, int tree) {
1.1047 + if (!_blossom_set->trivial(blossom)) {
1.1048 + _delta4->erase(blossom);
1.1049 + (*_blossom_data)[blossom].pot -=
1.1050 + 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1.1051 + }
1.1052 +
1.1053 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.1054 + n != INVALID; ++n) {
1.1055 + int ni = (*_node_index)[n];
1.1056 +
1.1057 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
1.1058 +
1.1059 + (*_node_data)[ni].heap.clear();
1.1060 + (*_node_data)[ni].heap_index.clear();
1.1061 + (*_node_data)[ni].pot +=
1.1062 + 2 * _delta_sum - (*_blossom_data)[blossom].offset;
1.1063 +
1.1064 + _delta1->push(n, (*_node_data)[ni].pot);
1.1065 +
1.1066 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.1067 + Node v = _graph.source(e);
1.1068 + int vb = _blossom_set->find(v);
1.1069 + int vi = (*_node_index)[v];
1.1070 +
1.1071 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.1072 + dualScale * _weight[e];
1.1073 +
1.1074 + if ((*_blossom_data)[vb].status == EVEN) {
1.1075 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1.1076 + _delta3->push(e, rw / 2);
1.1077 + }
1.1078 + } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1.1079 + if (_delta3->state(e) != _delta3->IN_HEAP) {
1.1080 + _delta3->push(e, rw);
1.1081 + }
1.1082 + } else {
1.1083 +
1.1084 + typename std::map<int, Arc>::iterator it =
1.1085 + (*_node_data)[vi].heap_index.find(tree);
1.1086 +
1.1087 + if (it != (*_node_data)[vi].heap_index.end()) {
1.1088 + if ((*_node_data)[vi].heap[it->second] > rw) {
1.1089 + (*_node_data)[vi].heap.replace(it->second, e);
1.1090 + (*_node_data)[vi].heap.decrease(e, rw);
1.1091 + it->second = e;
1.1092 + }
1.1093 + } else {
1.1094 + (*_node_data)[vi].heap.push(e, rw);
1.1095 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1.1096 + }
1.1097 +
1.1098 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1.1099 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1.1100 +
1.1101 + if ((*_blossom_data)[vb].status == MATCHED) {
1.1102 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
1.1103 + _delta2->push(vb, _blossom_set->classPrio(vb) -
1.1104 + (*_blossom_data)[vb].offset);
1.1105 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1.1106 + (*_blossom_data)[vb].offset) {
1.1107 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1.1108 + (*_blossom_data)[vb].offset);
1.1109 + }
1.1110 + }
1.1111 + }
1.1112 + }
1.1113 + }
1.1114 + }
1.1115 + (*_blossom_data)[blossom].offset = 0;
1.1116 + }
1.1117 +
1.1118 +
1.1119 + void matchedToUnmatched(int blossom) {
1.1120 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.1121 + _delta2->erase(blossom);
1.1122 + }
1.1123 +
1.1124 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.1125 + n != INVALID; ++n) {
1.1126 + int ni = (*_node_index)[n];
1.1127 +
1.1128 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
1.1129 +
1.1130 + (*_node_data)[ni].heap.clear();
1.1131 + (*_node_data)[ni].heap_index.clear();
1.1132 +
1.1133 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.1134 + Node v = _graph.target(e);
1.1135 + int vb = _blossom_set->find(v);
1.1136 + int vi = (*_node_index)[v];
1.1137 +
1.1138 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.1139 + dualScale * _weight[e];
1.1140 +
1.1141 + if ((*_blossom_data)[vb].status == EVEN) {
1.1142 + if (_delta3->state(e) != _delta3->IN_HEAP) {
1.1143 + _delta3->push(e, rw);
1.1144 + }
1.1145 + }
1.1146 + }
1.1147 + }
1.1148 + }
1.1149 +
1.1150 + void unmatchedToMatched(int blossom) {
1.1151 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.1152 + n != INVALID; ++n) {
1.1153 + int ni = (*_node_index)[n];
1.1154 +
1.1155 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.1156 + Node v = _graph.source(e);
1.1157 + int vb = _blossom_set->find(v);
1.1158 + int vi = (*_node_index)[v];
1.1159 +
1.1160 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.1161 + dualScale * _weight[e];
1.1162 +
1.1163 + if (vb == blossom) {
1.1164 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.1165 + _delta3->erase(e);
1.1166 + }
1.1167 + } else if ((*_blossom_data)[vb].status == EVEN) {
1.1168 +
1.1169 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.1170 + _delta3->erase(e);
1.1171 + }
1.1172 +
1.1173 + int vt = _tree_set->find(vb);
1.1174 +
1.1175 + Arc r = _graph.oppositeArc(e);
1.1176 +
1.1177 + typename std::map<int, Arc>::iterator it =
1.1178 + (*_node_data)[ni].heap_index.find(vt);
1.1179 +
1.1180 + if (it != (*_node_data)[ni].heap_index.end()) {
1.1181 + if ((*_node_data)[ni].heap[it->second] > rw) {
1.1182 + (*_node_data)[ni].heap.replace(it->second, r);
1.1183 + (*_node_data)[ni].heap.decrease(r, rw);
1.1184 + it->second = r;
1.1185 + }
1.1186 + } else {
1.1187 + (*_node_data)[ni].heap.push(r, rw);
1.1188 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1.1189 + }
1.1190 +
1.1191 + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1.1192 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1.1193 +
1.1194 + if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1.1195 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.1196 + (*_blossom_data)[blossom].offset);
1.1197 + } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
1.1198 + (*_blossom_data)[blossom].offset){
1.1199 + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1.1200 + (*_blossom_data)[blossom].offset);
1.1201 + }
1.1202 + }
1.1203 +
1.1204 + } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1.1205 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.1206 + _delta3->erase(e);
1.1207 + }
1.1208 + }
1.1209 + }
1.1210 + }
1.1211 + }
1.1212 +
1.1213 + void alternatePath(int even, int tree) {
1.1214 + int odd;
1.1215 +
1.1216 + evenToMatched(even, tree);
1.1217 + (*_blossom_data)[even].status = MATCHED;
1.1218 +
1.1219 + while ((*_blossom_data)[even].pred != INVALID) {
1.1220 + odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
1.1221 + (*_blossom_data)[odd].status = MATCHED;
1.1222 + oddToMatched(odd);
1.1223 + (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
1.1224 +
1.1225 + even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
1.1226 + (*_blossom_data)[even].status = MATCHED;
1.1227 + evenToMatched(even, tree);
1.1228 + (*_blossom_data)[even].next =
1.1229 + _graph.oppositeArc((*_blossom_data)[odd].pred);
1.1230 + }
1.1231 +
1.1232 + }
1.1233 +
1.1234 + void destroyTree(int tree) {
1.1235 + for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
1.1236 + if ((*_blossom_data)[b].status == EVEN) {
1.1237 + (*_blossom_data)[b].status = MATCHED;
1.1238 + evenToMatched(b, tree);
1.1239 + } else if ((*_blossom_data)[b].status == ODD) {
1.1240 + (*_blossom_data)[b].status = MATCHED;
1.1241 + oddToMatched(b);
1.1242 + }
1.1243 + }
1.1244 + _tree_set->eraseClass(tree);
1.1245 + }
1.1246 +
1.1247 +
1.1248 + void unmatchNode(const Node& node) {
1.1249 + int blossom = _blossom_set->find(node);
1.1250 + int tree = _tree_set->find(blossom);
1.1251 +
1.1252 + alternatePath(blossom, tree);
1.1253 + destroyTree(tree);
1.1254 +
1.1255 + (*_blossom_data)[blossom].status = UNMATCHED;
1.1256 + (*_blossom_data)[blossom].base = node;
1.1257 + matchedToUnmatched(blossom);
1.1258 + }
1.1259 +
1.1260 +
1.1261 + void augmentOnEdge(const Edge& edge) {
1.1262 +
1.1263 + int left = _blossom_set->find(_graph.u(edge));
1.1264 + int right = _blossom_set->find(_graph.v(edge));
1.1265 +
1.1266 + if ((*_blossom_data)[left].status == EVEN) {
1.1267 + int left_tree = _tree_set->find(left);
1.1268 + alternatePath(left, left_tree);
1.1269 + destroyTree(left_tree);
1.1270 + } else {
1.1271 + (*_blossom_data)[left].status = MATCHED;
1.1272 + unmatchedToMatched(left);
1.1273 + }
1.1274 +
1.1275 + if ((*_blossom_data)[right].status == EVEN) {
1.1276 + int right_tree = _tree_set->find(right);
1.1277 + alternatePath(right, right_tree);
1.1278 + destroyTree(right_tree);
1.1279 + } else {
1.1280 + (*_blossom_data)[right].status = MATCHED;
1.1281 + unmatchedToMatched(right);
1.1282 + }
1.1283 +
1.1284 + (*_blossom_data)[left].next = _graph.direct(edge, true);
1.1285 + (*_blossom_data)[right].next = _graph.direct(edge, false);
1.1286 + }
1.1287 +
1.1288 + void extendOnArc(const Arc& arc) {
1.1289 + int base = _blossom_set->find(_graph.target(arc));
1.1290 + int tree = _tree_set->find(base);
1.1291 +
1.1292 + int odd = _blossom_set->find(_graph.source(arc));
1.1293 + _tree_set->insert(odd, tree);
1.1294 + (*_blossom_data)[odd].status = ODD;
1.1295 + matchedToOdd(odd);
1.1296 + (*_blossom_data)[odd].pred = arc;
1.1297 +
1.1298 + int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
1.1299 + (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
1.1300 + _tree_set->insert(even, tree);
1.1301 + (*_blossom_data)[even].status = EVEN;
1.1302 + matchedToEven(even, tree);
1.1303 + }
1.1304 +
1.1305 + void shrinkOnEdge(const Edge& edge, int tree) {
1.1306 + int nca = -1;
1.1307 + std::vector<int> left_path, right_path;
1.1308 +
1.1309 + {
1.1310 + std::set<int> left_set, right_set;
1.1311 + int left = _blossom_set->find(_graph.u(edge));
1.1312 + left_path.push_back(left);
1.1313 + left_set.insert(left);
1.1314 +
1.1315 + int right = _blossom_set->find(_graph.v(edge));
1.1316 + right_path.push_back(right);
1.1317 + right_set.insert(right);
1.1318 +
1.1319 + while (true) {
1.1320 +
1.1321 + if ((*_blossom_data)[left].pred == INVALID) break;
1.1322 +
1.1323 + left =
1.1324 + _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1.1325 + left_path.push_back(left);
1.1326 + left =
1.1327 + _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1.1328 + left_path.push_back(left);
1.1329 +
1.1330 + left_set.insert(left);
1.1331 +
1.1332 + if (right_set.find(left) != right_set.end()) {
1.1333 + nca = left;
1.1334 + break;
1.1335 + }
1.1336 +
1.1337 + if ((*_blossom_data)[right].pred == INVALID) break;
1.1338 +
1.1339 + right =
1.1340 + _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1.1341 + right_path.push_back(right);
1.1342 + right =
1.1343 + _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1.1344 + right_path.push_back(right);
1.1345 +
1.1346 + right_set.insert(right);
1.1347 +
1.1348 + if (left_set.find(right) != left_set.end()) {
1.1349 + nca = right;
1.1350 + break;
1.1351 + }
1.1352 +
1.1353 + }
1.1354 +
1.1355 + if (nca == -1) {
1.1356 + if ((*_blossom_data)[left].pred == INVALID) {
1.1357 + nca = right;
1.1358 + while (left_set.find(nca) == left_set.end()) {
1.1359 + nca =
1.1360 + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1.1361 + right_path.push_back(nca);
1.1362 + nca =
1.1363 + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1.1364 + right_path.push_back(nca);
1.1365 + }
1.1366 + } else {
1.1367 + nca = left;
1.1368 + while (right_set.find(nca) == right_set.end()) {
1.1369 + nca =
1.1370 + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1.1371 + left_path.push_back(nca);
1.1372 + nca =
1.1373 + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1.1374 + left_path.push_back(nca);
1.1375 + }
1.1376 + }
1.1377 + }
1.1378 + }
1.1379 +
1.1380 + std::vector<int> subblossoms;
1.1381 + Arc prev;
1.1382 +
1.1383 + prev = _graph.direct(edge, true);
1.1384 + for (int i = 0; left_path[i] != nca; i += 2) {
1.1385 + subblossoms.push_back(left_path[i]);
1.1386 + (*_blossom_data)[left_path[i]].next = prev;
1.1387 + _tree_set->erase(left_path[i]);
1.1388 +
1.1389 + subblossoms.push_back(left_path[i + 1]);
1.1390 + (*_blossom_data)[left_path[i + 1]].status = EVEN;
1.1391 + oddToEven(left_path[i + 1], tree);
1.1392 + _tree_set->erase(left_path[i + 1]);
1.1393 + prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
1.1394 + }
1.1395 +
1.1396 + int k = 0;
1.1397 + while (right_path[k] != nca) ++k;
1.1398 +
1.1399 + subblossoms.push_back(nca);
1.1400 + (*_blossom_data)[nca].next = prev;
1.1401 +
1.1402 + for (int i = k - 2; i >= 0; i -= 2) {
1.1403 + subblossoms.push_back(right_path[i + 1]);
1.1404 + (*_blossom_data)[right_path[i + 1]].status = EVEN;
1.1405 + oddToEven(right_path[i + 1], tree);
1.1406 + _tree_set->erase(right_path[i + 1]);
1.1407 +
1.1408 + (*_blossom_data)[right_path[i + 1]].next =
1.1409 + (*_blossom_data)[right_path[i + 1]].pred;
1.1410 +
1.1411 + subblossoms.push_back(right_path[i]);
1.1412 + _tree_set->erase(right_path[i]);
1.1413 + }
1.1414 +
1.1415 + int surface =
1.1416 + _blossom_set->join(subblossoms.begin(), subblossoms.end());
1.1417 +
1.1418 + for (int i = 0; i < int(subblossoms.size()); ++i) {
1.1419 + if (!_blossom_set->trivial(subblossoms[i])) {
1.1420 + (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
1.1421 + }
1.1422 + (*_blossom_data)[subblossoms[i]].status = MATCHED;
1.1423 + }
1.1424 +
1.1425 + (*_blossom_data)[surface].pot = -2 * _delta_sum;
1.1426 + (*_blossom_data)[surface].offset = 0;
1.1427 + (*_blossom_data)[surface].status = EVEN;
1.1428 + (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
1.1429 + (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
1.1430 +
1.1431 + _tree_set->insert(surface, tree);
1.1432 + _tree_set->erase(nca);
1.1433 + }
1.1434 +
1.1435 + void splitBlossom(int blossom) {
1.1436 + Arc next = (*_blossom_data)[blossom].next;
1.1437 + Arc pred = (*_blossom_data)[blossom].pred;
1.1438 +
1.1439 + int tree = _tree_set->find(blossom);
1.1440 +
1.1441 + (*_blossom_data)[blossom].status = MATCHED;
1.1442 + oddToMatched(blossom);
1.1443 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.1444 + _delta2->erase(blossom);
1.1445 + }
1.1446 +
1.1447 + std::vector<int> subblossoms;
1.1448 + _blossom_set->split(blossom, std::back_inserter(subblossoms));
1.1449 +
1.1450 + Value offset = (*_blossom_data)[blossom].offset;
1.1451 + int b = _blossom_set->find(_graph.source(pred));
1.1452 + int d = _blossom_set->find(_graph.source(next));
1.1453 +
1.1454 + int ib = -1, id = -1;
1.1455 + for (int i = 0; i < int(subblossoms.size()); ++i) {
1.1456 + if (subblossoms[i] == b) ib = i;
1.1457 + if (subblossoms[i] == d) id = i;
1.1458 +
1.1459 + (*_blossom_data)[subblossoms[i]].offset = offset;
1.1460 + if (!_blossom_set->trivial(subblossoms[i])) {
1.1461 + (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
1.1462 + }
1.1463 + if (_blossom_set->classPrio(subblossoms[i]) !=
1.1464 + std::numeric_limits<Value>::max()) {
1.1465 + _delta2->push(subblossoms[i],
1.1466 + _blossom_set->classPrio(subblossoms[i]) -
1.1467 + (*_blossom_data)[subblossoms[i]].offset);
1.1468 + }
1.1469 + }
1.1470 +
1.1471 + if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
1.1472 + for (int i = (id + 1) % subblossoms.size();
1.1473 + i != ib; i = (i + 2) % subblossoms.size()) {
1.1474 + int sb = subblossoms[i];
1.1475 + int tb = subblossoms[(i + 1) % subblossoms.size()];
1.1476 + (*_blossom_data)[sb].next =
1.1477 + _graph.oppositeArc((*_blossom_data)[tb].next);
1.1478 + }
1.1479 +
1.1480 + for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
1.1481 + int sb = subblossoms[i];
1.1482 + int tb = subblossoms[(i + 1) % subblossoms.size()];
1.1483 + int ub = subblossoms[(i + 2) % subblossoms.size()];
1.1484 +
1.1485 + (*_blossom_data)[sb].status = ODD;
1.1486 + matchedToOdd(sb);
1.1487 + _tree_set->insert(sb, tree);
1.1488 + (*_blossom_data)[sb].pred = pred;
1.1489 + (*_blossom_data)[sb].next =
1.1490 + _graph.oppositeArc((*_blossom_data)[tb].next);
1.1491 +
1.1492 + pred = (*_blossom_data)[ub].next;
1.1493 +
1.1494 + (*_blossom_data)[tb].status = EVEN;
1.1495 + matchedToEven(tb, tree);
1.1496 + _tree_set->insert(tb, tree);
1.1497 + (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
1.1498 + }
1.1499 +
1.1500 + (*_blossom_data)[subblossoms[id]].status = ODD;
1.1501 + matchedToOdd(subblossoms[id]);
1.1502 + _tree_set->insert(subblossoms[id], tree);
1.1503 + (*_blossom_data)[subblossoms[id]].next = next;
1.1504 + (*_blossom_data)[subblossoms[id]].pred = pred;
1.1505 +
1.1506 + } else {
1.1507 +
1.1508 + for (int i = (ib + 1) % subblossoms.size();
1.1509 + i != id; i = (i + 2) % subblossoms.size()) {
1.1510 + int sb = subblossoms[i];
1.1511 + int tb = subblossoms[(i + 1) % subblossoms.size()];
1.1512 + (*_blossom_data)[sb].next =
1.1513 + _graph.oppositeArc((*_blossom_data)[tb].next);
1.1514 + }
1.1515 +
1.1516 + for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
1.1517 + int sb = subblossoms[i];
1.1518 + int tb = subblossoms[(i + 1) % subblossoms.size()];
1.1519 + int ub = subblossoms[(i + 2) % subblossoms.size()];
1.1520 +
1.1521 + (*_blossom_data)[sb].status = ODD;
1.1522 + matchedToOdd(sb);
1.1523 + _tree_set->insert(sb, tree);
1.1524 + (*_blossom_data)[sb].next = next;
1.1525 + (*_blossom_data)[sb].pred =
1.1526 + _graph.oppositeArc((*_blossom_data)[tb].next);
1.1527 +
1.1528 + (*_blossom_data)[tb].status = EVEN;
1.1529 + matchedToEven(tb, tree);
1.1530 + _tree_set->insert(tb, tree);
1.1531 + (*_blossom_data)[tb].pred =
1.1532 + (*_blossom_data)[tb].next =
1.1533 + _graph.oppositeArc((*_blossom_data)[ub].next);
1.1534 + next = (*_blossom_data)[ub].next;
1.1535 + }
1.1536 +
1.1537 + (*_blossom_data)[subblossoms[ib]].status = ODD;
1.1538 + matchedToOdd(subblossoms[ib]);
1.1539 + _tree_set->insert(subblossoms[ib], tree);
1.1540 + (*_blossom_data)[subblossoms[ib]].next = next;
1.1541 + (*_blossom_data)[subblossoms[ib]].pred = pred;
1.1542 + }
1.1543 + _tree_set->erase(blossom);
1.1544 + }
1.1545 +
1.1546 + void extractBlossom(int blossom, const Node& base, const Arc& matching) {
1.1547 + if (_blossom_set->trivial(blossom)) {
1.1548 + int bi = (*_node_index)[base];
1.1549 + Value pot = (*_node_data)[bi].pot;
1.1550 +
1.1551 + _matching->set(base, matching);
1.1552 + _blossom_node_list.push_back(base);
1.1553 + _node_potential->set(base, pot);
1.1554 + } else {
1.1555 +
1.1556 + Value pot = (*_blossom_data)[blossom].pot;
1.1557 + int bn = _blossom_node_list.size();
1.1558 +
1.1559 + std::vector<int> subblossoms;
1.1560 + _blossom_set->split(blossom, std::back_inserter(subblossoms));
1.1561 + int b = _blossom_set->find(base);
1.1562 + int ib = -1;
1.1563 + for (int i = 0; i < int(subblossoms.size()); ++i) {
1.1564 + if (subblossoms[i] == b) { ib = i; break; }
1.1565 + }
1.1566 +
1.1567 + for (int i = 1; i < int(subblossoms.size()); i += 2) {
1.1568 + int sb = subblossoms[(ib + i) % subblossoms.size()];
1.1569 + int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
1.1570 +
1.1571 + Arc m = (*_blossom_data)[tb].next;
1.1572 + extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
1.1573 + extractBlossom(tb, _graph.source(m), m);
1.1574 + }
1.1575 + extractBlossom(subblossoms[ib], base, matching);
1.1576 +
1.1577 + int en = _blossom_node_list.size();
1.1578 +
1.1579 + _blossom_potential.push_back(BlossomVariable(bn, en, pot));
1.1580 + }
1.1581 + }
1.1582 +
1.1583 + void extractMatching() {
1.1584 + std::vector<int> blossoms;
1.1585 + for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1.1586 + blossoms.push_back(c);
1.1587 + }
1.1588 +
1.1589 + for (int i = 0; i < int(blossoms.size()); ++i) {
1.1590 + if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
1.1591 +
1.1592 + Value offset = (*_blossom_data)[blossoms[i]].offset;
1.1593 + (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1.1594 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1.1595 + n != INVALID; ++n) {
1.1596 + (*_node_data)[(*_node_index)[n]].pot -= offset;
1.1597 + }
1.1598 +
1.1599 + Arc matching = (*_blossom_data)[blossoms[i]].next;
1.1600 + Node base = _graph.source(matching);
1.1601 + extractBlossom(blossoms[i], base, matching);
1.1602 + } else {
1.1603 + Node base = (*_blossom_data)[blossoms[i]].base;
1.1604 + extractBlossom(blossoms[i], base, INVALID);
1.1605 + }
1.1606 + }
1.1607 + }
1.1608 +
1.1609 + public:
1.1610 +
1.1611 + /// \brief Constructor
1.1612 + ///
1.1613 + /// Constructor.
1.1614 + MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
1.1615 + : _graph(graph), _weight(weight), _matching(0),
1.1616 + _node_potential(0), _blossom_potential(), _blossom_node_list(),
1.1617 + _node_num(0), _blossom_num(0),
1.1618 +
1.1619 + _blossom_index(0), _blossom_set(0), _blossom_data(0),
1.1620 + _node_index(0), _node_heap_index(0), _node_data(0),
1.1621 + _tree_set_index(0), _tree_set(0),
1.1622 +
1.1623 + _delta1_index(0), _delta1(0),
1.1624 + _delta2_index(0), _delta2(0),
1.1625 + _delta3_index(0), _delta3(0),
1.1626 + _delta4_index(0), _delta4(0),
1.1627 +
1.1628 + _delta_sum() {}
1.1629 +
1.1630 + ~MaxWeightedMatching() {
1.1631 + destroyStructures();
1.1632 + }
1.1633 +
1.1634 + /// \name Execution control
1.1635 + /// The simplest way to execute the algorithm is to use the
1.1636 + /// \c run() member function.
1.1637 +
1.1638 + ///@{
1.1639 +
1.1640 + /// \brief Initialize the algorithm
1.1641 + ///
1.1642 + /// Initialize the algorithm
1.1643 + void init() {
1.1644 + createStructures();
1.1645 +
1.1646 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.1647 + _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
1.1648 + }
1.1649 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1650 + _delta1_index->set(n, _delta1->PRE_HEAP);
1.1651 + }
1.1652 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.1653 + _delta3_index->set(e, _delta3->PRE_HEAP);
1.1654 + }
1.1655 + for (int i = 0; i < _blossom_num; ++i) {
1.1656 + _delta2_index->set(i, _delta2->PRE_HEAP);
1.1657 + _delta4_index->set(i, _delta4->PRE_HEAP);
1.1658 + }
1.1659 +
1.1660 + int index = 0;
1.1661 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1662 + Value max = 0;
1.1663 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.1664 + if (_graph.target(e) == n) continue;
1.1665 + if ((dualScale * _weight[e]) / 2 > max) {
1.1666 + max = (dualScale * _weight[e]) / 2;
1.1667 + }
1.1668 + }
1.1669 + _node_index->set(n, index);
1.1670 + (*_node_data)[index].pot = max;
1.1671 + _delta1->push(n, max);
1.1672 + int blossom =
1.1673 + _blossom_set->insert(n, std::numeric_limits<Value>::max());
1.1674 +
1.1675 + _tree_set->insert(blossom);
1.1676 +
1.1677 + (*_blossom_data)[blossom].status = EVEN;
1.1678 + (*_blossom_data)[blossom].pred = INVALID;
1.1679 + (*_blossom_data)[blossom].next = INVALID;
1.1680 + (*_blossom_data)[blossom].pot = 0;
1.1681 + (*_blossom_data)[blossom].offset = 0;
1.1682 + ++index;
1.1683 + }
1.1684 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.1685 + int si = (*_node_index)[_graph.u(e)];
1.1686 + int ti = (*_node_index)[_graph.v(e)];
1.1687 + if (_graph.u(e) != _graph.v(e)) {
1.1688 + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1.1689 + dualScale * _weight[e]) / 2);
1.1690 + }
1.1691 + }
1.1692 + }
1.1693 +
1.1694 + /// \brief Starts the algorithm
1.1695 + ///
1.1696 + /// Starts the algorithm
1.1697 + void start() {
1.1698 + enum OpType {
1.1699 + D1, D2, D3, D4
1.1700 + };
1.1701 +
1.1702 + int unmatched = _node_num;
1.1703 + while (unmatched > 0) {
1.1704 + Value d1 = !_delta1->empty() ?
1.1705 + _delta1->prio() : std::numeric_limits<Value>::max();
1.1706 +
1.1707 + Value d2 = !_delta2->empty() ?
1.1708 + _delta2->prio() : std::numeric_limits<Value>::max();
1.1709 +
1.1710 + Value d3 = !_delta3->empty() ?
1.1711 + _delta3->prio() : std::numeric_limits<Value>::max();
1.1712 +
1.1713 + Value d4 = !_delta4->empty() ?
1.1714 + _delta4->prio() : std::numeric_limits<Value>::max();
1.1715 +
1.1716 + _delta_sum = d1; OpType ot = D1;
1.1717 + if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1.1718 + if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
1.1719 + if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1.1720 +
1.1721 +
1.1722 + switch (ot) {
1.1723 + case D1:
1.1724 + {
1.1725 + Node n = _delta1->top();
1.1726 + unmatchNode(n);
1.1727 + --unmatched;
1.1728 + }
1.1729 + break;
1.1730 + case D2:
1.1731 + {
1.1732 + int blossom = _delta2->top();
1.1733 + Node n = _blossom_set->classTop(blossom);
1.1734 + Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
1.1735 + extendOnArc(e);
1.1736 + }
1.1737 + break;
1.1738 + case D3:
1.1739 + {
1.1740 + Edge e = _delta3->top();
1.1741 +
1.1742 + int left_blossom = _blossom_set->find(_graph.u(e));
1.1743 + int right_blossom = _blossom_set->find(_graph.v(e));
1.1744 +
1.1745 + if (left_blossom == right_blossom) {
1.1746 + _delta3->pop();
1.1747 + } else {
1.1748 + int left_tree;
1.1749 + if ((*_blossom_data)[left_blossom].status == EVEN) {
1.1750 + left_tree = _tree_set->find(left_blossom);
1.1751 + } else {
1.1752 + left_tree = -1;
1.1753 + ++unmatched;
1.1754 + }
1.1755 + int right_tree;
1.1756 + if ((*_blossom_data)[right_blossom].status == EVEN) {
1.1757 + right_tree = _tree_set->find(right_blossom);
1.1758 + } else {
1.1759 + right_tree = -1;
1.1760 + ++unmatched;
1.1761 + }
1.1762 +
1.1763 + if (left_tree == right_tree) {
1.1764 + shrinkOnEdge(e, left_tree);
1.1765 + } else {
1.1766 + augmentOnEdge(e);
1.1767 + unmatched -= 2;
1.1768 + }
1.1769 + }
1.1770 + } break;
1.1771 + case D4:
1.1772 + splitBlossom(_delta4->top());
1.1773 + break;
1.1774 + }
1.1775 + }
1.1776 + extractMatching();
1.1777 + }
1.1778 +
1.1779 + /// \brief Runs %MaxWeightedMatching algorithm.
1.1780 + ///
1.1781 + /// This method runs the %MaxWeightedMatching algorithm.
1.1782 + ///
1.1783 + /// \note mwm.run() is just a shortcut of the following code.
1.1784 + /// \code
1.1785 + /// mwm.init();
1.1786 + /// mwm.start();
1.1787 + /// \endcode
1.1788 + void run() {
1.1789 + init();
1.1790 + start();
1.1791 + }
1.1792 +
1.1793 + /// @}
1.1794 +
1.1795 + /// \name Primal solution
1.1796 + /// Functions to get the primal solution, ie. the matching.
1.1797 +
1.1798 + /// @{
1.1799 +
1.1800 + /// \brief Returns the weight of the matching.
1.1801 + ///
1.1802 + /// Returns the weight of the matching.
1.1803 + Value matchingValue() const {
1.1804 + Value sum = 0;
1.1805 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1806 + if ((*_matching)[n] != INVALID) {
1.1807 + sum += _weight[(*_matching)[n]];
1.1808 + }
1.1809 + }
1.1810 + return sum /= 2;
1.1811 + }
1.1812 +
1.1813 + /// \brief Returns the cardinality of the matching.
1.1814 + ///
1.1815 + /// Returns the cardinality of the matching.
1.1816 + int matchingSize() const {
1.1817 + int num = 0;
1.1818 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1819 + if ((*_matching)[n] != INVALID) {
1.1820 + ++num;
1.1821 + }
1.1822 + }
1.1823 + return num /= 2;
1.1824 + }
1.1825 +
1.1826 + /// \brief Returns true when the edge is in the matching.
1.1827 + ///
1.1828 + /// Returns true when the edge is in the matching.
1.1829 + bool matching(const Edge& edge) const {
1.1830 + return edge == (*_matching)[_graph.u(edge)];
1.1831 + }
1.1832 +
1.1833 + /// \brief Returns the incident matching arc.
1.1834 + ///
1.1835 + /// Returns the incident matching arc from given node. If the
1.1836 + /// node is not matched then it gives back \c INVALID.
1.1837 + Arc matching(const Node& node) const {
1.1838 + return (*_matching)[node];
1.1839 + }
1.1840 +
1.1841 + /// \brief Returns the mate of the node.
1.1842 + ///
1.1843 + /// Returns the adjancent node in a mathcing arc. If the node is
1.1844 + /// not matched then it gives back \c INVALID.
1.1845 + Node mate(const Node& node) const {
1.1846 + return (*_matching)[node] != INVALID ?
1.1847 + _graph.target((*_matching)[node]) : INVALID;
1.1848 + }
1.1849 +
1.1850 + /// @}
1.1851 +
1.1852 + /// \name Dual solution
1.1853 + /// Functions to get the dual solution.
1.1854 +
1.1855 + /// @{
1.1856 +
1.1857 + /// \brief Returns the value of the dual solution.
1.1858 + ///
1.1859 + /// Returns the value of the dual solution. It should be equal to
1.1860 + /// the primal value scaled by \ref dualScale "dual scale".
1.1861 + Value dualValue() const {
1.1862 + Value sum = 0;
1.1863 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1864 + sum += nodeValue(n);
1.1865 + }
1.1866 + for (int i = 0; i < blossomNum(); ++i) {
1.1867 + sum += blossomValue(i) * (blossomSize(i) / 2);
1.1868 + }
1.1869 + return sum;
1.1870 + }
1.1871 +
1.1872 + /// \brief Returns the value of the node.
1.1873 + ///
1.1874 + /// Returns the the value of the node.
1.1875 + Value nodeValue(const Node& n) const {
1.1876 + return (*_node_potential)[n];
1.1877 + }
1.1878 +
1.1879 + /// \brief Returns the number of the blossoms in the basis.
1.1880 + ///
1.1881 + /// Returns the number of the blossoms in the basis.
1.1882 + /// \see BlossomIt
1.1883 + int blossomNum() const {
1.1884 + return _blossom_potential.size();
1.1885 + }
1.1886 +
1.1887 +
1.1888 + /// \brief Returns the number of the nodes in the blossom.
1.1889 + ///
1.1890 + /// Returns the number of the nodes in the blossom.
1.1891 + int blossomSize(int k) const {
1.1892 + return _blossom_potential[k].end - _blossom_potential[k].begin;
1.1893 + }
1.1894 +
1.1895 + /// \brief Returns the value of the blossom.
1.1896 + ///
1.1897 + /// Returns the the value of the blossom.
1.1898 + /// \see BlossomIt
1.1899 + Value blossomValue(int k) const {
1.1900 + return _blossom_potential[k].value;
1.1901 + }
1.1902 +
1.1903 + /// \brief Iterator for obtaining the nodes of the blossom.
1.1904 + ///
1.1905 + /// Iterator for obtaining the nodes of the blossom. This class
1.1906 + /// provides a common lemon style iterator for listing a
1.1907 + /// subset of the nodes.
1.1908 + class BlossomIt {
1.1909 + public:
1.1910 +
1.1911 + /// \brief Constructor.
1.1912 + ///
1.1913 + /// Constructor to get the nodes of the variable.
1.1914 + BlossomIt(const MaxWeightedMatching& algorithm, int variable)
1.1915 + : _algorithm(&algorithm)
1.1916 + {
1.1917 + _index = _algorithm->_blossom_potential[variable].begin;
1.1918 + _last = _algorithm->_blossom_potential[variable].end;
1.1919 + }
1.1920 +
1.1921 + /// \brief Conversion to node.
1.1922 + ///
1.1923 + /// Conversion to node.
1.1924 + operator Node() const {
1.1925 + return _algorithm->_blossom_node_list[_index];
1.1926 + }
1.1927 +
1.1928 + /// \brief Increment operator.
1.1929 + ///
1.1930 + /// Increment operator.
1.1931 + BlossomIt& operator++() {
1.1932 + ++_index;
1.1933 + return *this;
1.1934 + }
1.1935 +
1.1936 + /// \brief Validity checking
1.1937 + ///
1.1938 + /// Checks whether the iterator is invalid.
1.1939 + bool operator==(Invalid) const { return _index == _last; }
1.1940 +
1.1941 + /// \brief Validity checking
1.1942 + ///
1.1943 + /// Checks whether the iterator is valid.
1.1944 + bool operator!=(Invalid) const { return _index != _last; }
1.1945 +
1.1946 + private:
1.1947 + const MaxWeightedMatching* _algorithm;
1.1948 + int _last;
1.1949 + int _index;
1.1950 + };
1.1951 +
1.1952 + /// @}
1.1953 +
1.1954 + };
1.1955 +
1.1956 + /// \ingroup matching
1.1957 + ///
1.1958 + /// \brief Weighted perfect matching in general graphs
1.1959 + ///
1.1960 + /// This class provides an efficient implementation of Edmond's
1.1961 + /// maximum weighted perfect matching algorithm. The implementation
1.1962 + /// is based on extensive use of priority queues and provides
1.1963 + /// \f$O(nm\log(n))\f$ time complexity.
1.1964 + ///
1.1965 + /// The maximum weighted matching problem is to find undirected
1.1966 + /// edges in the graph with maximum overall weight and no two of
1.1967 + /// them shares their ends and covers all nodes. The problem can be
1.1968 + /// formulated with the following linear program.
1.1969 + /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1.1970 + /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
1.1971 + \quad \forall B\in\mathcal{O}\f] */
1.1972 + /// \f[x_e \ge 0\quad \forall e\in E\f]
1.1973 + /// \f[\max \sum_{e\in E}x_ew_e\f]
1.1974 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1.1975 + /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
1.1976 + /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
1.1977 + /// subsets of the nodes.
1.1978 + ///
1.1979 + /// The algorithm calculates an optimal matching and a proof of the
1.1980 + /// optimality. The solution of the dual problem can be used to check
1.1981 + /// the result of the algorithm. The dual linear problem is the
1.1982 + /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
1.1983 + w_{uv} \quad \forall uv\in E\f] */
1.1984 + /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
1.1985 + /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
1.1986 + \frac{\vert B \vert - 1}{2}z_B\f] */
1.1987 + ///
1.1988 + /// The algorithm can be executed with \c run() or the \c init() and
1.1989 + /// then the \c start() member functions. After it the matching can
1.1990 + /// be asked with \c matching() or mate() functions. The dual
1.1991 + /// solution can be get with \c nodeValue(), \c blossomNum() and \c
1.1992 + /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
1.1993 + /// "BlossomIt" nested class which is able to iterate on the nodes
1.1994 + /// of a blossom. If the value type is integral then the dual
1.1995 + /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
1.1996 + template <typename _Graph,
1.1997 + typename _WeightMap = typename _Graph::template EdgeMap<int> >
1.1998 + class MaxWeightedPerfectMatching {
1.1999 + public:
1.2000 +
1.2001 + typedef _Graph Graph;
1.2002 + typedef _WeightMap WeightMap;
1.2003 + typedef typename WeightMap::Value Value;
1.2004 +
1.2005 + /// \brief Scaling factor for dual solution
1.2006 + ///
1.2007 + /// Scaling factor for dual solution, it is equal to 4 or 1
1.2008 + /// according to the value type.
1.2009 + static const int dualScale =
1.2010 + std::numeric_limits<Value>::is_integer ? 4 : 1;
1.2011 +
1.2012 + typedef typename Graph::template NodeMap<typename Graph::Arc>
1.2013 + MatchingMap;
1.2014 +
1.2015 + private:
1.2016 +
1.2017 + TEMPLATE_GRAPH_TYPEDEFS(Graph);
1.2018 +
1.2019 + typedef typename Graph::template NodeMap<Value> NodePotential;
1.2020 + typedef std::vector<Node> BlossomNodeList;
1.2021 +
1.2022 + struct BlossomVariable {
1.2023 + int begin, end;
1.2024 + Value value;
1.2025 +
1.2026 + BlossomVariable(int _begin, int _end, Value _value)
1.2027 + : begin(_begin), end(_end), value(_value) {}
1.2028 +
1.2029 + };
1.2030 +
1.2031 + typedef std::vector<BlossomVariable> BlossomPotential;
1.2032 +
1.2033 + const Graph& _graph;
1.2034 + const WeightMap& _weight;
1.2035 +
1.2036 + MatchingMap* _matching;
1.2037 +
1.2038 + NodePotential* _node_potential;
1.2039 +
1.2040 + BlossomPotential _blossom_potential;
1.2041 + BlossomNodeList _blossom_node_list;
1.2042 +
1.2043 + int _node_num;
1.2044 + int _blossom_num;
1.2045 +
1.2046 + typedef RangeMap<int> IntIntMap;
1.2047 +
1.2048 + enum Status {
1.2049 + EVEN = -1, MATCHED = 0, ODD = 1
1.2050 + };
1.2051 +
1.2052 + typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
1.2053 + struct BlossomData {
1.2054 + int tree;
1.2055 + Status status;
1.2056 + Arc pred, next;
1.2057 + Value pot, offset;
1.2058 + };
1.2059 +
1.2060 + IntNodeMap *_blossom_index;
1.2061 + BlossomSet *_blossom_set;
1.2062 + RangeMap<BlossomData>* _blossom_data;
1.2063 +
1.2064 + IntNodeMap *_node_index;
1.2065 + IntArcMap *_node_heap_index;
1.2066 +
1.2067 + struct NodeData {
1.2068 +
1.2069 + NodeData(IntArcMap& node_heap_index)
1.2070 + : heap(node_heap_index) {}
1.2071 +
1.2072 + int blossom;
1.2073 + Value pot;
1.2074 + BinHeap<Value, IntArcMap> heap;
1.2075 + std::map<int, Arc> heap_index;
1.2076 +
1.2077 + int tree;
1.2078 + };
1.2079 +
1.2080 + RangeMap<NodeData>* _node_data;
1.2081 +
1.2082 + typedef ExtendFindEnum<IntIntMap> TreeSet;
1.2083 +
1.2084 + IntIntMap *_tree_set_index;
1.2085 + TreeSet *_tree_set;
1.2086 +
1.2087 + IntIntMap *_delta2_index;
1.2088 + BinHeap<Value, IntIntMap> *_delta2;
1.2089 +
1.2090 + IntEdgeMap *_delta3_index;
1.2091 + BinHeap<Value, IntEdgeMap> *_delta3;
1.2092 +
1.2093 + IntIntMap *_delta4_index;
1.2094 + BinHeap<Value, IntIntMap> *_delta4;
1.2095 +
1.2096 + Value _delta_sum;
1.2097 +
1.2098 + void createStructures() {
1.2099 + _node_num = countNodes(_graph);
1.2100 + _blossom_num = _node_num * 3 / 2;
1.2101 +
1.2102 + if (!_matching) {
1.2103 + _matching = new MatchingMap(_graph);
1.2104 + }
1.2105 + if (!_node_potential) {
1.2106 + _node_potential = new NodePotential(_graph);
1.2107 + }
1.2108 + if (!_blossom_set) {
1.2109 + _blossom_index = new IntNodeMap(_graph);
1.2110 + _blossom_set = new BlossomSet(*_blossom_index);
1.2111 + _blossom_data = new RangeMap<BlossomData>(_blossom_num);
1.2112 + }
1.2113 +
1.2114 + if (!_node_index) {
1.2115 + _node_index = new IntNodeMap(_graph);
1.2116 + _node_heap_index = new IntArcMap(_graph);
1.2117 + _node_data = new RangeMap<NodeData>(_node_num,
1.2118 + NodeData(*_node_heap_index));
1.2119 + }
1.2120 +
1.2121 + if (!_tree_set) {
1.2122 + _tree_set_index = new IntIntMap(_blossom_num);
1.2123 + _tree_set = new TreeSet(*_tree_set_index);
1.2124 + }
1.2125 + if (!_delta2) {
1.2126 + _delta2_index = new IntIntMap(_blossom_num);
1.2127 + _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
1.2128 + }
1.2129 + if (!_delta3) {
1.2130 + _delta3_index = new IntEdgeMap(_graph);
1.2131 + _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
1.2132 + }
1.2133 + if (!_delta4) {
1.2134 + _delta4_index = new IntIntMap(_blossom_num);
1.2135 + _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
1.2136 + }
1.2137 + }
1.2138 +
1.2139 + void destroyStructures() {
1.2140 + _node_num = countNodes(_graph);
1.2141 + _blossom_num = _node_num * 3 / 2;
1.2142 +
1.2143 + if (_matching) {
1.2144 + delete _matching;
1.2145 + }
1.2146 + if (_node_potential) {
1.2147 + delete _node_potential;
1.2148 + }
1.2149 + if (_blossom_set) {
1.2150 + delete _blossom_index;
1.2151 + delete _blossom_set;
1.2152 + delete _blossom_data;
1.2153 + }
1.2154 +
1.2155 + if (_node_index) {
1.2156 + delete _node_index;
1.2157 + delete _node_heap_index;
1.2158 + delete _node_data;
1.2159 + }
1.2160 +
1.2161 + if (_tree_set) {
1.2162 + delete _tree_set_index;
1.2163 + delete _tree_set;
1.2164 + }
1.2165 + if (_delta2) {
1.2166 + delete _delta2_index;
1.2167 + delete _delta2;
1.2168 + }
1.2169 + if (_delta3) {
1.2170 + delete _delta3_index;
1.2171 + delete _delta3;
1.2172 + }
1.2173 + if (_delta4) {
1.2174 + delete _delta4_index;
1.2175 + delete _delta4;
1.2176 + }
1.2177 + }
1.2178 +
1.2179 + void matchedToEven(int blossom, int tree) {
1.2180 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.2181 + _delta2->erase(blossom);
1.2182 + }
1.2183 +
1.2184 + if (!_blossom_set->trivial(blossom)) {
1.2185 + (*_blossom_data)[blossom].pot -=
1.2186 + 2 * (_delta_sum - (*_blossom_data)[blossom].offset);
1.2187 + }
1.2188 +
1.2189 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.2190 + n != INVALID; ++n) {
1.2191 +
1.2192 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
1.2193 + int ni = (*_node_index)[n];
1.2194 +
1.2195 + (*_node_data)[ni].heap.clear();
1.2196 + (*_node_data)[ni].heap_index.clear();
1.2197 +
1.2198 + (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
1.2199 +
1.2200 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.2201 + Node v = _graph.source(e);
1.2202 + int vb = _blossom_set->find(v);
1.2203 + int vi = (*_node_index)[v];
1.2204 +
1.2205 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.2206 + dualScale * _weight[e];
1.2207 +
1.2208 + if ((*_blossom_data)[vb].status == EVEN) {
1.2209 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1.2210 + _delta3->push(e, rw / 2);
1.2211 + }
1.2212 + } else {
1.2213 + typename std::map<int, Arc>::iterator it =
1.2214 + (*_node_data)[vi].heap_index.find(tree);
1.2215 +
1.2216 + if (it != (*_node_data)[vi].heap_index.end()) {
1.2217 + if ((*_node_data)[vi].heap[it->second] > rw) {
1.2218 + (*_node_data)[vi].heap.replace(it->second, e);
1.2219 + (*_node_data)[vi].heap.decrease(e, rw);
1.2220 + it->second = e;
1.2221 + }
1.2222 + } else {
1.2223 + (*_node_data)[vi].heap.push(e, rw);
1.2224 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1.2225 + }
1.2226 +
1.2227 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1.2228 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1.2229 +
1.2230 + if ((*_blossom_data)[vb].status == MATCHED) {
1.2231 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
1.2232 + _delta2->push(vb, _blossom_set->classPrio(vb) -
1.2233 + (*_blossom_data)[vb].offset);
1.2234 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1.2235 + (*_blossom_data)[vb].offset){
1.2236 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1.2237 + (*_blossom_data)[vb].offset);
1.2238 + }
1.2239 + }
1.2240 + }
1.2241 + }
1.2242 + }
1.2243 + }
1.2244 + (*_blossom_data)[blossom].offset = 0;
1.2245 + }
1.2246 +
1.2247 + void matchedToOdd(int blossom) {
1.2248 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.2249 + _delta2->erase(blossom);
1.2250 + }
1.2251 + (*_blossom_data)[blossom].offset += _delta_sum;
1.2252 + if (!_blossom_set->trivial(blossom)) {
1.2253 + _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
1.2254 + (*_blossom_data)[blossom].offset);
1.2255 + }
1.2256 + }
1.2257 +
1.2258 + void evenToMatched(int blossom, int tree) {
1.2259 + if (!_blossom_set->trivial(blossom)) {
1.2260 + (*_blossom_data)[blossom].pot += 2 * _delta_sum;
1.2261 + }
1.2262 +
1.2263 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.2264 + n != INVALID; ++n) {
1.2265 + int ni = (*_node_index)[n];
1.2266 + (*_node_data)[ni].pot -= _delta_sum;
1.2267 +
1.2268 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.2269 + Node v = _graph.source(e);
1.2270 + int vb = _blossom_set->find(v);
1.2271 + int vi = (*_node_index)[v];
1.2272 +
1.2273 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.2274 + dualScale * _weight[e];
1.2275 +
1.2276 + if (vb == blossom) {
1.2277 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.2278 + _delta3->erase(e);
1.2279 + }
1.2280 + } else if ((*_blossom_data)[vb].status == EVEN) {
1.2281 +
1.2282 + if (_delta3->state(e) == _delta3->IN_HEAP) {
1.2283 + _delta3->erase(e);
1.2284 + }
1.2285 +
1.2286 + int vt = _tree_set->find(vb);
1.2287 +
1.2288 + if (vt != tree) {
1.2289 +
1.2290 + Arc r = _graph.oppositeArc(e);
1.2291 +
1.2292 + typename std::map<int, Arc>::iterator it =
1.2293 + (*_node_data)[ni].heap_index.find(vt);
1.2294 +
1.2295 + if (it != (*_node_data)[ni].heap_index.end()) {
1.2296 + if ((*_node_data)[ni].heap[it->second] > rw) {
1.2297 + (*_node_data)[ni].heap.replace(it->second, r);
1.2298 + (*_node_data)[ni].heap.decrease(r, rw);
1.2299 + it->second = r;
1.2300 + }
1.2301 + } else {
1.2302 + (*_node_data)[ni].heap.push(r, rw);
1.2303 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1.2304 + }
1.2305 +
1.2306 + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1.2307 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1.2308 +
1.2309 + if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1.2310 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.2311 + (*_blossom_data)[blossom].offset);
1.2312 + } else if ((*_delta2)[blossom] >
1.2313 + _blossom_set->classPrio(blossom) -
1.2314 + (*_blossom_data)[blossom].offset){
1.2315 + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1.2316 + (*_blossom_data)[blossom].offset);
1.2317 + }
1.2318 + }
1.2319 + }
1.2320 + } else {
1.2321 +
1.2322 + typename std::map<int, Arc>::iterator it =
1.2323 + (*_node_data)[vi].heap_index.find(tree);
1.2324 +
1.2325 + if (it != (*_node_data)[vi].heap_index.end()) {
1.2326 + (*_node_data)[vi].heap.erase(it->second);
1.2327 + (*_node_data)[vi].heap_index.erase(it);
1.2328 + if ((*_node_data)[vi].heap.empty()) {
1.2329 + _blossom_set->increase(v, std::numeric_limits<Value>::max());
1.2330 + } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
1.2331 + _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
1.2332 + }
1.2333 +
1.2334 + if ((*_blossom_data)[vb].status == MATCHED) {
1.2335 + if (_blossom_set->classPrio(vb) ==
1.2336 + std::numeric_limits<Value>::max()) {
1.2337 + _delta2->erase(vb);
1.2338 + } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
1.2339 + (*_blossom_data)[vb].offset) {
1.2340 + _delta2->increase(vb, _blossom_set->classPrio(vb) -
1.2341 + (*_blossom_data)[vb].offset);
1.2342 + }
1.2343 + }
1.2344 + }
1.2345 + }
1.2346 + }
1.2347 + }
1.2348 + }
1.2349 +
1.2350 + void oddToMatched(int blossom) {
1.2351 + (*_blossom_data)[blossom].offset -= _delta_sum;
1.2352 +
1.2353 + if (_blossom_set->classPrio(blossom) !=
1.2354 + std::numeric_limits<Value>::max()) {
1.2355 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1.2356 + (*_blossom_data)[blossom].offset);
1.2357 + }
1.2358 +
1.2359 + if (!_blossom_set->trivial(blossom)) {
1.2360 + _delta4->erase(blossom);
1.2361 + }
1.2362 + }
1.2363 +
1.2364 + void oddToEven(int blossom, int tree) {
1.2365 + if (!_blossom_set->trivial(blossom)) {
1.2366 + _delta4->erase(blossom);
1.2367 + (*_blossom_data)[blossom].pot -=
1.2368 + 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
1.2369 + }
1.2370 +
1.2371 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1.2372 + n != INVALID; ++n) {
1.2373 + int ni = (*_node_index)[n];
1.2374 +
1.2375 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
1.2376 +
1.2377 + (*_node_data)[ni].heap.clear();
1.2378 + (*_node_data)[ni].heap_index.clear();
1.2379 + (*_node_data)[ni].pot +=
1.2380 + 2 * _delta_sum - (*_blossom_data)[blossom].offset;
1.2381 +
1.2382 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.2383 + Node v = _graph.source(e);
1.2384 + int vb = _blossom_set->find(v);
1.2385 + int vi = (*_node_index)[v];
1.2386 +
1.2387 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1.2388 + dualScale * _weight[e];
1.2389 +
1.2390 + if ((*_blossom_data)[vb].status == EVEN) {
1.2391 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1.2392 + _delta3->push(e, rw / 2);
1.2393 + }
1.2394 + } else {
1.2395 +
1.2396 + typename std::map<int, Arc>::iterator it =
1.2397 + (*_node_data)[vi].heap_index.find(tree);
1.2398 +
1.2399 + if (it != (*_node_data)[vi].heap_index.end()) {
1.2400 + if ((*_node_data)[vi].heap[it->second] > rw) {
1.2401 + (*_node_data)[vi].heap.replace(it->second, e);
1.2402 + (*_node_data)[vi].heap.decrease(e, rw);
1.2403 + it->second = e;
1.2404 + }
1.2405 + } else {
1.2406 + (*_node_data)[vi].heap.push(e, rw);
1.2407 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
1.2408 + }
1.2409 +
1.2410 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
1.2411 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
1.2412 +
1.2413 + if ((*_blossom_data)[vb].status == MATCHED) {
1.2414 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
1.2415 + _delta2->push(vb, _blossom_set->classPrio(vb) -
1.2416 + (*_blossom_data)[vb].offset);
1.2417 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
1.2418 + (*_blossom_data)[vb].offset) {
1.2419 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
1.2420 + (*_blossom_data)[vb].offset);
1.2421 + }
1.2422 + }
1.2423 + }
1.2424 + }
1.2425 + }
1.2426 + }
1.2427 + (*_blossom_data)[blossom].offset = 0;
1.2428 + }
1.2429 +
1.2430 + void alternatePath(int even, int tree) {
1.2431 + int odd;
1.2432 +
1.2433 + evenToMatched(even, tree);
1.2434 + (*_blossom_data)[even].status = MATCHED;
1.2435 +
1.2436 + while ((*_blossom_data)[even].pred != INVALID) {
1.2437 + odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
1.2438 + (*_blossom_data)[odd].status = MATCHED;
1.2439 + oddToMatched(odd);
1.2440 + (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
1.2441 +
1.2442 + even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
1.2443 + (*_blossom_data)[even].status = MATCHED;
1.2444 + evenToMatched(even, tree);
1.2445 + (*_blossom_data)[even].next =
1.2446 + _graph.oppositeArc((*_blossom_data)[odd].pred);
1.2447 + }
1.2448 +
1.2449 + }
1.2450 +
1.2451 + void destroyTree(int tree) {
1.2452 + for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
1.2453 + if ((*_blossom_data)[b].status == EVEN) {
1.2454 + (*_blossom_data)[b].status = MATCHED;
1.2455 + evenToMatched(b, tree);
1.2456 + } else if ((*_blossom_data)[b].status == ODD) {
1.2457 + (*_blossom_data)[b].status = MATCHED;
1.2458 + oddToMatched(b);
1.2459 + }
1.2460 + }
1.2461 + _tree_set->eraseClass(tree);
1.2462 + }
1.2463 +
1.2464 + void augmentOnEdge(const Edge& edge) {
1.2465 +
1.2466 + int left = _blossom_set->find(_graph.u(edge));
1.2467 + int right = _blossom_set->find(_graph.v(edge));
1.2468 +
1.2469 + int left_tree = _tree_set->find(left);
1.2470 + alternatePath(left, left_tree);
1.2471 + destroyTree(left_tree);
1.2472 +
1.2473 + int right_tree = _tree_set->find(right);
1.2474 + alternatePath(right, right_tree);
1.2475 + destroyTree(right_tree);
1.2476 +
1.2477 + (*_blossom_data)[left].next = _graph.direct(edge, true);
1.2478 + (*_blossom_data)[right].next = _graph.direct(edge, false);
1.2479 + }
1.2480 +
1.2481 + void extendOnArc(const Arc& arc) {
1.2482 + int base = _blossom_set->find(_graph.target(arc));
1.2483 + int tree = _tree_set->find(base);
1.2484 +
1.2485 + int odd = _blossom_set->find(_graph.source(arc));
1.2486 + _tree_set->insert(odd, tree);
1.2487 + (*_blossom_data)[odd].status = ODD;
1.2488 + matchedToOdd(odd);
1.2489 + (*_blossom_data)[odd].pred = arc;
1.2490 +
1.2491 + int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
1.2492 + (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
1.2493 + _tree_set->insert(even, tree);
1.2494 + (*_blossom_data)[even].status = EVEN;
1.2495 + matchedToEven(even, tree);
1.2496 + }
1.2497 +
1.2498 + void shrinkOnEdge(const Edge& edge, int tree) {
1.2499 + int nca = -1;
1.2500 + std::vector<int> left_path, right_path;
1.2501 +
1.2502 + {
1.2503 + std::set<int> left_set, right_set;
1.2504 + int left = _blossom_set->find(_graph.u(edge));
1.2505 + left_path.push_back(left);
1.2506 + left_set.insert(left);
1.2507 +
1.2508 + int right = _blossom_set->find(_graph.v(edge));
1.2509 + right_path.push_back(right);
1.2510 + right_set.insert(right);
1.2511 +
1.2512 + while (true) {
1.2513 +
1.2514 + if ((*_blossom_data)[left].pred == INVALID) break;
1.2515 +
1.2516 + left =
1.2517 + _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1.2518 + left_path.push_back(left);
1.2519 + left =
1.2520 + _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
1.2521 + left_path.push_back(left);
1.2522 +
1.2523 + left_set.insert(left);
1.2524 +
1.2525 + if (right_set.find(left) != right_set.end()) {
1.2526 + nca = left;
1.2527 + break;
1.2528 + }
1.2529 +
1.2530 + if ((*_blossom_data)[right].pred == INVALID) break;
1.2531 +
1.2532 + right =
1.2533 + _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1.2534 + right_path.push_back(right);
1.2535 + right =
1.2536 + _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
1.2537 + right_path.push_back(right);
1.2538 +
1.2539 + right_set.insert(right);
1.2540 +
1.2541 + if (left_set.find(right) != left_set.end()) {
1.2542 + nca = right;
1.2543 + break;
1.2544 + }
1.2545 +
1.2546 + }
1.2547 +
1.2548 + if (nca == -1) {
1.2549 + if ((*_blossom_data)[left].pred == INVALID) {
1.2550 + nca = right;
1.2551 + while (left_set.find(nca) == left_set.end()) {
1.2552 + nca =
1.2553 + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1.2554 + right_path.push_back(nca);
1.2555 + nca =
1.2556 + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1.2557 + right_path.push_back(nca);
1.2558 + }
1.2559 + } else {
1.2560 + nca = left;
1.2561 + while (right_set.find(nca) == right_set.end()) {
1.2562 + nca =
1.2563 + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1.2564 + left_path.push_back(nca);
1.2565 + nca =
1.2566 + _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
1.2567 + left_path.push_back(nca);
1.2568 + }
1.2569 + }
1.2570 + }
1.2571 + }
1.2572 +
1.2573 + std::vector<int> subblossoms;
1.2574 + Arc prev;
1.2575 +
1.2576 + prev = _graph.direct(edge, true);
1.2577 + for (int i = 0; left_path[i] != nca; i += 2) {
1.2578 + subblossoms.push_back(left_path[i]);
1.2579 + (*_blossom_data)[left_path[i]].next = prev;
1.2580 + _tree_set->erase(left_path[i]);
1.2581 +
1.2582 + subblossoms.push_back(left_path[i + 1]);
1.2583 + (*_blossom_data)[left_path[i + 1]].status = EVEN;
1.2584 + oddToEven(left_path[i + 1], tree);
1.2585 + _tree_set->erase(left_path[i + 1]);
1.2586 + prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
1.2587 + }
1.2588 +
1.2589 + int k = 0;
1.2590 + while (right_path[k] != nca) ++k;
1.2591 +
1.2592 + subblossoms.push_back(nca);
1.2593 + (*_blossom_data)[nca].next = prev;
1.2594 +
1.2595 + for (int i = k - 2; i >= 0; i -= 2) {
1.2596 + subblossoms.push_back(right_path[i + 1]);
1.2597 + (*_blossom_data)[right_path[i + 1]].status = EVEN;
1.2598 + oddToEven(right_path[i + 1], tree);
1.2599 + _tree_set->erase(right_path[i + 1]);
1.2600 +
1.2601 + (*_blossom_data)[right_path[i + 1]].next =
1.2602 + (*_blossom_data)[right_path[i + 1]].pred;
1.2603 +
1.2604 + subblossoms.push_back(right_path[i]);
1.2605 + _tree_set->erase(right_path[i]);
1.2606 + }
1.2607 +
1.2608 + int surface =
1.2609 + _blossom_set->join(subblossoms.begin(), subblossoms.end());
1.2610 +
1.2611 + for (int i = 0; i < int(subblossoms.size()); ++i) {
1.2612 + if (!_blossom_set->trivial(subblossoms[i])) {
1.2613 + (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
1.2614 + }
1.2615 + (*_blossom_data)[subblossoms[i]].status = MATCHED;
1.2616 + }
1.2617 +
1.2618 + (*_blossom_data)[surface].pot = -2 * _delta_sum;
1.2619 + (*_blossom_data)[surface].offset = 0;
1.2620 + (*_blossom_data)[surface].status = EVEN;
1.2621 + (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
1.2622 + (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
1.2623 +
1.2624 + _tree_set->insert(surface, tree);
1.2625 + _tree_set->erase(nca);
1.2626 + }
1.2627 +
1.2628 + void splitBlossom(int blossom) {
1.2629 + Arc next = (*_blossom_data)[blossom].next;
1.2630 + Arc pred = (*_blossom_data)[blossom].pred;
1.2631 +
1.2632 + int tree = _tree_set->find(blossom);
1.2633 +
1.2634 + (*_blossom_data)[blossom].status = MATCHED;
1.2635 + oddToMatched(blossom);
1.2636 + if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1.2637 + _delta2->erase(blossom);
1.2638 + }
1.2639 +
1.2640 + std::vector<int> subblossoms;
1.2641 + _blossom_set->split(blossom, std::back_inserter(subblossoms));
1.2642 +
1.2643 + Value offset = (*_blossom_data)[blossom].offset;
1.2644 + int b = _blossom_set->find(_graph.source(pred));
1.2645 + int d = _blossom_set->find(_graph.source(next));
1.2646 +
1.2647 + int ib = -1, id = -1;
1.2648 + for (int i = 0; i < int(subblossoms.size()); ++i) {
1.2649 + if (subblossoms[i] == b) ib = i;
1.2650 + if (subblossoms[i] == d) id = i;
1.2651 +
1.2652 + (*_blossom_data)[subblossoms[i]].offset = offset;
1.2653 + if (!_blossom_set->trivial(subblossoms[i])) {
1.2654 + (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
1.2655 + }
1.2656 + if (_blossom_set->classPrio(subblossoms[i]) !=
1.2657 + std::numeric_limits<Value>::max()) {
1.2658 + _delta2->push(subblossoms[i],
1.2659 + _blossom_set->classPrio(subblossoms[i]) -
1.2660 + (*_blossom_data)[subblossoms[i]].offset);
1.2661 + }
1.2662 + }
1.2663 +
1.2664 + if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
1.2665 + for (int i = (id + 1) % subblossoms.size();
1.2666 + i != ib; i = (i + 2) % subblossoms.size()) {
1.2667 + int sb = subblossoms[i];
1.2668 + int tb = subblossoms[(i + 1) % subblossoms.size()];
1.2669 + (*_blossom_data)[sb].next =
1.2670 + _graph.oppositeArc((*_blossom_data)[tb].next);
1.2671 + }
1.2672 +
1.2673 + for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
1.2674 + int sb = subblossoms[i];
1.2675 + int tb = subblossoms[(i + 1) % subblossoms.size()];
1.2676 + int ub = subblossoms[(i + 2) % subblossoms.size()];
1.2677 +
1.2678 + (*_blossom_data)[sb].status = ODD;
1.2679 + matchedToOdd(sb);
1.2680 + _tree_set->insert(sb, tree);
1.2681 + (*_blossom_data)[sb].pred = pred;
1.2682 + (*_blossom_data)[sb].next =
1.2683 + _graph.oppositeArc((*_blossom_data)[tb].next);
1.2684 +
1.2685 + pred = (*_blossom_data)[ub].next;
1.2686 +
1.2687 + (*_blossom_data)[tb].status = EVEN;
1.2688 + matchedToEven(tb, tree);
1.2689 + _tree_set->insert(tb, tree);
1.2690 + (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
1.2691 + }
1.2692 +
1.2693 + (*_blossom_data)[subblossoms[id]].status = ODD;
1.2694 + matchedToOdd(subblossoms[id]);
1.2695 + _tree_set->insert(subblossoms[id], tree);
1.2696 + (*_blossom_data)[subblossoms[id]].next = next;
1.2697 + (*_blossom_data)[subblossoms[id]].pred = pred;
1.2698 +
1.2699 + } else {
1.2700 +
1.2701 + for (int i = (ib + 1) % subblossoms.size();
1.2702 + i != id; i = (i + 2) % subblossoms.size()) {
1.2703 + int sb = subblossoms[i];
1.2704 + int tb = subblossoms[(i + 1) % subblossoms.size()];
1.2705 + (*_blossom_data)[sb].next =
1.2706 + _graph.oppositeArc((*_blossom_data)[tb].next);
1.2707 + }
1.2708 +
1.2709 + for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
1.2710 + int sb = subblossoms[i];
1.2711 + int tb = subblossoms[(i + 1) % subblossoms.size()];
1.2712 + int ub = subblossoms[(i + 2) % subblossoms.size()];
1.2713 +
1.2714 + (*_blossom_data)[sb].status = ODD;
1.2715 + matchedToOdd(sb);
1.2716 + _tree_set->insert(sb, tree);
1.2717 + (*_blossom_data)[sb].next = next;
1.2718 + (*_blossom_data)[sb].pred =
1.2719 + _graph.oppositeArc((*_blossom_data)[tb].next);
1.2720 +
1.2721 + (*_blossom_data)[tb].status = EVEN;
1.2722 + matchedToEven(tb, tree);
1.2723 + _tree_set->insert(tb, tree);
1.2724 + (*_blossom_data)[tb].pred =
1.2725 + (*_blossom_data)[tb].next =
1.2726 + _graph.oppositeArc((*_blossom_data)[ub].next);
1.2727 + next = (*_blossom_data)[ub].next;
1.2728 + }
1.2729 +
1.2730 + (*_blossom_data)[subblossoms[ib]].status = ODD;
1.2731 + matchedToOdd(subblossoms[ib]);
1.2732 + _tree_set->insert(subblossoms[ib], tree);
1.2733 + (*_blossom_data)[subblossoms[ib]].next = next;
1.2734 + (*_blossom_data)[subblossoms[ib]].pred = pred;
1.2735 + }
1.2736 + _tree_set->erase(blossom);
1.2737 + }
1.2738 +
1.2739 + void extractBlossom(int blossom, const Node& base, const Arc& matching) {
1.2740 + if (_blossom_set->trivial(blossom)) {
1.2741 + int bi = (*_node_index)[base];
1.2742 + Value pot = (*_node_data)[bi].pot;
1.2743 +
1.2744 + _matching->set(base, matching);
1.2745 + _blossom_node_list.push_back(base);
1.2746 + _node_potential->set(base, pot);
1.2747 + } else {
1.2748 +
1.2749 + Value pot = (*_blossom_data)[blossom].pot;
1.2750 + int bn = _blossom_node_list.size();
1.2751 +
1.2752 + std::vector<int> subblossoms;
1.2753 + _blossom_set->split(blossom, std::back_inserter(subblossoms));
1.2754 + int b = _blossom_set->find(base);
1.2755 + int ib = -1;
1.2756 + for (int i = 0; i < int(subblossoms.size()); ++i) {
1.2757 + if (subblossoms[i] == b) { ib = i; break; }
1.2758 + }
1.2759 +
1.2760 + for (int i = 1; i < int(subblossoms.size()); i += 2) {
1.2761 + int sb = subblossoms[(ib + i) % subblossoms.size()];
1.2762 + int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
1.2763 +
1.2764 + Arc m = (*_blossom_data)[tb].next;
1.2765 + extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
1.2766 + extractBlossom(tb, _graph.source(m), m);
1.2767 + }
1.2768 + extractBlossom(subblossoms[ib], base, matching);
1.2769 +
1.2770 + int en = _blossom_node_list.size();
1.2771 +
1.2772 + _blossom_potential.push_back(BlossomVariable(bn, en, pot));
1.2773 + }
1.2774 + }
1.2775 +
1.2776 + void extractMatching() {
1.2777 + std::vector<int> blossoms;
1.2778 + for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1.2779 + blossoms.push_back(c);
1.2780 + }
1.2781 +
1.2782 + for (int i = 0; i < int(blossoms.size()); ++i) {
1.2783 +
1.2784 + Value offset = (*_blossom_data)[blossoms[i]].offset;
1.2785 + (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1.2786 + for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1.2787 + n != INVALID; ++n) {
1.2788 + (*_node_data)[(*_node_index)[n]].pot -= offset;
1.2789 + }
1.2790 +
1.2791 + Arc matching = (*_blossom_data)[blossoms[i]].next;
1.2792 + Node base = _graph.source(matching);
1.2793 + extractBlossom(blossoms[i], base, matching);
1.2794 + }
1.2795 + }
1.2796 +
1.2797 + public:
1.2798 +
1.2799 + /// \brief Constructor
1.2800 + ///
1.2801 + /// Constructor.
1.2802 + MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
1.2803 + : _graph(graph), _weight(weight), _matching(0),
1.2804 + _node_potential(0), _blossom_potential(), _blossom_node_list(),
1.2805 + _node_num(0), _blossom_num(0),
1.2806 +
1.2807 + _blossom_index(0), _blossom_set(0), _blossom_data(0),
1.2808 + _node_index(0), _node_heap_index(0), _node_data(0),
1.2809 + _tree_set_index(0), _tree_set(0),
1.2810 +
1.2811 + _delta2_index(0), _delta2(0),
1.2812 + _delta3_index(0), _delta3(0),
1.2813 + _delta4_index(0), _delta4(0),
1.2814 +
1.2815 + _delta_sum() {}
1.2816 +
1.2817 + ~MaxWeightedPerfectMatching() {
1.2818 + destroyStructures();
1.2819 + }
1.2820 +
1.2821 + /// \name Execution control
1.2822 + /// The simplest way to execute the algorithm is to use the
1.2823 + /// \c run() member function.
1.2824 +
1.2825 + ///@{
1.2826 +
1.2827 + /// \brief Initialize the algorithm
1.2828 + ///
1.2829 + /// Initialize the algorithm
1.2830 + void init() {
1.2831 + createStructures();
1.2832 +
1.2833 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.2834 + _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
1.2835 + }
1.2836 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.2837 + _delta3_index->set(e, _delta3->PRE_HEAP);
1.2838 + }
1.2839 + for (int i = 0; i < _blossom_num; ++i) {
1.2840 + _delta2_index->set(i, _delta2->PRE_HEAP);
1.2841 + _delta4_index->set(i, _delta4->PRE_HEAP);
1.2842 + }
1.2843 +
1.2844 + int index = 0;
1.2845 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.2846 + Value max = - std::numeric_limits<Value>::max();
1.2847 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.2848 + if (_graph.target(e) == n) continue;
1.2849 + if ((dualScale * _weight[e]) / 2 > max) {
1.2850 + max = (dualScale * _weight[e]) / 2;
1.2851 + }
1.2852 + }
1.2853 + _node_index->set(n, index);
1.2854 + (*_node_data)[index].pot = max;
1.2855 + int blossom =
1.2856 + _blossom_set->insert(n, std::numeric_limits<Value>::max());
1.2857 +
1.2858 + _tree_set->insert(blossom);
1.2859 +
1.2860 + (*_blossom_data)[blossom].status = EVEN;
1.2861 + (*_blossom_data)[blossom].pred = INVALID;
1.2862 + (*_blossom_data)[blossom].next = INVALID;
1.2863 + (*_blossom_data)[blossom].pot = 0;
1.2864 + (*_blossom_data)[blossom].offset = 0;
1.2865 + ++index;
1.2866 + }
1.2867 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.2868 + int si = (*_node_index)[_graph.u(e)];
1.2869 + int ti = (*_node_index)[_graph.v(e)];
1.2870 + if (_graph.u(e) != _graph.v(e)) {
1.2871 + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
1.2872 + dualScale * _weight[e]) / 2);
1.2873 + }
1.2874 + }
1.2875 + }
1.2876 +
1.2877 + /// \brief Starts the algorithm
1.2878 + ///
1.2879 + /// Starts the algorithm
1.2880 + bool start() {
1.2881 + enum OpType {
1.2882 + D2, D3, D4
1.2883 + };
1.2884 +
1.2885 + int unmatched = _node_num;
1.2886 + while (unmatched > 0) {
1.2887 + Value d2 = !_delta2->empty() ?
1.2888 + _delta2->prio() : std::numeric_limits<Value>::max();
1.2889 +
1.2890 + Value d3 = !_delta3->empty() ?
1.2891 + _delta3->prio() : std::numeric_limits<Value>::max();
1.2892 +
1.2893 + Value d4 = !_delta4->empty() ?
1.2894 + _delta4->prio() : std::numeric_limits<Value>::max();
1.2895 +
1.2896 + _delta_sum = d2; OpType ot = D2;
1.2897 + if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
1.2898 + if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1.2899 +
1.2900 + if (_delta_sum == std::numeric_limits<Value>::max()) {
1.2901 + return false;
1.2902 + }
1.2903 +
1.2904 + switch (ot) {
1.2905 + case D2:
1.2906 + {
1.2907 + int blossom = _delta2->top();
1.2908 + Node n = _blossom_set->classTop(blossom);
1.2909 + Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
1.2910 + extendOnArc(e);
1.2911 + }
1.2912 + break;
1.2913 + case D3:
1.2914 + {
1.2915 + Edge e = _delta3->top();
1.2916 +
1.2917 + int left_blossom = _blossom_set->find(_graph.u(e));
1.2918 + int right_blossom = _blossom_set->find(_graph.v(e));
1.2919 +
1.2920 + if (left_blossom == right_blossom) {
1.2921 + _delta3->pop();
1.2922 + } else {
1.2923 + int left_tree = _tree_set->find(left_blossom);
1.2924 + int right_tree = _tree_set->find(right_blossom);
1.2925 +
1.2926 + if (left_tree == right_tree) {
1.2927 + shrinkOnEdge(e, left_tree);
1.2928 + } else {
1.2929 + augmentOnEdge(e);
1.2930 + unmatched -= 2;
1.2931 + }
1.2932 + }
1.2933 + } break;
1.2934 + case D4:
1.2935 + splitBlossom(_delta4->top());
1.2936 + break;
1.2937 + }
1.2938 + }
1.2939 + extractMatching();
1.2940 + return true;
1.2941 + }
1.2942 +
1.2943 + /// \brief Runs %MaxWeightedPerfectMatching algorithm.
1.2944 + ///
1.2945 + /// This method runs the %MaxWeightedPerfectMatching algorithm.
1.2946 + ///
1.2947 + /// \note mwm.run() is just a shortcut of the following code.
1.2948 + /// \code
1.2949 + /// mwm.init();
1.2950 + /// mwm.start();
1.2951 + /// \endcode
1.2952 + bool run() {
1.2953 + init();
1.2954 + return start();
1.2955 + }
1.2956 +
1.2957 + /// @}
1.2958 +
1.2959 + /// \name Primal solution
1.2960 + /// Functions to get the primal solution, ie. the matching.
1.2961 +
1.2962 + /// @{
1.2963 +
1.2964 + /// \brief Returns the matching value.
1.2965 + ///
1.2966 + /// Returns the matching value.
1.2967 + Value matchingValue() const {
1.2968 + Value sum = 0;
1.2969 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.2970 + if ((*_matching)[n] != INVALID) {
1.2971 + sum += _weight[(*_matching)[n]];
1.2972 + }
1.2973 + }
1.2974 + return sum /= 2;
1.2975 + }
1.2976 +
1.2977 + /// \brief Returns true when the edge is in the matching.
1.2978 + ///
1.2979 + /// Returns true when the edge is in the matching.
1.2980 + bool matching(const Edge& edge) const {
1.2981 + return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
1.2982 + }
1.2983 +
1.2984 + /// \brief Returns the incident matching edge.
1.2985 + ///
1.2986 + /// Returns the incident matching arc from given edge.
1.2987 + Arc matching(const Node& node) const {
1.2988 + return (*_matching)[node];
1.2989 + }
1.2990 +
1.2991 + /// \brief Returns the mate of the node.
1.2992 + ///
1.2993 + /// Returns the adjancent node in a mathcing arc.
1.2994 + Node mate(const Node& node) const {
1.2995 + return _graph.target((*_matching)[node]);
1.2996 + }
1.2997 +
1.2998 + /// @}
1.2999 +
1.3000 + /// \name Dual solution
1.3001 + /// Functions to get the dual solution.
1.3002 +
1.3003 + /// @{
1.3004 +
1.3005 + /// \brief Returns the value of the dual solution.
1.3006 + ///
1.3007 + /// Returns the value of the dual solution. It should be equal to
1.3008 + /// the primal value scaled by \ref dualScale "dual scale".
1.3009 + Value dualValue() const {
1.3010 + Value sum = 0;
1.3011 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.3012 + sum += nodeValue(n);
1.3013 + }
1.3014 + for (int i = 0; i < blossomNum(); ++i) {
1.3015 + sum += blossomValue(i) * (blossomSize(i) / 2);
1.3016 + }
1.3017 + return sum;
1.3018 + }
1.3019 +
1.3020 + /// \brief Returns the value of the node.
1.3021 + ///
1.3022 + /// Returns the the value of the node.
1.3023 + Value nodeValue(const Node& n) const {
1.3024 + return (*_node_potential)[n];
1.3025 + }
1.3026 +
1.3027 + /// \brief Returns the number of the blossoms in the basis.
1.3028 + ///
1.3029 + /// Returns the number of the blossoms in the basis.
1.3030 + /// \see BlossomIt
1.3031 + int blossomNum() const {
1.3032 + return _blossom_potential.size();
1.3033 + }
1.3034 +
1.3035 +
1.3036 + /// \brief Returns the number of the nodes in the blossom.
1.3037 + ///
1.3038 + /// Returns the number of the nodes in the blossom.
1.3039 + int blossomSize(int k) const {
1.3040 + return _blossom_potential[k].end - _blossom_potential[k].begin;
1.3041 + }
1.3042 +
1.3043 + /// \brief Returns the value of the blossom.
1.3044 + ///
1.3045 + /// Returns the the value of the blossom.
1.3046 + /// \see BlossomIt
1.3047 + Value blossomValue(int k) const {
1.3048 + return _blossom_potential[k].value;
1.3049 + }
1.3050 +
1.3051 + /// \brief Iterator for obtaining the nodes of the blossom.
1.3052 + ///
1.3053 + /// Iterator for obtaining the nodes of the blossom. This class
1.3054 + /// provides a common lemon style iterator for listing a
1.3055 + /// subset of the nodes.
1.3056 + class BlossomIt {
1.3057 + public:
1.3058 +
1.3059 + /// \brief Constructor.
1.3060 + ///
1.3061 + /// Constructor to get the nodes of the variable.
1.3062 + BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
1.3063 + : _algorithm(&algorithm)
1.3064 + {
1.3065 + _index = _algorithm->_blossom_potential[variable].begin;
1.3066 + _last = _algorithm->_blossom_potential[variable].end;
1.3067 + }
1.3068 +
1.3069 + /// \brief Conversion to node.
1.3070 + ///
1.3071 + /// Conversion to node.
1.3072 + operator Node() const {
1.3073 + return _algorithm->_blossom_node_list[_index];
1.3074 + }
1.3075 +
1.3076 + /// \brief Increment operator.
1.3077 + ///
1.3078 + /// Increment operator.
1.3079 + BlossomIt& operator++() {
1.3080 + ++_index;
1.3081 + return *this;
1.3082 + }
1.3083 +
1.3084 + /// \brief Validity checking
1.3085 + ///
1.3086 + /// Checks whether the iterator is invalid.
1.3087 + bool operator==(Invalid) const { return _index == _last; }
1.3088 +
1.3089 + /// \brief Validity checking
1.3090 + ///
1.3091 + /// Checks whether the iterator is valid.
1.3092 + bool operator!=(Invalid) const { return _index != _last; }
1.3093 +
1.3094 + private:
1.3095 + const MaxWeightedPerfectMatching* _algorithm;
1.3096 + int _last;
1.3097 + int _index;
1.3098 + };
1.3099 +
1.3100 + /// @}
1.3101 +
1.3102 + };
1.3103 +
1.3104 +
1.3105 +} //END OF NAMESPACE LEMON
1.3106 +
1.3107 +#endif //LEMON_MAX_MATCHING_H