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