1.1 --- a/lemon/max_matching.h Fri Apr 17 09:54:14 2009 +0200
1.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
1.3 @@ -1,3244 +0,0 @@
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