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