1.1 --- a/lemon/max_matching.h Mon Oct 13 13:56:00 2008 +0200
1.2 +++ b/lemon/max_matching.h Mon Oct 13 14:00:11 2008 +0200
1.3 @@ -31,76 +31,405 @@
1.4
1.5 ///\ingroup matching
1.6 ///\file
1.7 -///\brief Maximum matching algorithms in graph.
1.8 +///\brief Maximum matching algorithms in general graphs.
1.9
1.10 namespace lemon {
1.11
1.12 - ///\ingroup matching
1.13 + /// \ingroup matching
1.14 ///
1.15 - ///\brief Edmonds' alternating forest maximum matching algorithm.
1.16 + /// \brief Edmonds' alternating forest maximum matching algorithm.
1.17 ///
1.18 - ///This class provides Edmonds' alternating forest matching
1.19 - ///algorithm. The starting matching (if any) can be passed to the
1.20 - ///algorithm using some of init functions.
1.21 + /// This class provides Edmonds' alternating forest matching
1.22 + /// algorithm. The starting matching (if any) can be passed to the
1.23 + /// algorithm using some of init functions.
1.24 ///
1.25 - ///The dual side of a matching is a map of the nodes to
1.26 - ///MaxMatching::DecompType, having values \c D, \c A and \c C
1.27 - ///showing the Gallai-Edmonds decomposition of the digraph. The nodes
1.28 - ///in \c D induce a digraph with factor-critical components, the nodes
1.29 - ///in \c A form the barrier, and the nodes in \c C induce a digraph
1.30 - ///having a perfect matching. This decomposition can be attained by
1.31 - ///calling \c decomposition() after running the algorithm.
1.32 + /// The dual side of a matching is a map of the nodes to
1.33 + /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c
1.34 + /// MATCHED/C showing the Gallai-Edmonds decomposition of the
1.35 + /// graph. The nodes in \c EVEN/D induce a graph with
1.36 + /// factor-critical components, the nodes in \c ODD/A form the
1.37 + /// barrier, and the nodes in \c MATCHED/C induce a graph having a
1.38 + /// perfect matching. The number of the fractor critical components
1.39 + /// minus the number of barrier nodes is a lower bound on the
1.40 + /// unmatched nodes, and if the matching is optimal this bound is
1.41 + /// tight. This decomposition can be attained by calling \c
1.42 + /// decomposition() after running the algorithm.
1.43 ///
1.44 - ///\param Digraph The graph type the algorithm runs on.
1.45 - template <typename Graph>
1.46 + /// \param _Graph The graph type the algorithm runs on.
1.47 + template <typename _Graph>
1.48 class MaxMatching {
1.49 -
1.50 - protected:
1.51 + public:
1.52 +
1.53 + typedef _Graph Graph;
1.54 + typedef typename Graph::template NodeMap<typename Graph::Arc>
1.55 + MatchingMap;
1.56 +
1.57 + ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
1.58 + ///
1.59 + ///Indicates the Gallai-Edmonds decomposition of the graph, which
1.60 + ///shows an upper bound on the size of a maximum matching. The
1.61 + ///nodes with Status \c EVEN/D induce a graph with factor-critical
1.62 + ///components, the nodes in \c ODD/A form the canonical barrier,
1.63 + ///and the nodes in \c MATCHED/C induce a graph having a perfect
1.64 + ///matching.
1.65 + enum Status {
1.66 + EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2
1.67 + };
1.68 +
1.69 + typedef typename Graph::template NodeMap<Status> StatusMap;
1.70 +
1.71 + private:
1.72
1.73 TEMPLATE_GRAPH_TYPEDEFS(Graph);
1.74
1.75 - typedef typename Graph::template NodeMap<int> UFECrossRef;
1.76 - typedef UnionFindEnum<UFECrossRef> UFE;
1.77 - typedef std::vector<Node> NV;
1.78 -
1.79 - typedef typename Graph::template NodeMap<int> EFECrossRef;
1.80 - typedef ExtendFindEnum<EFECrossRef> EFE;
1.81 + typedef UnionFindEnum<IntNodeMap> BlossomSet;
1.82 + typedef ExtendFindEnum<IntNodeMap> TreeSet;
1.83 + typedef RangeMap<Node> NodeIntMap;
1.84 + typedef MatchingMap EarMap;
1.85 + typedef std::vector<Node> NodeQueue;
1.86 +
1.87 + const Graph& _graph;
1.88 + MatchingMap* _matching;
1.89 + StatusMap* _status;
1.90 +
1.91 + EarMap* _ear;
1.92 +
1.93 + IntNodeMap* _blossom_set_index;
1.94 + BlossomSet* _blossom_set;
1.95 + NodeIntMap* _blossom_rep;
1.96 +
1.97 + IntNodeMap* _tree_set_index;
1.98 + TreeSet* _tree_set;
1.99 +
1.100 + NodeQueue _node_queue;
1.101 + int _process, _postpone, _last;
1.102 +
1.103 + int _node_num;
1.104 +
1.105 + private:
1.106 +
1.107 + void createStructures() {
1.108 + _node_num = countNodes(_graph);
1.109 + if (!_matching) {
1.110 + _matching = new MatchingMap(_graph);
1.111 + }
1.112 + if (!_status) {
1.113 + _status = new StatusMap(_graph);
1.114 + }
1.115 + if (!_ear) {
1.116 + _ear = new EarMap(_graph);
1.117 + }
1.118 + if (!_blossom_set) {
1.119 + _blossom_set_index = new IntNodeMap(_graph);
1.120 + _blossom_set = new BlossomSet(*_blossom_set_index);
1.121 + }
1.122 + if (!_blossom_rep) {
1.123 + _blossom_rep = new NodeIntMap(_node_num);
1.124 + }
1.125 + if (!_tree_set) {
1.126 + _tree_set_index = new IntNodeMap(_graph);
1.127 + _tree_set = new TreeSet(*_tree_set_index);
1.128 + }
1.129 + _node_queue.resize(_node_num);
1.130 + }
1.131 +
1.132 + void destroyStructures() {
1.133 + if (_matching) {
1.134 + delete _matching;
1.135 + }
1.136 + if (_status) {
1.137 + delete _status;
1.138 + }
1.139 + if (_ear) {
1.140 + delete _ear;
1.141 + }
1.142 + if (_blossom_set) {
1.143 + delete _blossom_set;
1.144 + delete _blossom_set_index;
1.145 + }
1.146 + if (_blossom_rep) {
1.147 + delete _blossom_rep;
1.148 + }
1.149 + if (_tree_set) {
1.150 + delete _tree_set_index;
1.151 + delete _tree_set;
1.152 + }
1.153 + }
1.154 +
1.155 + void processDense(const Node& n) {
1.156 + _process = _postpone = _last = 0;
1.157 + _node_queue[_last++] = n;
1.158 +
1.159 + while (_process != _last) {
1.160 + Node u = _node_queue[_process++];
1.161 + for (OutArcIt a(_graph, u); a != INVALID; ++a) {
1.162 + Node v = _graph.target(a);
1.163 + if ((*_status)[v] == MATCHED) {
1.164 + extendOnArc(a);
1.165 + } else if ((*_status)[v] == UNMATCHED) {
1.166 + augmentOnArc(a);
1.167 + return;
1.168 + }
1.169 + }
1.170 + }
1.171 +
1.172 + while (_postpone != _last) {
1.173 + Node u = _node_queue[_postpone++];
1.174 +
1.175 + for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
1.176 + Node v = _graph.target(a);
1.177 +
1.178 + if ((*_status)[v] == EVEN) {
1.179 + if (_blossom_set->find(u) != _blossom_set->find(v)) {
1.180 + shrinkOnEdge(a);
1.181 + }
1.182 + }
1.183 +
1.184 + while (_process != _last) {
1.185 + Node w = _node_queue[_process++];
1.186 + for (OutArcIt b(_graph, w); b != INVALID; ++b) {
1.187 + Node x = _graph.target(b);
1.188 + if ((*_status)[x] == MATCHED) {
1.189 + extendOnArc(b);
1.190 + } else if ((*_status)[x] == UNMATCHED) {
1.191 + augmentOnArc(b);
1.192 + return;
1.193 + }
1.194 + }
1.195 + }
1.196 + }
1.197 + }
1.198 + }
1.199 +
1.200 + void processSparse(const Node& n) {
1.201 + _process = _last = 0;
1.202 + _node_queue[_last++] = n;
1.203 + while (_process != _last) {
1.204 + Node u = _node_queue[_process++];
1.205 + for (OutArcIt a(_graph, u); a != INVALID; ++a) {
1.206 + Node v = _graph.target(a);
1.207 +
1.208 + if ((*_status)[v] == EVEN) {
1.209 + if (_blossom_set->find(u) != _blossom_set->find(v)) {
1.210 + shrinkOnEdge(a);
1.211 + }
1.212 + } else if ((*_status)[v] == MATCHED) {
1.213 + extendOnArc(a);
1.214 + } else if ((*_status)[v] == UNMATCHED) {
1.215 + augmentOnArc(a);
1.216 + return;
1.217 + }
1.218 + }
1.219 + }
1.220 + }
1.221 +
1.222 + void shrinkOnEdge(const Edge& e) {
1.223 + Node nca = INVALID;
1.224 +
1.225 + {
1.226 + std::set<Node> left_set, right_set;
1.227 +
1.228 + Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
1.229 + left_set.insert(left);
1.230 +
1.231 + Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
1.232 + right_set.insert(right);
1.233 +
1.234 + while (true) {
1.235 + if ((*_matching)[left] == INVALID) break;
1.236 + left = _graph.target((*_matching)[left]);
1.237 + left = (*_blossom_rep)[_blossom_set->
1.238 + find(_graph.target((*_ear)[left]))];
1.239 + if (right_set.find(left) != right_set.end()) {
1.240 + nca = left;
1.241 + break;
1.242 + }
1.243 + left_set.insert(left);
1.244 +
1.245 + if ((*_matching)[right] == INVALID) break;
1.246 + right = _graph.target((*_matching)[right]);
1.247 + right = (*_blossom_rep)[_blossom_set->
1.248 + find(_graph.target((*_ear)[right]))];
1.249 + if (left_set.find(right) != left_set.end()) {
1.250 + nca = right;
1.251 + break;
1.252 + }
1.253 + right_set.insert(right);
1.254 + }
1.255 +
1.256 + if (nca == INVALID) {
1.257 + if ((*_matching)[left] == INVALID) {
1.258 + nca = right;
1.259 + while (left_set.find(nca) == left_set.end()) {
1.260 + nca = _graph.target((*_matching)[nca]);
1.261 + nca =(*_blossom_rep)[_blossom_set->
1.262 + find(_graph.target((*_ear)[nca]))];
1.263 + }
1.264 + } else {
1.265 + nca = left;
1.266 + while (right_set.find(nca) == right_set.end()) {
1.267 + nca = _graph.target((*_matching)[nca]);
1.268 + nca = (*_blossom_rep)[_blossom_set->
1.269 + find(_graph.target((*_ear)[nca]))];
1.270 + }
1.271 + }
1.272 + }
1.273 + }
1.274 +
1.275 + {
1.276 +
1.277 + Node node = _graph.u(e);
1.278 + Arc arc = _graph.direct(e, true);
1.279 + Node base = (*_blossom_rep)[_blossom_set->find(node)];
1.280 +
1.281 + while (base != nca) {
1.282 + _ear->set(node, arc);
1.283 +
1.284 + Node n = node;
1.285 + while (n != base) {
1.286 + n = _graph.target((*_matching)[n]);
1.287 + Arc a = (*_ear)[n];
1.288 + n = _graph.target(a);
1.289 + _ear->set(n, _graph.oppositeArc(a));
1.290 + }
1.291 + node = _graph.target((*_matching)[base]);
1.292 + _tree_set->erase(base);
1.293 + _tree_set->erase(node);
1.294 + _blossom_set->insert(node, _blossom_set->find(base));
1.295 + _status->set(node, EVEN);
1.296 + _node_queue[_last++] = node;
1.297 + arc = _graph.oppositeArc((*_ear)[node]);
1.298 + node = _graph.target((*_ear)[node]);
1.299 + base = (*_blossom_rep)[_blossom_set->find(node)];
1.300 + _blossom_set->join(_graph.target(arc), base);
1.301 + }
1.302 + }
1.303 +
1.304 + _blossom_rep->set(_blossom_set->find(nca), nca);
1.305 +
1.306 + {
1.307 +
1.308 + Node node = _graph.v(e);
1.309 + Arc arc = _graph.direct(e, false);
1.310 + Node base = (*_blossom_rep)[_blossom_set->find(node)];
1.311 +
1.312 + while (base != nca) {
1.313 + _ear->set(node, arc);
1.314 +
1.315 + Node n = node;
1.316 + while (n != base) {
1.317 + n = _graph.target((*_matching)[n]);
1.318 + Arc a = (*_ear)[n];
1.319 + n = _graph.target(a);
1.320 + _ear->set(n, _graph.oppositeArc(a));
1.321 + }
1.322 + node = _graph.target((*_matching)[base]);
1.323 + _tree_set->erase(base);
1.324 + _tree_set->erase(node);
1.325 + _blossom_set->insert(node, _blossom_set->find(base));
1.326 + _status->set(node, EVEN);
1.327 + _node_queue[_last++] = node;
1.328 + arc = _graph.oppositeArc((*_ear)[node]);
1.329 + node = _graph.target((*_ear)[node]);
1.330 + base = (*_blossom_rep)[_blossom_set->find(node)];
1.331 + _blossom_set->join(_graph.target(arc), base);
1.332 + }
1.333 + }
1.334 +
1.335 + _blossom_rep->set(_blossom_set->find(nca), nca);
1.336 + }
1.337 +
1.338 +
1.339 +
1.340 + void extendOnArc(const Arc& a) {
1.341 + Node base = _graph.source(a);
1.342 + Node odd = _graph.target(a);
1.343 +
1.344 + _ear->set(odd, _graph.oppositeArc(a));
1.345 + Node even = _graph.target((*_matching)[odd]);
1.346 + _blossom_rep->set(_blossom_set->insert(even), even);
1.347 + _status->set(odd, ODD);
1.348 + _status->set(even, EVEN);
1.349 + int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
1.350 + _tree_set->insert(odd, tree);
1.351 + _tree_set->insert(even, tree);
1.352 + _node_queue[_last++] = even;
1.353 +
1.354 + }
1.355 +
1.356 + void augmentOnArc(const Arc& a) {
1.357 + Node even = _graph.source(a);
1.358 + Node odd = _graph.target(a);
1.359 +
1.360 + int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
1.361 +
1.362 + _matching->set(odd, _graph.oppositeArc(a));
1.363 + _status->set(odd, MATCHED);
1.364 +
1.365 + Arc arc = (*_matching)[even];
1.366 + _matching->set(even, a);
1.367 +
1.368 + while (arc != INVALID) {
1.369 + odd = _graph.target(arc);
1.370 + arc = (*_ear)[odd];
1.371 + even = _graph.target(arc);
1.372 + _matching->set(odd, arc);
1.373 + arc = (*_matching)[even];
1.374 + _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
1.375 + }
1.376 +
1.377 + for (typename TreeSet::ItemIt it(*_tree_set, tree);
1.378 + it != INVALID; ++it) {
1.379 + if ((*_status)[it] == ODD) {
1.380 + _status->set(it, MATCHED);
1.381 + } else {
1.382 + int blossom = _blossom_set->find(it);
1.383 + for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
1.384 + jt != INVALID; ++jt) {
1.385 + _status->set(jt, MATCHED);
1.386 + }
1.387 + _blossom_set->eraseClass(blossom);
1.388 + }
1.389 + }
1.390 + _tree_set->eraseClass(tree);
1.391 +
1.392 + }
1.393
1.394 public:
1.395
1.396 - ///\brief Indicates the Gallai-Edmonds decomposition of the digraph.
1.397 + /// \brief Constructor
1.398 ///
1.399 - ///Indicates the Gallai-Edmonds decomposition of the digraph, which
1.400 - ///shows an upper bound on the size of a maximum matching. The
1.401 - ///nodes with DecompType \c D induce a digraph with factor-critical
1.402 - ///components, the nodes in \c A form the canonical barrier, and the
1.403 - ///nodes in \c C induce a digraph having a perfect matching.
1.404 - enum DecompType {
1.405 - D=0,
1.406 - A=1,
1.407 - C=2
1.408 - };
1.409 -
1.410 - protected:
1.411 -
1.412 - static const int HEUR_density=2;
1.413 - const Graph& g;
1.414 - typename Graph::template NodeMap<Node> _mate;
1.415 - typename Graph::template NodeMap<DecompType> position;
1.416 -
1.417 - public:
1.418 -
1.419 - MaxMatching(const Graph& _g)
1.420 - : g(_g), _mate(_g), position(_g) {}
1.421 -
1.422 - ///\brief Sets the actual matching to the empty matching.
1.423 + /// Constructor.
1.424 + MaxMatching(const Graph& graph)
1.425 + : _graph(graph), _matching(0), _status(0), _ear(0),
1.426 + _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
1.427 + _tree_set_index(0), _tree_set(0) {}
1.428 +
1.429 + ~MaxMatching() {
1.430 + destroyStructures();
1.431 + }
1.432 +
1.433 + /// \name Execution control
1.434 + /// The simplest way to execute the algorithm is to use the member
1.435 + /// \c run() member function.
1.436 + /// \n
1.437 +
1.438 + /// If you need more control on the execution, first you must call
1.439 + /// \ref init(), \ref greedyInit() or \ref matchingInit()
1.440 + /// functions, then you can start the algorithm with the \ref
1.441 + /// startParse() or startDense() functions.
1.442 +
1.443 + ///@{
1.444 +
1.445 + /// \brief Sets the actual matching to the empty matching.
1.446 ///
1.447 - ///Sets the actual matching to the empty matching.
1.448 + /// Sets the actual matching to the empty matching.
1.449 ///
1.450 void init() {
1.451 - for(NodeIt v(g); v!=INVALID; ++v) {
1.452 - _mate.set(v,INVALID);
1.453 - position.set(v,C);
1.454 + createStructures();
1.455 + for(NodeIt n(_graph); n != INVALID; ++n) {
1.456 + _matching->set(n, INVALID);
1.457 + _status->set(n, UNMATCHED);
1.458 }
1.459 }
1.460
1.461 @@ -108,88 +437,96 @@
1.462 ///
1.463 ///For initial matchig it finds a maximal greedy matching.
1.464 void greedyInit() {
1.465 - for(NodeIt v(g); v!=INVALID; ++v) {
1.466 - _mate.set(v,INVALID);
1.467 - position.set(v,C);
1.468 + createStructures();
1.469 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.470 + _matching->set(n, INVALID);
1.471 + _status->set(n, UNMATCHED);
1.472 }
1.473 - for(NodeIt v(g); v!=INVALID; ++v)
1.474 - if ( _mate[v]==INVALID ) {
1.475 - for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) {
1.476 - Node y=g.runningNode(e);
1.477 - if ( _mate[y]==INVALID && y!=v ) {
1.478 - _mate.set(v,y);
1.479 - _mate.set(y,v);
1.480 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.481 + if ((*_matching)[n] == INVALID) {
1.482 + for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
1.483 + Node v = _graph.target(a);
1.484 + if ((*_matching)[v] == INVALID && v != n) {
1.485 + _matching->set(n, a);
1.486 + _status->set(n, MATCHED);
1.487 + _matching->set(v, _graph.oppositeArc(a));
1.488 + _status->set(v, MATCHED);
1.489 break;
1.490 }
1.491 }
1.492 }
1.493 - }
1.494 -
1.495 - ///\brief Initialize the matching from each nodes' mate.
1.496 - ///
1.497 - ///Initialize the matching from a \c Node valued \c Node map. This
1.498 - ///map must be \e symmetric, i.e. if \c map[u]==v then \c
1.499 - ///map[v]==u must hold, and \c uv will be an arc of the initial
1.500 - ///matching.
1.501 - template <typename MateMap>
1.502 - void mateMapInit(MateMap& map) {
1.503 - for(NodeIt v(g); v!=INVALID; ++v) {
1.504 - _mate.set(v,map[v]);
1.505 - position.set(v,C);
1.506 }
1.507 }
1.508
1.509 - ///\brief Initialize the matching from a node map with the
1.510 - ///incident matching arcs.
1.511 +
1.512 + /// \brief Initialize the matching from the map containing a matching.
1.513 ///
1.514 - ///Initialize the matching from an \c Edge valued \c Node map. \c
1.515 - ///map[v] must be an \c Edge incident to \c v. This map must have
1.516 - ///the property that if \c g.oppositeNode(u,map[u])==v then \c \c
1.517 - ///g.oppositeNode(v,map[v])==u holds, and now some arc joining \c
1.518 - ///u to \c v will be an arc of the matching.
1.519 - template<typename MatchingMap>
1.520 - void matchingMapInit(MatchingMap& map) {
1.521 - for(NodeIt v(g); v!=INVALID; ++v) {
1.522 - position.set(v,C);
1.523 - Edge e=map[v];
1.524 - if ( e!=INVALID )
1.525 - _mate.set(v,g.oppositeNode(v,e));
1.526 - else
1.527 - _mate.set(v,INVALID);
1.528 + /// Initialize the matching from a \c bool valued \c Edge map. This
1.529 + /// map must have the property that there are no two incident edges
1.530 + /// with true value, ie. it contains a matching.
1.531 + /// \return %True if the map contains a matching.
1.532 + template <typename MatchingMap>
1.533 + bool matchingInit(const MatchingMap& matching) {
1.534 + createStructures();
1.535 +
1.536 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.537 + _matching->set(n, INVALID);
1.538 + _status->set(n, UNMATCHED);
1.539 }
1.540 + for(EdgeIt e(_graph); e!=INVALID; ++e) {
1.541 + if (matching[e]) {
1.542 +
1.543 + Node u = _graph.u(e);
1.544 + if ((*_matching)[u] != INVALID) return false;
1.545 + _matching->set(u, _graph.direct(e, true));
1.546 + _status->set(u, MATCHED);
1.547 +
1.548 + Node v = _graph.v(e);
1.549 + if ((*_matching)[v] != INVALID) return false;
1.550 + _matching->set(v, _graph.direct(e, false));
1.551 + _status->set(v, MATCHED);
1.552 + }
1.553 + }
1.554 + return true;
1.555 }
1.556
1.557 - ///\brief Initialize the matching from the map containing the
1.558 - ///undirected matching arcs.
1.559 + /// \brief Starts Edmonds' algorithm
1.560 ///
1.561 - ///Initialize the matching from a \c bool valued \c Edge map. This
1.562 - ///map must have the property that there are no two incident arcs
1.563 - ///\c e, \c f with \c map[e]==map[f]==true. The arcs \c e with \c
1.564 - ///map[e]==true form the matching.
1.565 - template <typename MatchingMap>
1.566 - void matchingInit(MatchingMap& map) {
1.567 - for(NodeIt v(g); v!=INVALID; ++v) {
1.568 - _mate.set(v,INVALID);
1.569 - position.set(v,C);
1.570 - }
1.571 - for(EdgeIt e(g); e!=INVALID; ++e) {
1.572 - if ( map[e] ) {
1.573 - Node u=g.u(e);
1.574 - Node v=g.v(e);
1.575 - _mate.set(u,v);
1.576 - _mate.set(v,u);
1.577 + /// If runs the original Edmonds' algorithm.
1.578 + void startSparse() {
1.579 + for(NodeIt n(_graph); n != INVALID; ++n) {
1.580 + if ((*_status)[n] == UNMATCHED) {
1.581 + (*_blossom_rep)[_blossom_set->insert(n)] = n;
1.582 + _tree_set->insert(n);
1.583 + _status->set(n, EVEN);
1.584 + processSparse(n);
1.585 }
1.586 }
1.587 }
1.588
1.589 -
1.590 - ///\brief Runs Edmonds' algorithm.
1.591 + /// \brief Starts Edmonds' algorithm.
1.592 ///
1.593 - ///Runs Edmonds' algorithm for sparse digraphs (number of arcs <
1.594 - ///2*number of nodes), and a heuristical Edmonds' algorithm with a
1.595 - ///heuristic of postponing shrinks for dense digraphs.
1.596 + /// It runs Edmonds' algorithm with a heuristic of postponing
1.597 + /// shrinks, giving a faster algorithm for dense graphs.
1.598 + void startDense() {
1.599 + for(NodeIt n(_graph); n != INVALID; ++n) {
1.600 + if ((*_status)[n] == UNMATCHED) {
1.601 + (*_blossom_rep)[_blossom_set->insert(n)] = n;
1.602 + _tree_set->insert(n);
1.603 + _status->set(n, EVEN);
1.604 + processDense(n);
1.605 + }
1.606 + }
1.607 + }
1.608 +
1.609 +
1.610 + /// \brief Runs Edmonds' algorithm
1.611 + ///
1.612 + /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>)
1.613 + /// or Edmonds' algorithm with a heuristic of
1.614 + /// postponing shrinks for dense graphs.
1.615 void run() {
1.616 - if (countEdges(g) < HEUR_density * countNodes(g)) {
1.617 + if (countEdges(_graph) < 2 * countNodes(_graph)) {
1.618 greedyInit();
1.619 startSparse();
1.620 } else {
1.621 @@ -198,422 +535,75 @@
1.622 }
1.623 }
1.624
1.625 -
1.626 - ///\brief Starts Edmonds' algorithm.
1.627 - ///
1.628 - ///If runs the original Edmonds' algorithm.
1.629 - void startSparse() {
1.630 -
1.631 - typename Graph::template NodeMap<Node> ear(g,INVALID);
1.632 - //undefined for the base nodes of the blossoms (i.e. for the
1.633 - //representative elements of UFE blossom) and for the nodes in C
1.634 -
1.635 - UFECrossRef blossom_base(g);
1.636 - UFE blossom(blossom_base);
1.637 - NV rep(countNodes(g));
1.638 -
1.639 - EFECrossRef tree_base(g);
1.640 - EFE tree(tree_base);
1.641 -
1.642 - //If these UFE's would be members of the class then also
1.643 - //blossom_base and tree_base should be a member.
1.644 -
1.645 - //We build only one tree and the other vertices uncovered by the
1.646 - //matching belong to C. (They can be considered as singleton
1.647 - //trees.) If this tree can be augmented or no more
1.648 - //grow/augmentation/shrink is possible then we return to this
1.649 - //"for" cycle.
1.650 - for(NodeIt v(g); v!=INVALID; ++v) {
1.651 - if (position[v]==C && _mate[v]==INVALID) {
1.652 - rep[blossom.insert(v)] = v;
1.653 - tree.insert(v);
1.654 - position.set(v,D);
1.655 - normShrink(v, ear, blossom, rep, tree);
1.656 - }
1.657 - }
1.658 - }
1.659 -
1.660 - ///\brief Starts Edmonds' algorithm.
1.661 - ///
1.662 - ///It runs Edmonds' algorithm with a heuristic of postponing
1.663 - ///shrinks, giving a faster algorithm for dense digraphs.
1.664 - void startDense() {
1.665 -
1.666 - typename Graph::template NodeMap<Node> ear(g,INVALID);
1.667 - //undefined for the base nodes of the blossoms (i.e. for the
1.668 - //representative elements of UFE blossom) and for the nodes in C
1.669 -
1.670 - UFECrossRef blossom_base(g);
1.671 - UFE blossom(blossom_base);
1.672 - NV rep(countNodes(g));
1.673 -
1.674 - EFECrossRef tree_base(g);
1.675 - EFE tree(tree_base);
1.676 -
1.677 - //If these UFE's would be members of the class then also
1.678 - //blossom_base and tree_base should be a member.
1.679 -
1.680 - //We build only one tree and the other vertices uncovered by the
1.681 - //matching belong to C. (They can be considered as singleton
1.682 - //trees.) If this tree can be augmented or no more
1.683 - //grow/augmentation/shrink is possible then we return to this
1.684 - //"for" cycle.
1.685 - for(NodeIt v(g); v!=INVALID; ++v) {
1.686 - if ( position[v]==C && _mate[v]==INVALID ) {
1.687 - rep[blossom.insert(v)] = v;
1.688 - tree.insert(v);
1.689 - position.set(v,D);
1.690 - lateShrink(v, ear, blossom, rep, tree);
1.691 - }
1.692 - }
1.693 - }
1.694 -
1.695 -
1.696 + /// @}
1.697 +
1.698 + /// \name Primal solution
1.699 + /// Functions for get the primal solution, ie. the matching.
1.700 +
1.701 + /// @{
1.702
1.703 ///\brief Returns the size of the actual matching stored.
1.704 ///
1.705 ///Returns the size of the actual matching stored. After \ref
1.706 - ///run() it returns the size of a maximum matching in the digraph.
1.707 - int size() const {
1.708 - int s=0;
1.709 - for(NodeIt v(g); v!=INVALID; ++v) {
1.710 - if ( _mate[v]!=INVALID ) {
1.711 - ++s;
1.712 + ///run() it returns the size of the maximum matching in the graph.
1.713 + int matchingSize() const {
1.714 + int size = 0;
1.715 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.716 + if ((*_matching)[n] != INVALID) {
1.717 + ++size;
1.718 }
1.719 }
1.720 - return s/2;
1.721 + return size / 2;
1.722 }
1.723
1.724 + /// \brief Returns true when the edge is in the matching.
1.725 + ///
1.726 + /// Returns true when the edge is in the matching.
1.727 + bool matching(const Edge& edge) const {
1.728 + return edge == (*_matching)[_graph.u(edge)];
1.729 + }
1.730 +
1.731 + /// \brief Returns the matching edge incident to the given node.
1.732 + ///
1.733 + /// Returns the matching edge of a \c node in the actual matching or
1.734 + /// INVALID if the \c node is not covered by the actual matching.
1.735 + Arc matching(const Node& n) const {
1.736 + return (*_matching)[n];
1.737 + }
1.738
1.739 ///\brief Returns the mate of a node in the actual matching.
1.740 ///
1.741 - ///Returns the mate of a \c node in the actual matching.
1.742 - ///Returns INVALID if the \c node is not covered by the actual matching.
1.743 - Node mate(const Node& node) const {
1.744 - return _mate[node];
1.745 + ///Returns the mate of a \c node in the actual matching or
1.746 + ///INVALID if the \c node is not covered by the actual matching.
1.747 + Node mate(const Node& n) const {
1.748 + return (*_matching)[n] != INVALID ?
1.749 + _graph.target((*_matching)[n]) : INVALID;
1.750 }
1.751
1.752 - ///\brief Returns the matching arc incident to the given node.
1.753 - ///
1.754 - ///Returns the matching arc of a \c node in the actual matching.
1.755 - ///Returns INVALID if the \c node is not covered by the actual matching.
1.756 - Edge matchingArc(const Node& node) const {
1.757 - if (_mate[node] == INVALID) return INVALID;
1.758 - Node n = node < _mate[node] ? node : _mate[node];
1.759 - for (IncEdgeIt e(g, n); e != INVALID; ++e) {
1.760 - if (g.oppositeNode(n, e) == _mate[n]) {
1.761 - return e;
1.762 - }
1.763 - }
1.764 - return INVALID;
1.765 - }
1.766 + /// @}
1.767 +
1.768 + /// \name Dual solution
1.769 + /// Functions for get the dual solution, ie. the decomposition.
1.770 +
1.771 + /// @{
1.772
1.773 /// \brief Returns the class of the node in the Edmonds-Gallai
1.774 /// decomposition.
1.775 ///
1.776 /// Returns the class of the node in the Edmonds-Gallai
1.777 /// decomposition.
1.778 - DecompType decomposition(const Node& n) {
1.779 - return position[n] == A;
1.780 + Status decomposition(const Node& n) const {
1.781 + return (*_status)[n];
1.782 }
1.783
1.784 /// \brief Returns true when the node is in the barrier.
1.785 ///
1.786 /// Returns true when the node is in the barrier.
1.787 - bool barrier(const Node& n) {
1.788 - return position[n] == A;
1.789 + bool barrier(const Node& n) const {
1.790 + return (*_status)[n] == ODD;
1.791 }
1.792
1.793 - ///\brief Gives back the matching in a \c Node of mates.
1.794 - ///
1.795 - ///Writes the stored matching to a \c Node valued \c Node map. The
1.796 - ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c
1.797 - ///map[v]==u will hold, and now \c uv is an arc of the matching.
1.798 - template <typename MateMap>
1.799 - void mateMap(MateMap& map) const {
1.800 - for(NodeIt v(g); v!=INVALID; ++v) {
1.801 - map.set(v,_mate[v]);
1.802 - }
1.803 - }
1.804 -
1.805 - ///\brief Gives back the matching in an \c Edge valued \c Node
1.806 - ///map.
1.807 - ///
1.808 - ///Writes the stored matching to an \c Edge valued \c Node
1.809 - ///map. \c map[v] will be an \c Edge incident to \c v. This
1.810 - ///map will have the property that if \c g.oppositeNode(u,map[u])
1.811 - ///== v then \c map[u]==map[v] holds, and now this arc is an arc
1.812 - ///of the matching.
1.813 - template <typename MatchingMap>
1.814 - void matchingMap(MatchingMap& map) const {
1.815 - typename Graph::template NodeMap<bool> todo(g,true);
1.816 - for(NodeIt v(g); v!=INVALID; ++v) {
1.817 - if (_mate[v]!=INVALID && v < _mate[v]) {
1.818 - Node u=_mate[v];
1.819 - for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
1.820 - if ( g.runningNode(e) == u ) {
1.821 - map.set(u,e);
1.822 - map.set(v,e);
1.823 - todo.set(u,false);
1.824 - todo.set(v,false);
1.825 - break;
1.826 - }
1.827 - }
1.828 - }
1.829 - }
1.830 - }
1.831 -
1.832 -
1.833 - ///\brief Gives back the matching in a \c bool valued \c Edge
1.834 - ///map.
1.835 - ///
1.836 - ///Writes the matching stored to a \c bool valued \c Arc
1.837 - ///map. This map will have the property that there are no two
1.838 - ///incident arcs \c e, \c f with \c map[e]==map[f]==true. The
1.839 - ///arcs \c e with \c map[e]==true form the matching.
1.840 - template<typename MatchingMap>
1.841 - void matching(MatchingMap& map) const {
1.842 - for(EdgeIt e(g); e!=INVALID; ++e) map.set(e,false);
1.843 -
1.844 - typename Graph::template NodeMap<bool> todo(g,true);
1.845 - for(NodeIt v(g); v!=INVALID; ++v) {
1.846 - if ( todo[v] && _mate[v]!=INVALID ) {
1.847 - Node u=_mate[v];
1.848 - for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
1.849 - if ( g.runningNode(e) == u ) {
1.850 - map.set(e,true);
1.851 - todo.set(u,false);
1.852 - todo.set(v,false);
1.853 - break;
1.854 - }
1.855 - }
1.856 - }
1.857 - }
1.858 - }
1.859 -
1.860 -
1.861 - ///\brief Returns the canonical decomposition of the digraph after running
1.862 - ///the algorithm.
1.863 - ///
1.864 - ///After calling any run methods of the class, it writes the
1.865 - ///Gallai-Edmonds canonical decomposition of the digraph. \c map
1.866 - ///must be a node map of \ref DecompType 's.
1.867 - template <typename DecompositionMap>
1.868 - void decomposition(DecompositionMap& map) const {
1.869 - for(NodeIt v(g); v!=INVALID; ++v) map.set(v,position[v]);
1.870 - }
1.871 -
1.872 - ///\brief Returns a barrier on the nodes.
1.873 - ///
1.874 - ///After calling any run methods of the class, it writes a
1.875 - ///canonical barrier on the nodes. The odd component number of the
1.876 - ///remaining digraph minus the barrier size is a lower bound for the
1.877 - ///uncovered nodes in the digraph. The \c map must be a node map of
1.878 - ///bools.
1.879 - template <typename BarrierMap>
1.880 - void barrier(BarrierMap& barrier) {
1.881 - for(NodeIt v(g); v!=INVALID; ++v) barrier.set(v,position[v] == A);
1.882 - }
1.883 -
1.884 - private:
1.885 -
1.886 -
1.887 - void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear,
1.888 - UFE& blossom, NV& rep, EFE& tree) {
1.889 - //We have one tree which we grow, and also shrink but only if it
1.890 - //cannot be postponed. If we augment then we return to the "for"
1.891 - //cycle of runEdmonds().
1.892 -
1.893 - std::queue<Node> Q; //queue of the totally unscanned nodes
1.894 - Q.push(v);
1.895 - std::queue<Node> R;
1.896 - //queue of the nodes which must be scanned for a possible shrink
1.897 -
1.898 - while ( !Q.empty() ) {
1.899 - Node x=Q.front();
1.900 - Q.pop();
1.901 - for( IncEdgeIt e(g,x); e!= INVALID; ++e ) {
1.902 - Node y=g.runningNode(e);
1.903 - //growOrAugment grows if y is covered by the matching and
1.904 - //augments if not. In this latter case it returns 1.
1.905 - if (position[y]==C &&
1.906 - growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
1.907 - }
1.908 - R.push(x);
1.909 - }
1.910 -
1.911 - while ( !R.empty() ) {
1.912 - Node x=R.front();
1.913 - R.pop();
1.914 -
1.915 - for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) {
1.916 - Node y=g.runningNode(e);
1.917 -
1.918 - if ( position[y] == D && blossom.find(x) != blossom.find(y) )
1.919 - //Recall that we have only one tree.
1.920 - shrink( x, y, ear, blossom, rep, tree, Q);
1.921 -
1.922 - while ( !Q.empty() ) {
1.923 - Node z=Q.front();
1.924 - Q.pop();
1.925 - for( IncEdgeIt f(g,z); f!= INVALID; ++f ) {
1.926 - Node w=g.runningNode(f);
1.927 - //growOrAugment grows if y is covered by the matching and
1.928 - //augments if not. In this latter case it returns 1.
1.929 - if (position[w]==C &&
1.930 - growOrAugment(w, z, ear, blossom, rep, tree, Q)) return;
1.931 - }
1.932 - R.push(z);
1.933 - }
1.934 - } //for e
1.935 - } // while ( !R.empty() )
1.936 - }
1.937 -
1.938 - void normShrink(Node v, typename Graph::template NodeMap<Node>& ear,
1.939 - UFE& blossom, NV& rep, EFE& tree) {
1.940 - //We have one tree, which we grow and shrink. If we augment then we
1.941 - //return to the "for" cycle of runEdmonds().
1.942 -
1.943 - std::queue<Node> Q; //queue of the unscanned nodes
1.944 - Q.push(v);
1.945 - while ( !Q.empty() ) {
1.946 -
1.947 - Node x=Q.front();
1.948 - Q.pop();
1.949 -
1.950 - for( IncEdgeIt e(g,x); e!=INVALID; ++e ) {
1.951 - Node y=g.runningNode(e);
1.952 -
1.953 - switch ( position[y] ) {
1.954 - case D: //x and y must be in the same tree
1.955 - if ( blossom.find(x) != blossom.find(y))
1.956 - //x and y are in the same tree
1.957 - shrink(x, y, ear, blossom, rep, tree, Q);
1.958 - break;
1.959 - case C:
1.960 - //growOrAugment grows if y is covered by the matching and
1.961 - //augments if not. In this latter case it returns 1.
1.962 - if (growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
1.963 - break;
1.964 - default: break;
1.965 - }
1.966 - }
1.967 - }
1.968 - }
1.969 -
1.970 - void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear,
1.971 - UFE& blossom, NV& rep, EFE& tree,std::queue<Node>& Q) {
1.972 - //x and y are the two adjacent vertices in two blossoms.
1.973 -
1.974 - typename Graph::template NodeMap<bool> path(g,false);
1.975 -
1.976 - Node b=rep[blossom.find(x)];
1.977 - path.set(b,true);
1.978 - b=_mate[b];
1.979 - while ( b!=INVALID ) {
1.980 - b=rep[blossom.find(ear[b])];
1.981 - path.set(b,true);
1.982 - b=_mate[b];
1.983 - } //we go until the root through bases of blossoms and odd vertices
1.984 -
1.985 - Node top=y;
1.986 - Node middle=rep[blossom.find(top)];
1.987 - Node bottom=x;
1.988 - while ( !path[middle] )
1.989 - shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
1.990 - //Until we arrive to a node on the path, we update blossom, tree
1.991 - //and the positions of the odd nodes.
1.992 -
1.993 - Node base=middle;
1.994 - top=x;
1.995 - middle=rep[blossom.find(top)];
1.996 - bottom=y;
1.997 - Node blossom_base=rep[blossom.find(base)];
1.998 - while ( middle!=blossom_base )
1.999 - shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
1.1000 - //Until we arrive to a node on the path, we update blossom, tree
1.1001 - //and the positions of the odd nodes.
1.1002 -
1.1003 - rep[blossom.find(base)] = base;
1.1004 - }
1.1005 -
1.1006 - void shrinkStep(Node& top, Node& middle, Node& bottom,
1.1007 - typename Graph::template NodeMap<Node>& ear,
1.1008 - UFE& blossom, NV& rep, EFE& tree, std::queue<Node>& Q) {
1.1009 - //We traverse a blossom and update everything.
1.1010 -
1.1011 - ear.set(top,bottom);
1.1012 - Node t=top;
1.1013 - while ( t!=middle ) {
1.1014 - Node u=_mate[t];
1.1015 - t=ear[u];
1.1016 - ear.set(t,u);
1.1017 - }
1.1018 - bottom=_mate[middle];
1.1019 - position.set(bottom,D);
1.1020 - Q.push(bottom);
1.1021 - top=ear[bottom];
1.1022 - Node oldmiddle=middle;
1.1023 - middle=rep[blossom.find(top)];
1.1024 - tree.erase(bottom);
1.1025 - tree.erase(oldmiddle);
1.1026 - blossom.insert(bottom);
1.1027 - blossom.join(bottom, oldmiddle);
1.1028 - blossom.join(top, oldmiddle);
1.1029 - }
1.1030 -
1.1031 -
1.1032 -
1.1033 - bool growOrAugment(Node& y, Node& x, typename Graph::template
1.1034 - NodeMap<Node>& ear, UFE& blossom, NV& rep, EFE& tree,
1.1035 - std::queue<Node>& Q) {
1.1036 - //x is in a blossom in the tree, y is outside. If y is covered by
1.1037 - //the matching we grow, otherwise we augment. In this case we
1.1038 - //return 1.
1.1039 -
1.1040 - if ( _mate[y]!=INVALID ) { //grow
1.1041 - ear.set(y,x);
1.1042 - Node w=_mate[y];
1.1043 - rep[blossom.insert(w)] = w;
1.1044 - position.set(y,A);
1.1045 - position.set(w,D);
1.1046 - int t = tree.find(rep[blossom.find(x)]);
1.1047 - tree.insert(y,t);
1.1048 - tree.insert(w,t);
1.1049 - Q.push(w);
1.1050 - } else { //augment
1.1051 - augment(x, ear, blossom, rep, tree);
1.1052 - _mate.set(x,y);
1.1053 - _mate.set(y,x);
1.1054 - return true;
1.1055 - }
1.1056 - return false;
1.1057 - }
1.1058 -
1.1059 - void augment(Node x, typename Graph::template NodeMap<Node>& ear,
1.1060 - UFE& blossom, NV& rep, EFE& tree) {
1.1061 - Node v=_mate[x];
1.1062 - while ( v!=INVALID ) {
1.1063 -
1.1064 - Node u=ear[v];
1.1065 - _mate.set(v,u);
1.1066 - Node tmp=v;
1.1067 - v=_mate[u];
1.1068 - _mate.set(u,tmp);
1.1069 - }
1.1070 - int y = tree.find(rep[blossom.find(x)]);
1.1071 - for (typename EFE::ItemIt tit(tree, y); tit != INVALID; ++tit) {
1.1072 - if ( position[tit] == D ) {
1.1073 - int b = blossom.find(tit);
1.1074 - for (typename UFE::ItemIt bit(blossom, b); bit != INVALID; ++bit) {
1.1075 - position.set(bit, C);
1.1076 - }
1.1077 - blossom.eraseClass(b);
1.1078 - } else position.set(tit, C);
1.1079 - }
1.1080 - tree.eraseClass(y);
1.1081 -
1.1082 - }
1.1083 + /// @}
1.1084
1.1085 };
1.1086
1.1087 @@ -627,25 +617,28 @@
1.1088 /// \f$O(nm\log(n))\f$ time complexity.
1.1089 ///
1.1090 /// The maximum weighted matching problem is to find undirected
1.1091 - /// arcs in the digraph with maximum overall weight and no two of
1.1092 - /// them shares their endpoints. The problem can be formulated with
1.1093 - /// the next linear program:
1.1094 + /// edges in the graph with maximum overall weight and no two of
1.1095 + /// them shares their ends. The problem can be formulated with the
1.1096 + /// following linear program.
1.1097 /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
1.1098 - ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
1.1099 + /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
1.1100 + \quad \forall B\in\mathcal{O}\f] */
1.1101 /// \f[x_e \ge 0\quad \forall e\in E\f]
1.1102 /// \f[\max \sum_{e\in E}x_ew_e\f]
1.1103 - /// where \f$\delta(X)\f$ is the set of arcs incident to a node in
1.1104 - /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in
1.1105 - /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
1.1106 - /// the nodes.
1.1107 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1.1108 + /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
1.1109 + /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
1.1110 + /// subsets of the nodes.
1.1111 ///
1.1112 /// The algorithm calculates an optimal matching and a proof of the
1.1113 /// optimality. The solution of the dual problem can be used to check
1.1114 - /// the result of the algorithm. The dual linear problem is the next:
1.1115 - /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f]
1.1116 + /// the result of the algorithm. The dual linear problem is the
1.1117 + /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
1.1118 + z_B \ge w_{uv} \quad \forall uv\in E\f] */
1.1119 /// \f[y_u \ge 0 \quad \forall u \in V\f]
1.1120 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
1.1121 - /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
1.1122 + /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
1.1123 + \frac{\vert B \vert - 1}{2}z_B\f] */
1.1124 ///
1.1125 /// The algorithm can be executed with \c run() or the \c init() and
1.1126 /// then the \c start() member functions. After it the matching can
1.1127 @@ -705,16 +698,13 @@
1.1128 int _node_num;
1.1129 int _blossom_num;
1.1130
1.1131 - typedef typename Graph::template NodeMap<int> NodeIntMap;
1.1132 - typedef typename Graph::template ArcMap<int> ArcIntMap;
1.1133 - typedef typename Graph::template EdgeMap<int> EdgeIntMap;
1.1134 typedef RangeMap<int> IntIntMap;
1.1135
1.1136 enum Status {
1.1137 EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
1.1138 };
1.1139
1.1140 - typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
1.1141 + typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
1.1142 struct BlossomData {
1.1143 int tree;
1.1144 Status status;
1.1145 @@ -723,21 +713,21 @@
1.1146 Node base;
1.1147 };
1.1148
1.1149 - NodeIntMap *_blossom_index;
1.1150 + IntNodeMap *_blossom_index;
1.1151 BlossomSet *_blossom_set;
1.1152 RangeMap<BlossomData>* _blossom_data;
1.1153
1.1154 - NodeIntMap *_node_index;
1.1155 - ArcIntMap *_node_heap_index;
1.1156 + IntNodeMap *_node_index;
1.1157 + IntArcMap *_node_heap_index;
1.1158
1.1159 struct NodeData {
1.1160
1.1161 - NodeData(ArcIntMap& node_heap_index)
1.1162 + NodeData(IntArcMap& node_heap_index)
1.1163 : heap(node_heap_index) {}
1.1164
1.1165 int blossom;
1.1166 Value pot;
1.1167 - BinHeap<Value, ArcIntMap> heap;
1.1168 + BinHeap<Value, IntArcMap> heap;
1.1169 std::map<int, Arc> heap_index;
1.1170
1.1171 int tree;
1.1172 @@ -750,14 +740,14 @@
1.1173 IntIntMap *_tree_set_index;
1.1174 TreeSet *_tree_set;
1.1175
1.1176 - NodeIntMap *_delta1_index;
1.1177 - BinHeap<Value, NodeIntMap> *_delta1;
1.1178 + IntNodeMap *_delta1_index;
1.1179 + BinHeap<Value, IntNodeMap> *_delta1;
1.1180
1.1181 IntIntMap *_delta2_index;
1.1182 BinHeap<Value, IntIntMap> *_delta2;
1.1183
1.1184 - EdgeIntMap *_delta3_index;
1.1185 - BinHeap<Value, EdgeIntMap> *_delta3;
1.1186 + IntEdgeMap *_delta3_index;
1.1187 + BinHeap<Value, IntEdgeMap> *_delta3;
1.1188
1.1189 IntIntMap *_delta4_index;
1.1190 BinHeap<Value, IntIntMap> *_delta4;
1.1191 @@ -775,14 +765,14 @@
1.1192 _node_potential = new NodePotential(_graph);
1.1193 }
1.1194 if (!_blossom_set) {
1.1195 - _blossom_index = new NodeIntMap(_graph);
1.1196 + _blossom_index = new IntNodeMap(_graph);
1.1197 _blossom_set = new BlossomSet(*_blossom_index);
1.1198 _blossom_data = new RangeMap<BlossomData>(_blossom_num);
1.1199 }
1.1200
1.1201 if (!_node_index) {
1.1202 - _node_index = new NodeIntMap(_graph);
1.1203 - _node_heap_index = new ArcIntMap(_graph);
1.1204 + _node_index = new IntNodeMap(_graph);
1.1205 + _node_heap_index = new IntArcMap(_graph);
1.1206 _node_data = new RangeMap<NodeData>(_node_num,
1.1207 NodeData(*_node_heap_index));
1.1208 }
1.1209 @@ -792,16 +782,16 @@
1.1210 _tree_set = new TreeSet(*_tree_set_index);
1.1211 }
1.1212 if (!_delta1) {
1.1213 - _delta1_index = new NodeIntMap(_graph);
1.1214 - _delta1 = new BinHeap<Value, NodeIntMap>(*_delta1_index);
1.1215 + _delta1_index = new IntNodeMap(_graph);
1.1216 + _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
1.1217 }
1.1218 if (!_delta2) {
1.1219 _delta2_index = new IntIntMap(_blossom_num);
1.1220 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
1.1221 }
1.1222 if (!_delta3) {
1.1223 - _delta3_index = new EdgeIntMap(_graph);
1.1224 - _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index);
1.1225 + _delta3_index = new IntEdgeMap(_graph);
1.1226 + _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
1.1227 }
1.1228 if (!_delta4) {
1.1229 _delta4_index = new IntIntMap(_blossom_num);
1.1230 @@ -1266,10 +1256,10 @@
1.1231 }
1.1232
1.1233
1.1234 - void augmentOnArc(const Edge& arc) {
1.1235 -
1.1236 - int left = _blossom_set->find(_graph.u(arc));
1.1237 - int right = _blossom_set->find(_graph.v(arc));
1.1238 + void augmentOnEdge(const Edge& edge) {
1.1239 +
1.1240 + int left = _blossom_set->find(_graph.u(edge));
1.1241 + int right = _blossom_set->find(_graph.v(edge));
1.1242
1.1243 if ((*_blossom_data)[left].status == EVEN) {
1.1244 int left_tree = _tree_set->find(left);
1.1245 @@ -1289,8 +1279,8 @@
1.1246 unmatchedToMatched(right);
1.1247 }
1.1248
1.1249 - (*_blossom_data)[left].next = _graph.direct(arc, true);
1.1250 - (*_blossom_data)[right].next = _graph.direct(arc, false);
1.1251 + (*_blossom_data)[left].next = _graph.direct(edge, true);
1.1252 + (*_blossom_data)[right].next = _graph.direct(edge, false);
1.1253 }
1.1254
1.1255 void extendOnArc(const Arc& arc) {
1.1256 @@ -1310,7 +1300,7 @@
1.1257 matchedToEven(even, tree);
1.1258 }
1.1259
1.1260 - void shrinkOnArc(const Edge& edge, int tree) {
1.1261 + void shrinkOnEdge(const Edge& edge, int tree) {
1.1262 int nca = -1;
1.1263 std::vector<int> left_path, right_path;
1.1264
1.1265 @@ -1652,7 +1642,7 @@
1.1266 createStructures();
1.1267
1.1268 for (ArcIt e(_graph); e != INVALID; ++e) {
1.1269 - _node_heap_index->set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP);
1.1270 + _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
1.1271 }
1.1272 for (NodeIt n(_graph); n != INVALID; ++n) {
1.1273 _delta1_index->set(n, _delta1->PRE_HEAP);
1.1274 @@ -1769,9 +1759,9 @@
1.1275 }
1.1276
1.1277 if (left_tree == right_tree) {
1.1278 - shrinkOnArc(e, left_tree);
1.1279 + shrinkOnEdge(e, left_tree);
1.1280 } else {
1.1281 - augmentOnArc(e);
1.1282 + augmentOnEdge(e);
1.1283 unmatched -= 2;
1.1284 }
1.1285 }
1.1286 @@ -1818,11 +1808,24 @@
1.1287 return sum /= 2;
1.1288 }
1.1289
1.1290 - /// \brief Returns true when the arc is in the matching.
1.1291 + /// \brief Returns the cardinality of the matching.
1.1292 ///
1.1293 - /// Returns true when the arc is in the matching.
1.1294 - bool matching(const Edge& arc) const {
1.1295 - return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true);
1.1296 + /// Returns the cardinality of the matching.
1.1297 + int matchingSize() const {
1.1298 + int num = 0;
1.1299 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1300 + if ((*_matching)[n] != INVALID) {
1.1301 + ++num;
1.1302 + }
1.1303 + }
1.1304 + return num /= 2;
1.1305 + }
1.1306 +
1.1307 + /// \brief Returns true when the edge is in the matching.
1.1308 + ///
1.1309 + /// Returns true when the edge is in the matching.
1.1310 + bool matching(const Edge& edge) const {
1.1311 + return edge == (*_matching)[_graph.u(edge)];
1.1312 }
1.1313
1.1314 /// \brief Returns the incident matching arc.
1.1315 @@ -1913,16 +1916,11 @@
1.1316 _last = _algorithm->_blossom_potential[variable].end;
1.1317 }
1.1318
1.1319 - /// \brief Invalid constructor.
1.1320 - ///
1.1321 - /// Invalid constructor.
1.1322 - BlossomIt(Invalid) : _index(-1) {}
1.1323 -
1.1324 /// \brief Conversion to node.
1.1325 ///
1.1326 /// Conversion to node.
1.1327 operator Node() const {
1.1328 - return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
1.1329 + return _algorithm->_blossom_node_list[_index];
1.1330 }
1.1331
1.1332 /// \brief Increment operator.
1.1333 @@ -1930,18 +1928,18 @@
1.1334 /// Increment operator.
1.1335 BlossomIt& operator++() {
1.1336 ++_index;
1.1337 - if (_index == _last) {
1.1338 - _index = -1;
1.1339 - }
1.1340 return *this;
1.1341 }
1.1342
1.1343 - bool operator==(const BlossomIt& it) const {
1.1344 - return _index == it._index;
1.1345 - }
1.1346 - bool operator!=(const BlossomIt& it) const {
1.1347 - return _index != it._index;
1.1348 - }
1.1349 + /// \brief Validity checking
1.1350 + ///
1.1351 + /// Checks whether the iterator is invalid.
1.1352 + bool operator==(Invalid) const { return _index == _last; }
1.1353 +
1.1354 + /// \brief Validity checking
1.1355 + ///
1.1356 + /// Checks whether the iterator is valid.
1.1357 + bool operator!=(Invalid) const { return _index != _last; }
1.1358
1.1359 private:
1.1360 const MaxWeightedMatching* _algorithm;
1.1361 @@ -1958,29 +1956,32 @@
1.1362 /// \brief Weighted perfect matching in general graphs
1.1363 ///
1.1364 /// This class provides an efficient implementation of Edmond's
1.1365 - /// maximum weighted perfecr matching algorithm. The implementation
1.1366 + /// maximum weighted perfect matching algorithm. The implementation
1.1367 /// is based on extensive use of priority queues and provides
1.1368 /// \f$O(nm\log(n))\f$ time complexity.
1.1369 ///
1.1370 /// The maximum weighted matching problem is to find undirected
1.1371 - /// arcs in the digraph with maximum overall weight and no two of
1.1372 - /// them shares their endpoints and covers all nodes. The problem
1.1373 - /// can be formulated with the next linear program:
1.1374 + /// edges in the graph with maximum overall weight and no two of
1.1375 + /// them shares their ends and covers all nodes. The problem can be
1.1376 + /// formulated with the following linear program.
1.1377 /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1.1378 - ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
1.1379 + /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
1.1380 + \quad \forall B\in\mathcal{O}\f] */
1.1381 /// \f[x_e \ge 0\quad \forall e\in E\f]
1.1382 /// \f[\max \sum_{e\in E}x_ew_e\f]
1.1383 - /// where \f$\delta(X)\f$ is the set of arcs incident to a node in
1.1384 - /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in
1.1385 - /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
1.1386 - /// the nodes.
1.1387 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1.1388 + /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
1.1389 + /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
1.1390 + /// subsets of the nodes.
1.1391 ///
1.1392 /// The algorithm calculates an optimal matching and a proof of the
1.1393 /// optimality. The solution of the dual problem can be used to check
1.1394 - /// the result of the algorithm. The dual linear problem is the next:
1.1395 - /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f]
1.1396 + /// the result of the algorithm. The dual linear problem is the
1.1397 + /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
1.1398 + w_{uv} \quad \forall uv\in E\f] */
1.1399 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
1.1400 - /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
1.1401 + /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
1.1402 + \frac{\vert B \vert - 1}{2}z_B\f] */
1.1403 ///
1.1404 /// The algorithm can be executed with \c run() or the \c init() and
1.1405 /// then the \c start() member functions. After it the matching can
1.1406 @@ -2040,16 +2041,13 @@
1.1407 int _node_num;
1.1408 int _blossom_num;
1.1409
1.1410 - typedef typename Graph::template NodeMap<int> NodeIntMap;
1.1411 - typedef typename Graph::template ArcMap<int> ArcIntMap;
1.1412 - typedef typename Graph::template EdgeMap<int> EdgeIntMap;
1.1413 typedef RangeMap<int> IntIntMap;
1.1414
1.1415 enum Status {
1.1416 EVEN = -1, MATCHED = 0, ODD = 1
1.1417 };
1.1418
1.1419 - typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
1.1420 + typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
1.1421 struct BlossomData {
1.1422 int tree;
1.1423 Status status;
1.1424 @@ -2057,21 +2055,21 @@
1.1425 Value pot, offset;
1.1426 };
1.1427
1.1428 - NodeIntMap *_blossom_index;
1.1429 + IntNodeMap *_blossom_index;
1.1430 BlossomSet *_blossom_set;
1.1431 RangeMap<BlossomData>* _blossom_data;
1.1432
1.1433 - NodeIntMap *_node_index;
1.1434 - ArcIntMap *_node_heap_index;
1.1435 + IntNodeMap *_node_index;
1.1436 + IntArcMap *_node_heap_index;
1.1437
1.1438 struct NodeData {
1.1439
1.1440 - NodeData(ArcIntMap& node_heap_index)
1.1441 + NodeData(IntArcMap& node_heap_index)
1.1442 : heap(node_heap_index) {}
1.1443
1.1444 int blossom;
1.1445 Value pot;
1.1446 - BinHeap<Value, ArcIntMap> heap;
1.1447 + BinHeap<Value, IntArcMap> heap;
1.1448 std::map<int, Arc> heap_index;
1.1449
1.1450 int tree;
1.1451 @@ -2087,8 +2085,8 @@
1.1452 IntIntMap *_delta2_index;
1.1453 BinHeap<Value, IntIntMap> *_delta2;
1.1454
1.1455 - EdgeIntMap *_delta3_index;
1.1456 - BinHeap<Value, EdgeIntMap> *_delta3;
1.1457 + IntEdgeMap *_delta3_index;
1.1458 + BinHeap<Value, IntEdgeMap> *_delta3;
1.1459
1.1460 IntIntMap *_delta4_index;
1.1461 BinHeap<Value, IntIntMap> *_delta4;
1.1462 @@ -2106,16 +2104,16 @@
1.1463 _node_potential = new NodePotential(_graph);
1.1464 }
1.1465 if (!_blossom_set) {
1.1466 - _blossom_index = new NodeIntMap(_graph);
1.1467 + _blossom_index = new IntNodeMap(_graph);
1.1468 _blossom_set = new BlossomSet(*_blossom_index);
1.1469 _blossom_data = new RangeMap<BlossomData>(_blossom_num);
1.1470 }
1.1471
1.1472 if (!_node_index) {
1.1473 - _node_index = new NodeIntMap(_graph);
1.1474 - _node_heap_index = new ArcIntMap(_graph);
1.1475 + _node_index = new IntNodeMap(_graph);
1.1476 + _node_heap_index = new IntArcMap(_graph);
1.1477 _node_data = new RangeMap<NodeData>(_node_num,
1.1478 - NodeData(*_node_heap_index));
1.1479 + NodeData(*_node_heap_index));
1.1480 }
1.1481
1.1482 if (!_tree_set) {
1.1483 @@ -2127,8 +2125,8 @@
1.1484 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
1.1485 }
1.1486 if (!_delta3) {
1.1487 - _delta3_index = new EdgeIntMap(_graph);
1.1488 - _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index);
1.1489 + _delta3_index = new IntEdgeMap(_graph);
1.1490 + _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
1.1491 }
1.1492 if (!_delta4) {
1.1493 _delta4_index = new IntIntMap(_blossom_num);
1.1494 @@ -2461,10 +2459,10 @@
1.1495 _tree_set->eraseClass(tree);
1.1496 }
1.1497
1.1498 - void augmentOnArc(const Edge& arc) {
1.1499 -
1.1500 - int left = _blossom_set->find(_graph.u(arc));
1.1501 - int right = _blossom_set->find(_graph.v(arc));
1.1502 + void augmentOnEdge(const Edge& edge) {
1.1503 +
1.1504 + int left = _blossom_set->find(_graph.u(edge));
1.1505 + int right = _blossom_set->find(_graph.v(edge));
1.1506
1.1507 int left_tree = _tree_set->find(left);
1.1508 alternatePath(left, left_tree);
1.1509 @@ -2474,8 +2472,8 @@
1.1510 alternatePath(right, right_tree);
1.1511 destroyTree(right_tree);
1.1512
1.1513 - (*_blossom_data)[left].next = _graph.direct(arc, true);
1.1514 - (*_blossom_data)[right].next = _graph.direct(arc, false);
1.1515 + (*_blossom_data)[left].next = _graph.direct(edge, true);
1.1516 + (*_blossom_data)[right].next = _graph.direct(edge, false);
1.1517 }
1.1518
1.1519 void extendOnArc(const Arc& arc) {
1.1520 @@ -2495,7 +2493,7 @@
1.1521 matchedToEven(even, tree);
1.1522 }
1.1523
1.1524 - void shrinkOnArc(const Edge& edge, int tree) {
1.1525 + void shrinkOnEdge(const Edge& edge, int tree) {
1.1526 int nca = -1;
1.1527 std::vector<int> left_path, right_path;
1.1528
1.1529 @@ -2831,7 +2829,7 @@
1.1530 createStructures();
1.1531
1.1532 for (ArcIt e(_graph); e != INVALID; ++e) {
1.1533 - _node_heap_index->set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP);
1.1534 + _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
1.1535 }
1.1536 for (EdgeIt e(_graph); e != INVALID; ++e) {
1.1537 _delta3_index->set(e, _delta3->PRE_HEAP);
1.1538 @@ -2924,9 +2922,9 @@
1.1539 int right_tree = _tree_set->find(right_blossom);
1.1540
1.1541 if (left_tree == right_tree) {
1.1542 - shrinkOnArc(e, left_tree);
1.1543 + shrinkOnEdge(e, left_tree);
1.1544 } else {
1.1545 - augmentOnArc(e);
1.1546 + augmentOnEdge(e);
1.1547 unmatched -= 2;
1.1548 }
1.1549 }
1.1550 @@ -2974,16 +2972,16 @@
1.1551 return sum /= 2;
1.1552 }
1.1553
1.1554 - /// \brief Returns true when the arc is in the matching.
1.1555 + /// \brief Returns true when the edge is in the matching.
1.1556 ///
1.1557 - /// Returns true when the arc is in the matching.
1.1558 - bool matching(const Edge& arc) const {
1.1559 - return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true);
1.1560 + /// Returns true when the edge is in the matching.
1.1561 + bool matching(const Edge& edge) const {
1.1562 + return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
1.1563 }
1.1564
1.1565 - /// \brief Returns the incident matching arc.
1.1566 + /// \brief Returns the incident matching edge.
1.1567 ///
1.1568 - /// Returns the incident matching arc from given node.
1.1569 + /// Returns the incident matching arc from given edge.
1.1570 Arc matching(const Node& node) const {
1.1571 return (*_matching)[node];
1.1572 }
1.1573 @@ -3066,16 +3064,11 @@
1.1574 _last = _algorithm->_blossom_potential[variable].end;
1.1575 }
1.1576
1.1577 - /// \brief Invalid constructor.
1.1578 - ///
1.1579 - /// Invalid constructor.
1.1580 - BlossomIt(Invalid) : _index(-1) {}
1.1581 -
1.1582 /// \brief Conversion to node.
1.1583 ///
1.1584 /// Conversion to node.
1.1585 operator Node() const {
1.1586 - return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
1.1587 + return _algorithm->_blossom_node_list[_index];
1.1588 }
1.1589
1.1590 /// \brief Increment operator.
1.1591 @@ -3083,18 +3076,18 @@
1.1592 /// Increment operator.
1.1593 BlossomIt& operator++() {
1.1594 ++_index;
1.1595 - if (_index == _last) {
1.1596 - _index = -1;
1.1597 - }
1.1598 return *this;
1.1599 }
1.1600
1.1601 - bool operator==(const BlossomIt& it) const {
1.1602 - return _index == it._index;
1.1603 - }
1.1604 - bool operator!=(const BlossomIt& it) const {
1.1605 - return _index != it._index;
1.1606 - }
1.1607 + /// \brief Validity checking
1.1608 + ///
1.1609 + /// Checks whether the iterator is invalid.
1.1610 + bool operator==(Invalid) const { return _index == _last; }
1.1611 +
1.1612 + /// \brief Validity checking
1.1613 + ///
1.1614 + /// Checks whether the iterator is valid.
1.1615 + bool operator!=(Invalid) const { return _index != _last; }
1.1616
1.1617 private:
1.1618 const MaxWeightedPerfectMatching* _algorithm;
2.1 --- a/test/CMakeLists.txt Mon Oct 13 13:56:00 2008 +0200
2.2 +++ b/test/CMakeLists.txt Mon Oct 13 14:00:11 2008 +0200
2.3 @@ -17,7 +17,6 @@
2.4 kruskal_test
2.5 maps_test
2.6 max_matching_test
2.7 - max_weighted_matching_test
2.8 random_test
2.9 path_test
2.10 time_measure_test
3.1 --- a/test/Makefile.am Mon Oct 13 13:56:00 2008 +0200
3.2 +++ b/test/Makefile.am Mon Oct 13 14:00:11 2008 +0200
3.3 @@ -20,7 +20,6 @@
3.4 test/kruskal_test \
3.5 test/maps_test \
3.6 test/max_matching_test \
3.7 - test/max_weighted_matching_test \
3.8 test/random_test \
3.9 test/path_test \
3.10 test/test_tools_fail \
3.11 @@ -45,7 +44,6 @@
3.12 test_kruskal_test_SOURCES = test/kruskal_test.cc
3.13 test_maps_test_SOURCES = test/maps_test.cc
3.14 test_max_matching_test_SOURCES = test/max_matching_test.cc
3.15 -test_max_weighted_matching_test_SOURCES = test/max_weighted_matching_test.cc
3.16 test_path_test_SOURCES = test/path_test.cc
3.17 test_random_test_SOURCES = test/random_test.cc
3.18 test_test_tools_fail_SOURCES = test/test_tools_fail.cc
4.1 --- a/test/max_matching_test.cc Mon Oct 13 13:56:00 2008 +0200
4.2 +++ b/test/max_matching_test.cc Mon Oct 13 14:00:11 2008 +0200
4.3 @@ -17,165 +17,294 @@
4.4 */
4.5
4.6 #include <iostream>
4.7 +#include <sstream>
4.8 #include <vector>
4.9 #include <queue>
4.10 #include <lemon/math.h>
4.11 #include <cstdlib>
4.12
4.13 +#include <lemon/max_matching.h>
4.14 +#include <lemon/smart_graph.h>
4.15 +#include <lemon/lgf_reader.h>
4.16 +
4.17 #include "test_tools.h"
4.18 -#include <lemon/list_graph.h>
4.19 -#include <lemon/max_matching.h>
4.20
4.21 using namespace std;
4.22 using namespace lemon;
4.23
4.24 +GRAPH_TYPEDEFS(SmartGraph);
4.25 +
4.26 +
4.27 +const int lgfn = 3;
4.28 +const std::string lgf[lgfn] = {
4.29 + "@nodes\n"
4.30 + "label\n"
4.31 + "0\n"
4.32 + "1\n"
4.33 + "2\n"
4.34 + "3\n"
4.35 + "4\n"
4.36 + "5\n"
4.37 + "6\n"
4.38 + "7\n"
4.39 + "@edges\n"
4.40 + " label weight\n"
4.41 + "7 4 0 984\n"
4.42 + "0 7 1 73\n"
4.43 + "7 1 2 204\n"
4.44 + "2 3 3 583\n"
4.45 + "2 7 4 565\n"
4.46 + "2 1 5 582\n"
4.47 + "0 4 6 551\n"
4.48 + "2 5 7 385\n"
4.49 + "1 5 8 561\n"
4.50 + "5 3 9 484\n"
4.51 + "7 5 10 904\n"
4.52 + "3 6 11 47\n"
4.53 + "7 6 12 888\n"
4.54 + "3 0 13 747\n"
4.55 + "6 1 14 310\n",
4.56 +
4.57 + "@nodes\n"
4.58 + "label\n"
4.59 + "0\n"
4.60 + "1\n"
4.61 + "2\n"
4.62 + "3\n"
4.63 + "4\n"
4.64 + "5\n"
4.65 + "6\n"
4.66 + "7\n"
4.67 + "@edges\n"
4.68 + " label weight\n"
4.69 + "2 5 0 710\n"
4.70 + "0 5 1 241\n"
4.71 + "2 4 2 856\n"
4.72 + "2 6 3 762\n"
4.73 + "4 1 4 747\n"
4.74 + "6 1 5 962\n"
4.75 + "4 7 6 723\n"
4.76 + "1 7 7 661\n"
4.77 + "2 3 8 376\n"
4.78 + "1 0 9 416\n"
4.79 + "6 7 10 391\n",
4.80 +
4.81 + "@nodes\n"
4.82 + "label\n"
4.83 + "0\n"
4.84 + "1\n"
4.85 + "2\n"
4.86 + "3\n"
4.87 + "4\n"
4.88 + "5\n"
4.89 + "6\n"
4.90 + "7\n"
4.91 + "@edges\n"
4.92 + " label weight\n"
4.93 + "6 2 0 553\n"
4.94 + "0 7 1 653\n"
4.95 + "6 3 2 22\n"
4.96 + "4 7 3 846\n"
4.97 + "7 2 4 981\n"
4.98 + "7 6 5 250\n"
4.99 + "5 2 6 539\n",
4.100 +};
4.101 +
4.102 +void checkMatching(const SmartGraph& graph,
4.103 + const MaxMatching<SmartGraph>& mm) {
4.104 + int num = 0;
4.105 +
4.106 + IntNodeMap comp_index(graph);
4.107 + UnionFind<IntNodeMap> comp(comp_index);
4.108 +
4.109 + int barrier_num = 0;
4.110 +
4.111 + for (NodeIt n(graph); n != INVALID; ++n) {
4.112 + check(mm.decomposition(n) == MaxMatching<SmartGraph>::EVEN ||
4.113 + mm.matching(n) != INVALID, "Wrong Gallai-Edmonds decomposition");
4.114 + if (mm.decomposition(n) == MaxMatching<SmartGraph>::ODD) {
4.115 + ++barrier_num;
4.116 + } else {
4.117 + comp.insert(n);
4.118 + }
4.119 + }
4.120 +
4.121 + for (EdgeIt e(graph); e != INVALID; ++e) {
4.122 + if (mm.matching(e)) {
4.123 + check(e == mm.matching(graph.u(e)), "Wrong matching");
4.124 + check(e == mm.matching(graph.v(e)), "Wrong matching");
4.125 + ++num;
4.126 + }
4.127 + check(mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::EVEN ||
4.128 + mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::MATCHED,
4.129 + "Wrong Gallai-Edmonds decomposition");
4.130 +
4.131 + check(mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::EVEN ||
4.132 + mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::MATCHED,
4.133 + "Wrong Gallai-Edmonds decomposition");
4.134 +
4.135 + if (mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::ODD &&
4.136 + mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::ODD) {
4.137 + comp.join(graph.u(e), graph.v(e));
4.138 + }
4.139 + }
4.140 +
4.141 + std::set<int> comp_root;
4.142 + int odd_comp_num = 0;
4.143 + for (NodeIt n(graph); n != INVALID; ++n) {
4.144 + if (mm.decomposition(n) != MaxMatching<SmartGraph>::ODD) {
4.145 + int root = comp.find(n);
4.146 + if (comp_root.find(root) == comp_root.end()) {
4.147 + comp_root.insert(root);
4.148 + if (comp.size(n) % 2 == 1) {
4.149 + ++odd_comp_num;
4.150 + }
4.151 + }
4.152 + }
4.153 + }
4.154 +
4.155 + check(mm.matchingSize() == num, "Wrong matching");
4.156 + check(2 * num == countNodes(graph) - (odd_comp_num - barrier_num),
4.157 + "Wrong matching");
4.158 + return;
4.159 +}
4.160 +
4.161 +void checkWeightedMatching(const SmartGraph& graph,
4.162 + const SmartGraph::EdgeMap<int>& weight,
4.163 + const MaxWeightedMatching<SmartGraph>& mwm) {
4.164 + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
4.165 + if (graph.u(e) == graph.v(e)) continue;
4.166 + int rw = mwm.nodeValue(graph.u(e)) + mwm.nodeValue(graph.v(e));
4.167 +
4.168 + for (int i = 0; i < mwm.blossomNum(); ++i) {
4.169 + bool s = false, t = false;
4.170 + for (MaxWeightedMatching<SmartGraph>::BlossomIt n(mwm, i);
4.171 + n != INVALID; ++n) {
4.172 + if (graph.u(e) == n) s = true;
4.173 + if (graph.v(e) == n) t = true;
4.174 + }
4.175 + if (s == true && t == true) {
4.176 + rw += mwm.blossomValue(i);
4.177 + }
4.178 + }
4.179 + rw -= weight[e] * mwm.dualScale;
4.180 +
4.181 + check(rw >= 0, "Negative reduced weight");
4.182 + check(rw == 0 || !mwm.matching(e),
4.183 + "Non-zero reduced weight on matching edge");
4.184 + }
4.185 +
4.186 + int pv = 0;
4.187 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
4.188 + if (mwm.matching(n) != INVALID) {
4.189 + check(mwm.nodeValue(n) >= 0, "Invalid node value");
4.190 + pv += weight[mwm.matching(n)];
4.191 + SmartGraph::Node o = graph.target(mwm.matching(n));
4.192 + check(mwm.mate(n) == o, "Invalid matching");
4.193 + check(mwm.matching(n) == graph.oppositeArc(mwm.matching(o)),
4.194 + "Invalid matching");
4.195 + } else {
4.196 + check(mwm.mate(n) == INVALID, "Invalid matching");
4.197 + check(mwm.nodeValue(n) == 0, "Invalid matching");
4.198 + }
4.199 + }
4.200 +
4.201 + int dv = 0;
4.202 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
4.203 + dv += mwm.nodeValue(n);
4.204 + }
4.205 +
4.206 + for (int i = 0; i < mwm.blossomNum(); ++i) {
4.207 + check(mwm.blossomValue(i) >= 0, "Invalid blossom value");
4.208 + check(mwm.blossomSize(i) % 2 == 1, "Even blossom size");
4.209 + dv += mwm.blossomValue(i) * ((mwm.blossomSize(i) - 1) / 2);
4.210 + }
4.211 +
4.212 + check(pv * mwm.dualScale == dv * 2, "Wrong duality");
4.213 +
4.214 + return;
4.215 +}
4.216 +
4.217 +void checkWeightedPerfectMatching(const SmartGraph& graph,
4.218 + const SmartGraph::EdgeMap<int>& weight,
4.219 + const MaxWeightedPerfectMatching<SmartGraph>& mwpm) {
4.220 + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
4.221 + if (graph.u(e) == graph.v(e)) continue;
4.222 + int rw = mwpm.nodeValue(graph.u(e)) + mwpm.nodeValue(graph.v(e));
4.223 +
4.224 + for (int i = 0; i < mwpm.blossomNum(); ++i) {
4.225 + bool s = false, t = false;
4.226 + for (MaxWeightedPerfectMatching<SmartGraph>::BlossomIt n(mwpm, i);
4.227 + n != INVALID; ++n) {
4.228 + if (graph.u(e) == n) s = true;
4.229 + if (graph.v(e) == n) t = true;
4.230 + }
4.231 + if (s == true && t == true) {
4.232 + rw += mwpm.blossomValue(i);
4.233 + }
4.234 + }
4.235 + rw -= weight[e] * mwpm.dualScale;
4.236 +
4.237 + check(rw >= 0, "Negative reduced weight");
4.238 + check(rw == 0 || !mwpm.matching(e),
4.239 + "Non-zero reduced weight on matching edge");
4.240 + }
4.241 +
4.242 + int pv = 0;
4.243 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
4.244 + check(mwpm.matching(n) != INVALID, "Non perfect");
4.245 + pv += weight[mwpm.matching(n)];
4.246 + SmartGraph::Node o = graph.target(mwpm.matching(n));
4.247 + check(mwpm.mate(n) == o, "Invalid matching");
4.248 + check(mwpm.matching(n) == graph.oppositeArc(mwpm.matching(o)),
4.249 + "Invalid matching");
4.250 + }
4.251 +
4.252 + int dv = 0;
4.253 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
4.254 + dv += mwpm.nodeValue(n);
4.255 + }
4.256 +
4.257 + for (int i = 0; i < mwpm.blossomNum(); ++i) {
4.258 + check(mwpm.blossomValue(i) >= 0, "Invalid blossom value");
4.259 + check(mwpm.blossomSize(i) % 2 == 1, "Even blossom size");
4.260 + dv += mwpm.blossomValue(i) * ((mwpm.blossomSize(i) - 1) / 2);
4.261 + }
4.262 +
4.263 + check(pv * mwpm.dualScale == dv * 2, "Wrong duality");
4.264 +
4.265 + return;
4.266 +}
4.267 +
4.268 +
4.269 int main() {
4.270
4.271 - typedef ListGraph Graph;
4.272 + for (int i = 0; i < lgfn; ++i) {
4.273 + SmartGraph graph;
4.274 + SmartGraph::EdgeMap<int> weight(graph);
4.275
4.276 - GRAPH_TYPEDEFS(Graph);
4.277 + istringstream lgfs(lgf[i]);
4.278 + graphReader(graph, lgfs).
4.279 + edgeMap("weight", weight).run();
4.280
4.281 - Graph g;
4.282 - g.clear();
4.283 + MaxMatching<SmartGraph> mm(graph);
4.284 + mm.run();
4.285 + checkMatching(graph, mm);
4.286
4.287 - std::vector<Graph::Node> nodes;
4.288 - for (int i=0; i<13; ++i)
4.289 - nodes.push_back(g.addNode());
4.290 + MaxWeightedMatching<SmartGraph> mwm(graph, weight);
4.291 + mwm.run();
4.292 + checkWeightedMatching(graph, weight, mwm);
4.293
4.294 - g.addEdge(nodes[0], nodes[0]);
4.295 - g.addEdge(nodes[6], nodes[10]);
4.296 - g.addEdge(nodes[5], nodes[10]);
4.297 - g.addEdge(nodes[4], nodes[10]);
4.298 - g.addEdge(nodes[3], nodes[11]);
4.299 - g.addEdge(nodes[1], nodes[6]);
4.300 - g.addEdge(nodes[4], nodes[7]);
4.301 - g.addEdge(nodes[1], nodes[8]);
4.302 - g.addEdge(nodes[0], nodes[8]);
4.303 - g.addEdge(nodes[3], nodes[12]);
4.304 - g.addEdge(nodes[6], nodes[9]);
4.305 - g.addEdge(nodes[9], nodes[11]);
4.306 - g.addEdge(nodes[2], nodes[10]);
4.307 - g.addEdge(nodes[10], nodes[8]);
4.308 - g.addEdge(nodes[5], nodes[8]);
4.309 - g.addEdge(nodes[6], nodes[3]);
4.310 - g.addEdge(nodes[0], nodes[5]);
4.311 - g.addEdge(nodes[6], nodes[12]);
4.312 + MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
4.313 + bool perfect = mwpm.run();
4.314
4.315 - MaxMatching<Graph> max_matching(g);
4.316 - max_matching.init();
4.317 - max_matching.startDense();
4.318 + check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
4.319 + "Perfect matching found");
4.320
4.321 - int s=0;
4.322 - Graph::NodeMap<Node> mate(g,INVALID);
4.323 - max_matching.mateMap(mate);
4.324 - for(NodeIt v(g); v!=INVALID; ++v) {
4.325 - if ( mate[v]!=INVALID ) ++s;
4.326 - }
4.327 - int size=int(s/2); //size will be used as the size of a maxmatching
4.328 -
4.329 - for(NodeIt v(g); v!=INVALID; ++v) {
4.330 - max_matching.mate(v);
4.331 - }
4.332 -
4.333 - check ( size == max_matching.size(), "mate() returns a different size matching than max_matching.size()" );
4.334 -
4.335 - Graph::NodeMap<MaxMatching<Graph>::DecompType> pos0(g);
4.336 - max_matching.decomposition(pos0);
4.337 -
4.338 - max_matching.init();
4.339 - max_matching.startSparse();
4.340 - s=0;
4.341 - max_matching.mateMap(mate);
4.342 - for(NodeIt v(g); v!=INVALID; ++v) {
4.343 - if ( mate[v]!=INVALID ) ++s;
4.344 - }
4.345 - check ( int(s/2) == size, "The size does not equal!" );
4.346 -
4.347 - Graph::NodeMap<MaxMatching<Graph>::DecompType> pos1(g);
4.348 - max_matching.decomposition(pos1);
4.349 -
4.350 - max_matching.run();
4.351 - s=0;
4.352 - max_matching.mateMap(mate);
4.353 - for(NodeIt v(g); v!=INVALID; ++v) {
4.354 - if ( mate[v]!=INVALID ) ++s;
4.355 - }
4.356 - check ( int(s/2) == size, "The size does not equal!" );
4.357 -
4.358 - Graph::NodeMap<MaxMatching<Graph>::DecompType> pos2(g);
4.359 - max_matching.decomposition(pos2);
4.360 -
4.361 - max_matching.run();
4.362 - s=0;
4.363 - max_matching.mateMap(mate);
4.364 - for(NodeIt v(g); v!=INVALID; ++v) {
4.365 - if ( mate[v]!=INVALID ) ++s;
4.366 - }
4.367 - check ( int(s/2) == size, "The size does not equal!" );
4.368 -
4.369 - Graph::NodeMap<MaxMatching<Graph>::DecompType> pos(g);
4.370 - max_matching.decomposition(pos);
4.371 -
4.372 - bool ismatching=true;
4.373 - for(NodeIt v(g); v!=INVALID; ++v) {
4.374 - if ( mate[v]!=INVALID ) {
4.375 - Node u=mate[v];
4.376 - if (mate[u]!=v) ismatching=false;
4.377 + if (perfect) {
4.378 + checkWeightedPerfectMatching(graph, weight, mwpm);
4.379 }
4.380 }
4.381 - check ( ismatching, "It is not a matching!" );
4.382 -
4.383 - bool coincide=true;
4.384 - for(NodeIt v(g); v!=INVALID; ++v) {
4.385 - if ( pos0[v] != pos1[v] || pos1[v]!=pos2[v] || pos2[v]!=pos[v] ) {
4.386 - coincide=false;
4.387 - }
4.388 - }
4.389 - check ( coincide, "The decompositions do not coincide! " );
4.390 -
4.391 - bool noarc=true;
4.392 - for(EdgeIt e(g); e!=INVALID; ++e) {
4.393 - if ( (pos[g.v(e)]==max_matching.C &&
4.394 - pos[g.u(e)]==max_matching.D) ||
4.395 - (pos[g.v(e)]==max_matching.D &&
4.396 - pos[g.u(e)]==max_matching.C) )
4.397 - noarc=false;
4.398 - }
4.399 - check ( noarc, "There are arcs between D and C!" );
4.400 -
4.401 - bool oddcomp=true;
4.402 - Graph::NodeMap<bool> todo(g,true);
4.403 - int num_comp=0;
4.404 - for(NodeIt v(g); v!=INVALID; ++v) {
4.405 - if ( pos[v]==max_matching.D && todo[v] ) {
4.406 - int comp_size=1;
4.407 - ++num_comp;
4.408 - std::queue<Node> Q;
4.409 - Q.push(v);
4.410 - todo.set(v,false);
4.411 - while (!Q.empty()) {
4.412 - Node w=Q.front();
4.413 - Q.pop();
4.414 - for(IncEdgeIt e(g,w); e!=INVALID; ++e) {
4.415 - Node u=g.runningNode(e);
4.416 - if ( pos[u]==max_matching.D && todo[u] ) {
4.417 - ++comp_size;
4.418 - Q.push(u);
4.419 - todo.set(u,false);
4.420 - }
4.421 - }
4.422 - }
4.423 - if ( !(comp_size % 2) ) oddcomp=false;
4.424 - }
4.425 - }
4.426 - check ( oddcomp, "A component of g[D] is not odd." );
4.427 -
4.428 - int barrier=0;
4.429 - for(NodeIt v(g); v!=INVALID; ++v) {
4.430 - if ( pos[v]==max_matching.A ) ++barrier;
4.431 - }
4.432 - int expected_size=int( countNodes(g)-num_comp+barrier)/2;
4.433 - check ( size==expected_size, "The size of the matching is wrong." );
4.434
4.435 return 0;
4.436 }
5.1 --- a/test/max_weighted_matching_test.cc Mon Oct 13 13:56:00 2008 +0200
5.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
5.3 @@ -1,250 +0,0 @@
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 <sstream>
5.24 -#include <vector>
5.25 -#include <queue>
5.26 -#include <lemon/math.h>
5.27 -#include <cstdlib>
5.28 -
5.29 -#include "test_tools.h"
5.30 -#include <lemon/smart_graph.h>
5.31 -#include <lemon/max_matching.h>
5.32 -#include <lemon/lgf_reader.h>
5.33 -
5.34 -using namespace std;
5.35 -using namespace lemon;
5.36 -
5.37 -GRAPH_TYPEDEFS(SmartGraph);
5.38 -
5.39 -const int lgfn = 3;
5.40 -const std::string lgf[lgfn] = {
5.41 - "@nodes\n"
5.42 - "label\n"
5.43 - "0\n"
5.44 - "1\n"
5.45 - "2\n"
5.46 - "3\n"
5.47 - "4\n"
5.48 - "5\n"
5.49 - "6\n"
5.50 - "7\n"
5.51 - "@edges\n"
5.52 - "label weight\n"
5.53 - "7 4 0 984\n"
5.54 - "0 7 1 73\n"
5.55 - "7 1 2 204\n"
5.56 - "2 3 3 583\n"
5.57 - "2 7 4 565\n"
5.58 - "2 1 5 582\n"
5.59 - "0 4 6 551\n"
5.60 - "2 5 7 385\n"
5.61 - "1 5 8 561\n"
5.62 - "5 3 9 484\n"
5.63 - "7 5 10 904\n"
5.64 - "3 6 11 47\n"
5.65 - "7 6 12 888\n"
5.66 - "3 0 13 747\n"
5.67 - "6 1 14 310\n",
5.68 -
5.69 - "@nodes\n"
5.70 - "label\n"
5.71 - "0\n"
5.72 - "1\n"
5.73 - "2\n"
5.74 - "3\n"
5.75 - "4\n"
5.76 - "5\n"
5.77 - "6\n"
5.78 - "7\n"
5.79 - "@edges\n"
5.80 - "label weight\n"
5.81 - "2 5 0 710\n"
5.82 - "0 5 1 241\n"
5.83 - "2 4 2 856\n"
5.84 - "2 6 3 762\n"
5.85 - "4 1 4 747\n"
5.86 - "6 1 5 962\n"
5.87 - "4 7 6 723\n"
5.88 - "1 7 7 661\n"
5.89 - "2 3 8 376\n"
5.90 - "1 0 9 416\n"
5.91 - "6 7 10 391\n",
5.92 -
5.93 - "@nodes\n"
5.94 - "label\n"
5.95 - "0\n"
5.96 - "1\n"
5.97 - "2\n"
5.98 - "3\n"
5.99 - "4\n"
5.100 - "5\n"
5.101 - "6\n"
5.102 - "7\n"
5.103 - "@edges\n"
5.104 - "label weight\n"
5.105 - "6 2 0 553\n"
5.106 - "0 7 1 653\n"
5.107 - "6 3 2 22\n"
5.108 - "4 7 3 846\n"
5.109 - "7 2 4 981\n"
5.110 - "7 6 5 250\n"
5.111 - "5 2 6 539\n"
5.112 -};
5.113 -
5.114 -void checkMatching(const SmartGraph& graph,
5.115 - const SmartGraph::EdgeMap<int>& weight,
5.116 - const MaxWeightedMatching<SmartGraph>& mwm) {
5.117 - for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
5.118 - if (graph.u(e) == graph.v(e)) continue;
5.119 - int rw =
5.120 - mwm.nodeValue(graph.u(e)) + mwm.nodeValue(graph.v(e));
5.121 -
5.122 - for (int i = 0; i < mwm.blossomNum(); ++i) {
5.123 - bool s = false, t = false;
5.124 - for (MaxWeightedMatching<SmartGraph>::BlossomIt n(mwm, i);
5.125 - n != INVALID; ++n) {
5.126 - if (graph.u(e) == n) s = true;
5.127 - if (graph.v(e) == n) t = true;
5.128 - }
5.129 - if (s == true && t == true) {
5.130 - rw += mwm.blossomValue(i);
5.131 - }
5.132 - }
5.133 - rw -= weight[e] * mwm.dualScale;
5.134 -
5.135 - check(rw >= 0, "Negative reduced weight");
5.136 - check(rw == 0 || !mwm.matching(e),
5.137 - "Non-zero reduced weight on matching arc");
5.138 - }
5.139 -
5.140 - int pv = 0;
5.141 - for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
5.142 - if (mwm.matching(n) != INVALID) {
5.143 - check(mwm.nodeValue(n) >= 0, "Invalid node value");
5.144 - pv += weight[mwm.matching(n)];
5.145 - SmartGraph::Node o = graph.target(mwm.matching(n));
5.146 - check(mwm.mate(n) == o, "Invalid matching");
5.147 - check(mwm.matching(n) == graph.oppositeArc(mwm.matching(o)),
5.148 - "Invalid matching");
5.149 - } else {
5.150 - check(mwm.mate(n) == INVALID, "Invalid matching");
5.151 - check(mwm.nodeValue(n) == 0, "Invalid matching");
5.152 - }
5.153 - }
5.154 -
5.155 - int dv = 0;
5.156 - for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
5.157 - dv += mwm.nodeValue(n);
5.158 - }
5.159 -
5.160 - for (int i = 0; i < mwm.blossomNum(); ++i) {
5.161 - check(mwm.blossomValue(i) >= 0, "Invalid blossom value");
5.162 - check(mwm.blossomSize(i) % 2 == 1, "Even blossom size");
5.163 - dv += mwm.blossomValue(i) * ((mwm.blossomSize(i) - 1) / 2);
5.164 - }
5.165 -
5.166 - check(pv * mwm.dualScale == dv * 2, "Wrong duality");
5.167 -
5.168 - return;
5.169 -}
5.170 -
5.171 -void checkPerfectMatching(const SmartGraph& graph,
5.172 - const SmartGraph::EdgeMap<int>& weight,
5.173 - const MaxWeightedPerfectMatching<SmartGraph>& mwpm) {
5.174 - for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
5.175 - if (graph.u(e) == graph.v(e)) continue;
5.176 - int rw =
5.177 - mwpm.nodeValue(graph.u(e)) + mwpm.nodeValue(graph.v(e));
5.178 -
5.179 - for (int i = 0; i < mwpm.blossomNum(); ++i) {
5.180 - bool s = false, t = false;
5.181 - for (MaxWeightedPerfectMatching<SmartGraph>::BlossomIt n(mwpm, i);
5.182 - n != INVALID; ++n) {
5.183 - if (graph.u(e) == n) s = true;
5.184 - if (graph.v(e) == n) t = true;
5.185 - }
5.186 - if (s == true && t == true) {
5.187 - rw += mwpm.blossomValue(i);
5.188 - }
5.189 - }
5.190 - rw -= weight[e] * mwpm.dualScale;
5.191 -
5.192 - check(rw >= 0, "Negative reduced weight");
5.193 - check(rw == 0 || !mwpm.matching(e),
5.194 - "Non-zero reduced weight on matching arc");
5.195 - }
5.196 -
5.197 - int pv = 0;
5.198 - for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
5.199 - check(mwpm.matching(n) != INVALID, "Non perfect");
5.200 - pv += weight[mwpm.matching(n)];
5.201 - SmartGraph::Node o = graph.target(mwpm.matching(n));
5.202 - check(mwpm.mate(n) == o, "Invalid matching");
5.203 - check(mwpm.matching(n) == graph.oppositeArc(mwpm.matching(o)),
5.204 - "Invalid matching");
5.205 - }
5.206 -
5.207 - int dv = 0;
5.208 - for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
5.209 - dv += mwpm.nodeValue(n);
5.210 - }
5.211 -
5.212 - for (int i = 0; i < mwpm.blossomNum(); ++i) {
5.213 - check(mwpm.blossomValue(i) >= 0, "Invalid blossom value");
5.214 - check(mwpm.blossomSize(i) % 2 == 1, "Even blossom size");
5.215 - dv += mwpm.blossomValue(i) * ((mwpm.blossomSize(i) - 1) / 2);
5.216 - }
5.217 -
5.218 - check(pv * mwpm.dualScale == dv * 2, "Wrong duality");
5.219 -
5.220 - return;
5.221 -}
5.222 -
5.223 -
5.224 -int main() {
5.225 -
5.226 - for (int i = 0; i < lgfn; ++i) {
5.227 - SmartGraph graph;
5.228 - SmartGraph::EdgeMap<int> weight(graph);
5.229 -
5.230 - istringstream lgfs(lgf[i]);
5.231 - graphReader(graph, lgfs).
5.232 - edgeMap("weight", weight).run();
5.233 -
5.234 - MaxWeightedMatching<SmartGraph> mwm(graph, weight);
5.235 - mwm.run();
5.236 - checkMatching(graph, weight, mwm);
5.237 -
5.238 - MaxMatching<SmartGraph> mm(graph);
5.239 - mm.run();
5.240 -
5.241 - MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
5.242 - bool perfect = mwpm.run();
5.243 -
5.244 - check(perfect == (mm.size() * 2 == countNodes(graph)),
5.245 - "Perfect matching found");
5.246 -
5.247 - if (perfect) {
5.248 - checkPerfectMatching(graph, weight, mwpm);
5.249 - }
5.250 - }
5.251 -
5.252 - return 0;
5.253 -}