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