1 /* -*- mode: C++; indent-tabs-mode: nil; -*- 2 * 3 * This file is a part of LEMON, a generic C++ optimization library. 4 * 5 * Copyright (C) 2003-2013 6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport 7 * (Egervary Research Group on Combinatorial Optimization, EGRES). 8 * 9 * Permission to use, modify and distribute this software is granted 10 * provided that this copyright notice appears in all copies. For 11 * precise terms see the accompanying LICENSE file. 12 * 13 * This software is provided "AS IS" with no warranty of any kind, 14 * express or implied, and with no claim as to its suitability for any 15 * purpose. 16 * 17 */ 18 19 #ifndef LEMON_MATCHING_H 20 #define LEMON_MATCHING_H 21 22 #include <vector> 23 #include <queue> 24 #include <set> 25 #include <limits> 26 27 #include <lemon/core.h> 28 #include <lemon/unionfind.h> 29 #include <lemon/bin_heap.h> 30 #include <lemon/maps.h> 31 #include <lemon/fractional_matching.h> 32 33 ///\ingroup matching 34 ///\file 35 ///\brief Maximum matching algorithms in general graphs. 36 37 namespace lemon { 38 39 /// \ingroup matching 40 /// 41 /// \brief Maximum cardinality matching in general graphs 42 /// 43 /// This class implements Edmonds' alternating forest matching algorithm 44 /// for finding a maximum cardinality matching in a general undirected graph. 45 /// It can be started from an arbitrary initial matching 46 /// (the default is the empty one). 47 /// 48 /// The dual solution of the problem is a map of the nodes to 49 /// \ref MaxMatching::Status "Status", having values \c EVEN (or \c D), 50 /// \c ODD (or \c A) and \c MATCHED (or \c C) defining the Gallai-Edmonds 51 /// decomposition of the graph. The nodes in \c EVEN/D induce a subgraph 52 /// with factor-critical components, the nodes in \c ODD/A form the 53 /// canonical barrier, and the nodes in \c MATCHED/C induce a graph having 54 /// a perfect matching. The number of the factor-critical components 55 /// minus the number of barrier nodes is a lower bound on the 56 /// unmatched nodes, and the matching is optimal if and only if this bound is 57 /// tight. This decomposition can be obtained using \ref status() or 58 /// \ref statusMap() after running the algorithm. 59 /// 60 /// \tparam GR The undirected graph type the algorithm runs on. 61 template <typename GR> 62 class MaxMatching { 63 public: 64 65 /// The graph type of the algorithm 66 typedef GR Graph; 67 /// The type of the matching map 68 typedef typename Graph::template NodeMap<typename Graph::Arc> 69 MatchingMap; 70 71 ///\brief Status constants for Gallai-Edmonds decomposition. 72 /// 73 ///These constants are used for indicating the Gallai-Edmonds 74 ///decomposition of a graph. The nodes with status \c EVEN (or \c D) 75 ///induce a subgraph with factor-critical components, the nodes with 76 ///status \c ODD (or \c A) form the canonical barrier, and the nodes 77 ///with status \c MATCHED (or \c C) induce a subgraph having a 78 ///perfect matching. 79 enum Status { 80 EVEN = 1, ///< = 1. (\c D is an alias for \c EVEN.) 81 D = 1, 82 MATCHED = 0, ///< = 0. (\c C is an alias for \c MATCHED.) 83 C = 0, 84 ODD = -1, ///< = -1. (\c A is an alias for \c ODD.) 85 A = -1, 86 UNMATCHED = -2 ///< = -2. 87 }; 88 89 /// The type of the status map 90 typedef typename Graph::template NodeMap<Status> StatusMap; 91 92 private: 93 94 TEMPLATE_GRAPH_TYPEDEFS(Graph); 95 96 typedef UnionFindEnum<IntNodeMap> BlossomSet; 97 typedef ExtendFindEnum<IntNodeMap> TreeSet; 98 typedef RangeMap<Node> NodeIntMap; 99 typedef MatchingMap EarMap; 100 typedef std::vector<Node> NodeQueue; 101 102 const Graph& _graph; 103 MatchingMap* _matching; 104 StatusMap* _status; 105 106 EarMap* _ear; 107 108 IntNodeMap* _blossom_set_index; 109 BlossomSet* _blossom_set; 110 NodeIntMap* _blossom_rep; 111 112 IntNodeMap* _tree_set_index; 113 TreeSet* _tree_set; 114 115 NodeQueue _node_queue; 116 int _process, _postpone, _last; 117 118 int _node_num; 119 120 private: 121 createStructures()122 void createStructures() { 123 _node_num = countNodes(_graph); 124 if (!_matching) { 125 _matching = new MatchingMap(_graph); 126 } 127 if (!_status) { 128 _status = new StatusMap(_graph); 129 } 130 if (!_ear) { 131 _ear = new EarMap(_graph); 132 } 133 if (!_blossom_set) { 134 _blossom_set_index = new IntNodeMap(_graph); 135 _blossom_set = new BlossomSet(*_blossom_set_index); 136 } 137 if (!_blossom_rep) { 138 _blossom_rep = new NodeIntMap(_node_num); 139 } 140 if (!_tree_set) { 141 _tree_set_index = new IntNodeMap(_graph); 142 _tree_set = new TreeSet(*_tree_set_index); 143 } 144 _node_queue.resize(_node_num); 145 } 146 destroyStructures()147 void destroyStructures() { 148 if (_matching) { 149 delete _matching; 150 } 151 if (_status) { 152 delete _status; 153 } 154 if (_ear) { 155 delete _ear; 156 } 157 if (_blossom_set) { 158 delete _blossom_set; 159 delete _blossom_set_index; 160 } 161 if (_blossom_rep) { 162 delete _blossom_rep; 163 } 164 if (_tree_set) { 165 delete _tree_set_index; 166 delete _tree_set; 167 } 168 } 169 processDense(const Node & n)170 void processDense(const Node& n) { 171 _process = _postpone = _last = 0; 172 _node_queue[_last++] = n; 173 174 while (_process != _last) { 175 Node u = _node_queue[_process++]; 176 for (OutArcIt a(_graph, u); a != INVALID; ++a) { 177 Node v = _graph.target(a); 178 if ((*_status)[v] == MATCHED) { 179 extendOnArc(a); 180 } else if ((*_status)[v] == UNMATCHED) { 181 augmentOnArc(a); 182 return; 183 } 184 } 185 } 186 187 while (_postpone != _last) { 188 Node u = _node_queue[_postpone++]; 189 190 for (OutArcIt a(_graph, u); a != INVALID ; ++a) { 191 Node v = _graph.target(a); 192 193 if ((*_status)[v] == EVEN) { 194 if (_blossom_set->find(u) != _blossom_set->find(v)) { 195 shrinkOnEdge(a); 196 } 197 } 198 199 while (_process != _last) { 200 Node w = _node_queue[_process++]; 201 for (OutArcIt b(_graph, w); b != INVALID; ++b) { 202 Node x = _graph.target(b); 203 if ((*_status)[x] == MATCHED) { 204 extendOnArc(b); 205 } else if ((*_status)[x] == UNMATCHED) { 206 augmentOnArc(b); 207 return; 208 } 209 } 210 } 211 } 212 } 213 } 214 processSparse(const Node & n)215 void processSparse(const Node& n) { 216 _process = _last = 0; 217 _node_queue[_last++] = n; 218 while (_process != _last) { 219 Node u = _node_queue[_process++]; 220 for (OutArcIt a(_graph, u); a != INVALID; ++a) { 221 Node v = _graph.target(a); 222 223 if ((*_status)[v] == EVEN) { 224 if (_blossom_set->find(u) != _blossom_set->find(v)) { 225 shrinkOnEdge(a); 226 } 227 } else if ((*_status)[v] == MATCHED) { 228 extendOnArc(a); 229 } else if ((*_status)[v] == UNMATCHED) { 230 augmentOnArc(a); 231 return; 232 } 233 } 234 } 235 } 236 shrinkOnEdge(const Edge & e)237 void shrinkOnEdge(const Edge& e) { 238 Node nca = INVALID; 239 240 { 241 std::set<Node> left_set, right_set; 242 243 Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))]; 244 left_set.insert(left); 245 246 Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))]; 247 right_set.insert(right); 248 249 while (true) { 250 if ((*_matching)[left] == INVALID) break; 251 left = _graph.target((*_matching)[left]); 252 left = (*_blossom_rep)[_blossom_set-> 253 find(_graph.target((*_ear)[left]))]; 254 if (right_set.find(left) != right_set.end()) { 255 nca = left; 256 break; 257 } 258 left_set.insert(left); 259 260 if ((*_matching)[right] == INVALID) break; 261 right = _graph.target((*_matching)[right]); 262 right = (*_blossom_rep)[_blossom_set-> 263 find(_graph.target((*_ear)[right]))]; 264 if (left_set.find(right) != left_set.end()) { 265 nca = right; 266 break; 267 } 268 right_set.insert(right); 269 } 270 271 if (nca == INVALID) { 272 if ((*_matching)[left] == INVALID) { 273 nca = right; 274 while (left_set.find(nca) == left_set.end()) { 275 nca = _graph.target((*_matching)[nca]); 276 nca =(*_blossom_rep)[_blossom_set-> 277 find(_graph.target((*_ear)[nca]))]; 278 } 279 } else { 280 nca = left; 281 while (right_set.find(nca) == right_set.end()) { 282 nca = _graph.target((*_matching)[nca]); 283 nca = (*_blossom_rep)[_blossom_set-> 284 find(_graph.target((*_ear)[nca]))]; 285 } 286 } 287 } 288 } 289 290 { 291 292 Node node = _graph.u(e); 293 Arc arc = _graph.direct(e, true); 294 Node base = (*_blossom_rep)[_blossom_set->find(node)]; 295 296 while (base != nca) { 297 (*_ear)[node] = arc; 298 299 Node n = node; 300 while (n != base) { 301 n = _graph.target((*_matching)[n]); 302 Arc a = (*_ear)[n]; 303 n = _graph.target(a); 304 (*_ear)[n] = _graph.oppositeArc(a); 305 } 306 node = _graph.target((*_matching)[base]); 307 _tree_set->erase(base); 308 _tree_set->erase(node); 309 _blossom_set->insert(node, _blossom_set->find(base)); 310 (*_status)[node] = EVEN; 311 _node_queue[_last++] = node; 312 arc = _graph.oppositeArc((*_ear)[node]); 313 node = _graph.target((*_ear)[node]); 314 base = (*_blossom_rep)[_blossom_set->find(node)]; 315 _blossom_set->join(_graph.target(arc), base); 316 } 317 } 318 319 (*_blossom_rep)[_blossom_set->find(nca)] = nca; 320 321 { 322 323 Node node = _graph.v(e); 324 Arc arc = _graph.direct(e, false); 325 Node base = (*_blossom_rep)[_blossom_set->find(node)]; 326 327 while (base != nca) { 328 (*_ear)[node] = arc; 329 330 Node n = node; 331 while (n != base) { 332 n = _graph.target((*_matching)[n]); 333 Arc a = (*_ear)[n]; 334 n = _graph.target(a); 335 (*_ear)[n] = _graph.oppositeArc(a); 336 } 337 node = _graph.target((*_matching)[base]); 338 _tree_set->erase(base); 339 _tree_set->erase(node); 340 _blossom_set->insert(node, _blossom_set->find(base)); 341 (*_status)[node] = EVEN; 342 _node_queue[_last++] = node; 343 arc = _graph.oppositeArc((*_ear)[node]); 344 node = _graph.target((*_ear)[node]); 345 base = (*_blossom_rep)[_blossom_set->find(node)]; 346 _blossom_set->join(_graph.target(arc), base); 347 } 348 } 349 350 (*_blossom_rep)[_blossom_set->find(nca)] = nca; 351 } 352 extendOnArc(const Arc & a)353 void extendOnArc(const Arc& a) { 354 Node base = _graph.source(a); 355 Node odd = _graph.target(a); 356 357 (*_ear)[odd] = _graph.oppositeArc(a); 358 Node even = _graph.target((*_matching)[odd]); 359 (*_blossom_rep)[_blossom_set->insert(even)] = even; 360 (*_status)[odd] = ODD; 361 (*_status)[even] = EVEN; 362 int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]); 363 _tree_set->insert(odd, tree); 364 _tree_set->insert(even, tree); 365 _node_queue[_last++] = even; 366 367 } 368 augmentOnArc(const Arc & a)369 void augmentOnArc(const Arc& a) { 370 Node even = _graph.source(a); 371 Node odd = _graph.target(a); 372 373 int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]); 374 375 (*_matching)[odd] = _graph.oppositeArc(a); 376 (*_status)[odd] = MATCHED; 377 378 Arc arc = (*_matching)[even]; 379 (*_matching)[even] = a; 380 381 while (arc != INVALID) { 382 odd = _graph.target(arc); 383 arc = (*_ear)[odd]; 384 even = _graph.target(arc); 385 (*_matching)[odd] = arc; 386 arc = (*_matching)[even]; 387 (*_matching)[even] = _graph.oppositeArc((*_matching)[odd]); 388 } 389 390 for (typename TreeSet::ItemIt it(*_tree_set, tree); 391 it != INVALID; ++it) { 392 if ((*_status)[it] == ODD) { 393 (*_status)[it] = MATCHED; 394 } else { 395 int blossom = _blossom_set->find(it); 396 for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom); 397 jt != INVALID; ++jt) { 398 (*_status)[jt] = MATCHED; 399 } 400 _blossom_set->eraseClass(blossom); 401 } 402 } 403 _tree_set->eraseClass(tree); 404 405 } 406 407 public: 408 409 /// \brief Constructor 410 /// 411 /// Constructor. MaxMatching(const Graph & graph)412 MaxMatching(const Graph& graph) 413 : _graph(graph), _matching(0), _status(0), _ear(0), 414 _blossom_set_index(0), _blossom_set(0), _blossom_rep(0), 415 _tree_set_index(0), _tree_set(0) {} 416 ~MaxMatching()417 ~MaxMatching() { 418 destroyStructures(); 419 } 420 421 /// \name Execution Control 422 /// The simplest way to execute the algorithm is to use the 423 /// \c run() member function.\n 424 /// If you need better control on the execution, you have to call 425 /// one of the functions \ref init(), \ref greedyInit() or 426 /// \ref matchingInit() first, then you can start the algorithm with 427 /// \ref startSparse() or \ref startDense(). 428 429 ///@{ 430 431 /// \brief Set the initial matching to the empty matching. 432 /// 433 /// This function sets the initial matching to the empty matching. init()434 void init() { 435 createStructures(); 436 for(NodeIt n(_graph); n != INVALID; ++n) { 437 (*_matching)[n] = INVALID; 438 (*_status)[n] = UNMATCHED; 439 } 440 } 441 442 /// \brief Find an initial matching in a greedy way. 443 /// 444 /// This function finds an initial matching in a greedy way. greedyInit()445 void greedyInit() { 446 createStructures(); 447 for (NodeIt n(_graph); n != INVALID; ++n) { 448 (*_matching)[n] = INVALID; 449 (*_status)[n] = UNMATCHED; 450 } 451 for (NodeIt n(_graph); n != INVALID; ++n) { 452 if ((*_matching)[n] == INVALID) { 453 for (OutArcIt a(_graph, n); a != INVALID ; ++a) { 454 Node v = _graph.target(a); 455 if ((*_matching)[v] == INVALID && v != n) { 456 (*_matching)[n] = a; 457 (*_status)[n] = MATCHED; 458 (*_matching)[v] = _graph.oppositeArc(a); 459 (*_status)[v] = MATCHED; 460 break; 461 } 462 } 463 } 464 } 465 } 466 467 468 /// \brief Initialize the matching from a map. 469 /// 470 /// This function initializes the matching from a \c bool valued edge 471 /// map. This map should have the property that there are no two incident 472 /// edges with \c true value, i.e. it really contains a matching. 473 /// \return \c true if the map contains a matching. 474 template <typename MatchingMap> matchingInit(const MatchingMap & matching)475 bool matchingInit(const MatchingMap& matching) { 476 createStructures(); 477 478 for (NodeIt n(_graph); n != INVALID; ++n) { 479 (*_matching)[n] = INVALID; 480 (*_status)[n] = UNMATCHED; 481 } 482 for(EdgeIt e(_graph); e!=INVALID; ++e) { 483 if (matching[e]) { 484 485 Node u = _graph.u(e); 486 if ((*_matching)[u] != INVALID) return false; 487 (*_matching)[u] = _graph.direct(e, true); 488 (*_status)[u] = MATCHED; 489 490 Node v = _graph.v(e); 491 if ((*_matching)[v] != INVALID) return false; 492 (*_matching)[v] = _graph.direct(e, false); 493 (*_status)[v] = MATCHED; 494 } 495 } 496 return true; 497 } 498 499 /// \brief Start Edmonds' algorithm 500 /// 501 /// This function runs the original Edmonds' algorithm. 502 /// 503 /// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be 504 /// called before using this function. startSparse()505 void startSparse() { 506 for(NodeIt n(_graph); n != INVALID; ++n) { 507 if ((*_status)[n] == UNMATCHED) { 508 (*_blossom_rep)[_blossom_set->insert(n)] = n; 509 _tree_set->insert(n); 510 (*_status)[n] = EVEN; 511 processSparse(n); 512 } 513 } 514 } 515 516 /// \brief Start Edmonds' algorithm with a heuristic improvement 517 /// for dense graphs 518 /// 519 /// This function runs Edmonds' algorithm with a heuristic of postponing 520 /// shrinks, therefore resulting in a faster algorithm for dense graphs. 521 /// 522 /// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be 523 /// called before using this function. startDense()524 void startDense() { 525 for(NodeIt n(_graph); n != INVALID; ++n) { 526 if ((*_status)[n] == UNMATCHED) { 527 (*_blossom_rep)[_blossom_set->insert(n)] = n; 528 _tree_set->insert(n); 529 (*_status)[n] = EVEN; 530 processDense(n); 531 } 532 } 533 } 534 535 536 /// \brief Run Edmonds' algorithm 537 /// 538 /// This function runs Edmonds' algorithm. An additional heuristic of 539 /// postponing shrinks is used for relatively dense graphs 540 /// (for which <tt>m>=2*n</tt> holds). run()541 void run() { 542 if (countEdges(_graph) < 2 * countNodes(_graph)) { 543 greedyInit(); 544 startSparse(); 545 } else { 546 init(); 547 startDense(); 548 } 549 } 550 551 /// @} 552 553 /// \name Primal Solution 554 /// Functions to get the primal solution, i.e. the maximum matching. 555 556 /// @{ 557 558 /// \brief Return the size (cardinality) of the matching. 559 /// 560 /// This function returns the size (cardinality) of the current matching. 561 /// After run() it returns the size of the maximum matching in the graph. matchingSize()562 int matchingSize() const { 563 int size = 0; 564 for (NodeIt n(_graph); n != INVALID; ++n) { 565 if ((*_matching)[n] != INVALID) { 566 ++size; 567 } 568 } 569 return size / 2; 570 } 571 572 /// \brief Return \c true if the given edge is in the matching. 573 /// 574 /// This function returns \c true if the given edge is in the current 575 /// matching. matching(const Edge & edge)576 bool matching(const Edge& edge) const { 577 return edge == (*_matching)[_graph.u(edge)]; 578 } 579 580 /// \brief Return the matching arc (or edge) incident to the given node. 581 /// 582 /// This function returns the matching arc (or edge) incident to the 583 /// given node in the current matching or \c INVALID if the node is 584 /// not covered by the matching. matching(const Node & n)585 Arc matching(const Node& n) const { 586 return (*_matching)[n]; 587 } 588 589 /// \brief Return a const reference to the matching map. 590 /// 591 /// This function returns a const reference to a node map that stores 592 /// the matching arc (or edge) incident to each node. matchingMap()593 const MatchingMap& matchingMap() const { 594 return *_matching; 595 } 596 597 /// \brief Return the mate of the given node. 598 /// 599 /// This function returns the mate of the given node in the current 600 /// matching or \c INVALID if the node is not covered by the matching. mate(const Node & n)601 Node mate(const Node& n) const { 602 return (*_matching)[n] != INVALID ? 603 _graph.target((*_matching)[n]) : INVALID; 604 } 605 606 /// @} 607 608 /// \name Dual Solution 609 /// Functions to get the dual solution, i.e. the Gallai-Edmonds 610 /// decomposition. 611 612 /// @{ 613 614 /// \brief Return the status of the given node in the Edmonds-Gallai 615 /// decomposition. 616 /// 617 /// This function returns the \ref Status "status" of the given node 618 /// in the Edmonds-Gallai decomposition. status(const Node & n)619 Status status(const Node& n) const { 620 return (*_status)[n]; 621 } 622 623 /// \brief Return a const reference to the status map, which stores 624 /// the Edmonds-Gallai decomposition. 625 /// 626 /// This function returns a const reference to a node map that stores the 627 /// \ref Status "status" of each node in the Edmonds-Gallai decomposition. statusMap()628 const StatusMap& statusMap() const { 629 return *_status; 630 } 631 632 /// \brief Return \c true if the given node is in the barrier. 633 /// 634 /// This function returns \c true if the given node is in the barrier. barrier(const Node & n)635 bool barrier(const Node& n) const { 636 return (*_status)[n] == ODD; 637 } 638 639 /// @} 640 641 }; 642 643 /// \ingroup matching 644 /// 645 /// \brief Weighted matching in general graphs 646 /// 647 /// This class provides an efficient implementation of Edmond's 648 /// maximum weighted matching algorithm. The implementation is based 649 /// on extensive use of priority queues and provides 650 /// \f$O(nm\log n)\f$ time complexity. 651 /// 652 /// The maximum weighted matching problem is to find a subset of the 653 /// edges in an undirected graph with maximum overall weight for which 654 /// each node has at most one incident edge. 655 /// It can be formulated with the following linear program. 656 /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f] 657 /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} 658 \quad \forall B\in\mathcal{O}\f] */ 659 /// \f[x_e \ge 0\quad \forall e\in E\f] 660 /// \f[\max \sum_{e\in E}x_ew_e\f] 661 /// where \f$\delta(X)\f$ is the set of edges incident to a node in 662 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in 663 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality 664 /// subsets of the nodes. 665 /// 666 /// The algorithm calculates an optimal matching and a proof of the 667 /// optimality. The solution of the dual problem can be used to check 668 /// the result of the algorithm. The dual linear problem is the 669 /// following. 670 /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)} 671 z_B \ge w_{uv} \quad \forall uv\in E\f] */ 672 /// \f[y_u \ge 0 \quad \forall u \in V\f] 673 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] 674 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} 675 \frac{\vert B \vert - 1}{2}z_B\f] */ 676 /// 677 /// The algorithm can be executed with the run() function. 678 /// After it the matching (the primal solution) and the dual solution 679 /// can be obtained using the query functions and the 680 /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, 681 /// which is able to iterate on the nodes of a blossom. 682 /// If the value type is integer, then the dual solution is multiplied 683 /// by \ref MaxWeightedMatching::dualScale "4". 684 /// 685 /// \tparam GR The undirected graph type the algorithm runs on. 686 /// \tparam WM The type edge weight map. The default type is 687 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>". 688 #ifdef DOXYGEN 689 template <typename GR, typename WM> 690 #else 691 template <typename GR, 692 typename WM = typename GR::template EdgeMap<int> > 693 #endif 694 class MaxWeightedMatching { 695 public: 696 697 /// The graph type of the algorithm 698 typedef GR Graph; 699 /// The type of the edge weight map 700 typedef WM WeightMap; 701 /// The value type of the edge weights 702 typedef typename WeightMap::Value Value; 703 704 /// The type of the matching map 705 typedef typename Graph::template NodeMap<typename Graph::Arc> 706 MatchingMap; 707 708 /// \brief Scaling factor for dual solution 709 /// 710 /// Scaling factor for dual solution. It is equal to 4 or 1 711 /// according to the value type. 712 static const int dualScale = 713 std::numeric_limits<Value>::is_integer ? 4 : 1; 714 715 private: 716 717 TEMPLATE_GRAPH_TYPEDEFS(Graph); 718 719 typedef typename Graph::template NodeMap<Value> NodePotential; 720 typedef std::vector<Node> BlossomNodeList; 721 722 struct BlossomVariable { 723 int begin, end; 724 Value value; 725 BlossomVariableBlossomVariable726 BlossomVariable(int _begin, int _end, Value _value) 727 : begin(_begin), end(_end), value(_value) {} 728 729 }; 730 731 typedef std::vector<BlossomVariable> BlossomPotential; 732 733 const Graph& _graph; 734 const WeightMap& _weight; 735 736 MatchingMap* _matching; 737 738 NodePotential* _node_potential; 739 740 BlossomPotential _blossom_potential; 741 BlossomNodeList _blossom_node_list; 742 743 int _node_num; 744 int _blossom_num; 745 746 typedef RangeMap<int> IntIntMap; 747 748 enum Status { 749 EVEN = -1, MATCHED = 0, ODD = 1 750 }; 751 752 typedef HeapUnionFind<Value, IntNodeMap> BlossomSet; 753 struct BlossomData { 754 int tree; 755 Status status; 756 Arc pred, next; 757 Value pot, offset; 758 Node base; 759 }; 760 761 IntNodeMap *_blossom_index; 762 BlossomSet *_blossom_set; 763 RangeMap<BlossomData>* _blossom_data; 764 765 IntNodeMap *_node_index; 766 IntArcMap *_node_heap_index; 767 768 struct NodeData { 769 NodeDataNodeData770 NodeData(IntArcMap& node_heap_index) 771 : heap(node_heap_index) {} 772 773 int blossom; 774 Value pot; 775 BinHeap<Value, IntArcMap> heap; 776 std::map<int, Arc> heap_index; 777 778 int tree; 779 }; 780 781 RangeMap<NodeData>* _node_data; 782 783 typedef ExtendFindEnum<IntIntMap> TreeSet; 784 785 IntIntMap *_tree_set_index; 786 TreeSet *_tree_set; 787 788 IntNodeMap *_delta1_index; 789 BinHeap<Value, IntNodeMap> *_delta1; 790 791 IntIntMap *_delta2_index; 792 BinHeap<Value, IntIntMap> *_delta2; 793 794 IntEdgeMap *_delta3_index; 795 BinHeap<Value, IntEdgeMap> *_delta3; 796 797 IntIntMap *_delta4_index; 798 BinHeap<Value, IntIntMap> *_delta4; 799 800 Value _delta_sum; 801 int _unmatched; 802 803 typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching; 804 FractionalMatching *_fractional; 805 createStructures()806 void createStructures() { 807 _node_num = countNodes(_graph); 808 _blossom_num = _node_num * 3 / 2; 809 810 if (!_matching) { 811 _matching = new MatchingMap(_graph); 812 } 813 814 if (!_node_potential) { 815 _node_potential = new NodePotential(_graph); 816 } 817 818 if (!_blossom_set) { 819 _blossom_index = new IntNodeMap(_graph); 820 _blossom_set = new BlossomSet(*_blossom_index); 821 _blossom_data = new RangeMap<BlossomData>(_blossom_num); 822 } else if (_blossom_data->size() != _blossom_num) { 823 delete _blossom_data; 824 _blossom_data = new RangeMap<BlossomData>(_blossom_num); 825 } 826 827 if (!_node_index) { 828 _node_index = new IntNodeMap(_graph); 829 _node_heap_index = new IntArcMap(_graph); 830 _node_data = new RangeMap<NodeData>(_node_num, 831 NodeData(*_node_heap_index)); 832 } else { 833 delete _node_data; 834 _node_data = new RangeMap<NodeData>(_node_num, 835 NodeData(*_node_heap_index)); 836 } 837 838 if (!_tree_set) { 839 _tree_set_index = new IntIntMap(_blossom_num); 840 _tree_set = new TreeSet(*_tree_set_index); 841 } else { 842 _tree_set_index->resize(_blossom_num); 843 } 844 845 if (!_delta1) { 846 _delta1_index = new IntNodeMap(_graph); 847 _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index); 848 } 849 850 if (!_delta2) { 851 _delta2_index = new IntIntMap(_blossom_num); 852 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index); 853 } else { 854 _delta2_index->resize(_blossom_num); 855 } 856 857 if (!_delta3) { 858 _delta3_index = new IntEdgeMap(_graph); 859 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index); 860 } 861 862 if (!_delta4) { 863 _delta4_index = new IntIntMap(_blossom_num); 864 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index); 865 } else { 866 _delta4_index->resize(_blossom_num); 867 } 868 } 869 destroyStructures()870 void destroyStructures() { 871 if (_matching) { 872 delete _matching; 873 } 874 if (_node_potential) { 875 delete _node_potential; 876 } 877 if (_blossom_set) { 878 delete _blossom_index; 879 delete _blossom_set; 880 delete _blossom_data; 881 } 882 883 if (_node_index) { 884 delete _node_index; 885 delete _node_heap_index; 886 delete _node_data; 887 } 888 889 if (_tree_set) { 890 delete _tree_set_index; 891 delete _tree_set; 892 } 893 if (_delta1) { 894 delete _delta1_index; 895 delete _delta1; 896 } 897 if (_delta2) { 898 delete _delta2_index; 899 delete _delta2; 900 } 901 if (_delta3) { 902 delete _delta3_index; 903 delete _delta3; 904 } 905 if (_delta4) { 906 delete _delta4_index; 907 delete _delta4; 908 } 909 } 910 matchedToEven(int blossom,int tree)911 void matchedToEven(int blossom, int tree) { 912 if (_delta2->state(blossom) == _delta2->IN_HEAP) { 913 _delta2->erase(blossom); 914 } 915 916 if (!_blossom_set->trivial(blossom)) { 917 (*_blossom_data)[blossom].pot -= 918 2 * (_delta_sum - (*_blossom_data)[blossom].offset); 919 } 920 921 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 922 n != INVALID; ++n) { 923 924 _blossom_set->increase(n, std::numeric_limits<Value>::max()); 925 int ni = (*_node_index)[n]; 926 927 (*_node_data)[ni].heap.clear(); 928 (*_node_data)[ni].heap_index.clear(); 929 930 (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset; 931 932 _delta1->push(n, (*_node_data)[ni].pot); 933 934 for (InArcIt e(_graph, n); e != INVALID; ++e) { 935 Node v = _graph.source(e); 936 int vb = _blossom_set->find(v); 937 int vi = (*_node_index)[v]; 938 939 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 940 dualScale * _weight[e]; 941 942 if ((*_blossom_data)[vb].status == EVEN) { 943 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { 944 _delta3->push(e, rw / 2); 945 } 946 } else { 947 typename std::map<int, Arc>::iterator it = 948 (*_node_data)[vi].heap_index.find(tree); 949 950 if (it != (*_node_data)[vi].heap_index.end()) { 951 if ((*_node_data)[vi].heap[it->second] > rw) { 952 (*_node_data)[vi].heap.replace(it->second, e); 953 (*_node_data)[vi].heap.decrease(e, rw); 954 it->second = e; 955 } 956 } else { 957 (*_node_data)[vi].heap.push(e, rw); 958 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); 959 } 960 961 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { 962 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); 963 964 if ((*_blossom_data)[vb].status == MATCHED) { 965 if (_delta2->state(vb) != _delta2->IN_HEAP) { 966 _delta2->push(vb, _blossom_set->classPrio(vb) - 967 (*_blossom_data)[vb].offset); 968 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - 969 (*_blossom_data)[vb].offset) { 970 _delta2->decrease(vb, _blossom_set->classPrio(vb) - 971 (*_blossom_data)[vb].offset); 972 } 973 } 974 } 975 } 976 } 977 } 978 (*_blossom_data)[blossom].offset = 0; 979 } 980 matchedToOdd(int blossom)981 void matchedToOdd(int blossom) { 982 if (_delta2->state(blossom) == _delta2->IN_HEAP) { 983 _delta2->erase(blossom); 984 } 985 (*_blossom_data)[blossom].offset += _delta_sum; 986 if (!_blossom_set->trivial(blossom)) { 987 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + 988 (*_blossom_data)[blossom].offset); 989 } 990 } 991 evenToMatched(int blossom,int tree)992 void evenToMatched(int blossom, int tree) { 993 if (!_blossom_set->trivial(blossom)) { 994 (*_blossom_data)[blossom].pot += 2 * _delta_sum; 995 } 996 997 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 998 n != INVALID; ++n) { 999 int ni = (*_node_index)[n]; 1000 (*_node_data)[ni].pot -= _delta_sum; 1001 1002 _delta1->erase(n); 1003 1004 for (InArcIt e(_graph, n); e != INVALID; ++e) { 1005 Node v = _graph.source(e); 1006 int vb = _blossom_set->find(v); 1007 int vi = (*_node_index)[v]; 1008 1009 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 1010 dualScale * _weight[e]; 1011 1012 if (vb == blossom) { 1013 if (_delta3->state(e) == _delta3->IN_HEAP) { 1014 _delta3->erase(e); 1015 } 1016 } else if ((*_blossom_data)[vb].status == EVEN) { 1017 1018 if (_delta3->state(e) == _delta3->IN_HEAP) { 1019 _delta3->erase(e); 1020 } 1021 1022 int vt = _tree_set->find(vb); 1023 1024 if (vt != tree) { 1025 1026 Arc r = _graph.oppositeArc(e); 1027 1028 typename std::map<int, Arc>::iterator it = 1029 (*_node_data)[ni].heap_index.find(vt); 1030 1031 if (it != (*_node_data)[ni].heap_index.end()) { 1032 if ((*_node_data)[ni].heap[it->second] > rw) { 1033 (*_node_data)[ni].heap.replace(it->second, r); 1034 (*_node_data)[ni].heap.decrease(r, rw); 1035 it->second = r; 1036 } 1037 } else { 1038 (*_node_data)[ni].heap.push(r, rw); 1039 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); 1040 } 1041 1042 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { 1043 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); 1044 1045 if (_delta2->state(blossom) != _delta2->IN_HEAP) { 1046 _delta2->push(blossom, _blossom_set->classPrio(blossom) - 1047 (*_blossom_data)[blossom].offset); 1048 } else if ((*_delta2)[blossom] > 1049 _blossom_set->classPrio(blossom) - 1050 (*_blossom_data)[blossom].offset){ 1051 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - 1052 (*_blossom_data)[blossom].offset); 1053 } 1054 } 1055 } 1056 } else { 1057 1058 typename std::map<int, Arc>::iterator it = 1059 (*_node_data)[vi].heap_index.find(tree); 1060 1061 if (it != (*_node_data)[vi].heap_index.end()) { 1062 (*_node_data)[vi].heap.erase(it->second); 1063 (*_node_data)[vi].heap_index.erase(it); 1064 if ((*_node_data)[vi].heap.empty()) { 1065 _blossom_set->increase(v, std::numeric_limits<Value>::max()); 1066 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { 1067 _blossom_set->increase(v, (*_node_data)[vi].heap.prio()); 1068 } 1069 1070 if ((*_blossom_data)[vb].status == MATCHED) { 1071 if (_blossom_set->classPrio(vb) == 1072 std::numeric_limits<Value>::max()) { 1073 _delta2->erase(vb); 1074 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - 1075 (*_blossom_data)[vb].offset) { 1076 _delta2->increase(vb, _blossom_set->classPrio(vb) - 1077 (*_blossom_data)[vb].offset); 1078 } 1079 } 1080 } 1081 } 1082 } 1083 } 1084 } 1085 oddToMatched(int blossom)1086 void oddToMatched(int blossom) { 1087 (*_blossom_data)[blossom].offset -= _delta_sum; 1088 1089 if (_blossom_set->classPrio(blossom) != 1090 std::numeric_limits<Value>::max()) { 1091 _delta2->push(blossom, _blossom_set->classPrio(blossom) - 1092 (*_blossom_data)[blossom].offset); 1093 } 1094 1095 if (!_blossom_set->trivial(blossom)) { 1096 _delta4->erase(blossom); 1097 } 1098 } 1099 oddToEven(int blossom,int tree)1100 void oddToEven(int blossom, int tree) { 1101 if (!_blossom_set->trivial(blossom)) { 1102 _delta4->erase(blossom); 1103 (*_blossom_data)[blossom].pot -= 1104 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); 1105 } 1106 1107 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 1108 n != INVALID; ++n) { 1109 int ni = (*_node_index)[n]; 1110 1111 _blossom_set->increase(n, std::numeric_limits<Value>::max()); 1112 1113 (*_node_data)[ni].heap.clear(); 1114 (*_node_data)[ni].heap_index.clear(); 1115 (*_node_data)[ni].pot += 1116 2 * _delta_sum - (*_blossom_data)[blossom].offset; 1117 1118 _delta1->push(n, (*_node_data)[ni].pot); 1119 1120 for (InArcIt e(_graph, n); e != INVALID; ++e) { 1121 Node v = _graph.source(e); 1122 int vb = _blossom_set->find(v); 1123 int vi = (*_node_index)[v]; 1124 1125 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 1126 dualScale * _weight[e]; 1127 1128 if ((*_blossom_data)[vb].status == EVEN) { 1129 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { 1130 _delta3->push(e, rw / 2); 1131 } 1132 } else { 1133 1134 typename std::map<int, Arc>::iterator it = 1135 (*_node_data)[vi].heap_index.find(tree); 1136 1137 if (it != (*_node_data)[vi].heap_index.end()) { 1138 if ((*_node_data)[vi].heap[it->second] > rw) { 1139 (*_node_data)[vi].heap.replace(it->second, e); 1140 (*_node_data)[vi].heap.decrease(e, rw); 1141 it->second = e; 1142 } 1143 } else { 1144 (*_node_data)[vi].heap.push(e, rw); 1145 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); 1146 } 1147 1148 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { 1149 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); 1150 1151 if ((*_blossom_data)[vb].status == MATCHED) { 1152 if (_delta2->state(vb) != _delta2->IN_HEAP) { 1153 _delta2->push(vb, _blossom_set->classPrio(vb) - 1154 (*_blossom_data)[vb].offset); 1155 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - 1156 (*_blossom_data)[vb].offset) { 1157 _delta2->decrease(vb, _blossom_set->classPrio(vb) - 1158 (*_blossom_data)[vb].offset); 1159 } 1160 } 1161 } 1162 } 1163 } 1164 } 1165 (*_blossom_data)[blossom].offset = 0; 1166 } 1167 alternatePath(int even,int tree)1168 void alternatePath(int even, int tree) { 1169 int odd; 1170 1171 evenToMatched(even, tree); 1172 (*_blossom_data)[even].status = MATCHED; 1173 1174 while ((*_blossom_data)[even].pred != INVALID) { 1175 odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred)); 1176 (*_blossom_data)[odd].status = MATCHED; 1177 oddToMatched(odd); 1178 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred; 1179 1180 even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred)); 1181 (*_blossom_data)[even].status = MATCHED; 1182 evenToMatched(even, tree); 1183 (*_blossom_data)[even].next = 1184 _graph.oppositeArc((*_blossom_data)[odd].pred); 1185 } 1186 1187 } 1188 destroyTree(int tree)1189 void destroyTree(int tree) { 1190 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) { 1191 if ((*_blossom_data)[b].status == EVEN) { 1192 (*_blossom_data)[b].status = MATCHED; 1193 evenToMatched(b, tree); 1194 } else if ((*_blossom_data)[b].status == ODD) { 1195 (*_blossom_data)[b].status = MATCHED; 1196 oddToMatched(b); 1197 } 1198 } 1199 _tree_set->eraseClass(tree); 1200 } 1201 1202 unmatchNode(const Node & node)1203 void unmatchNode(const Node& node) { 1204 int blossom = _blossom_set->find(node); 1205 int tree = _tree_set->find(blossom); 1206 1207 alternatePath(blossom, tree); 1208 destroyTree(tree); 1209 1210 (*_blossom_data)[blossom].base = node; 1211 (*_blossom_data)[blossom].next = INVALID; 1212 } 1213 augmentOnEdge(const Edge & edge)1214 void augmentOnEdge(const Edge& edge) { 1215 1216 int left = _blossom_set->find(_graph.u(edge)); 1217 int right = _blossom_set->find(_graph.v(edge)); 1218 1219 int left_tree = _tree_set->find(left); 1220 alternatePath(left, left_tree); 1221 destroyTree(left_tree); 1222 1223 int right_tree = _tree_set->find(right); 1224 alternatePath(right, right_tree); 1225 destroyTree(right_tree); 1226 1227 (*_blossom_data)[left].next = _graph.direct(edge, true); 1228 (*_blossom_data)[right].next = _graph.direct(edge, false); 1229 } 1230 augmentOnArc(const Arc & arc)1231 void augmentOnArc(const Arc& arc) { 1232 1233 int left = _blossom_set->find(_graph.source(arc)); 1234 int right = _blossom_set->find(_graph.target(arc)); 1235 1236 (*_blossom_data)[left].status = MATCHED; 1237 1238 int right_tree = _tree_set->find(right); 1239 alternatePath(right, right_tree); 1240 destroyTree(right_tree); 1241 1242 (*_blossom_data)[left].next = arc; 1243 (*_blossom_data)[right].next = _graph.oppositeArc(arc); 1244 } 1245 extendOnArc(const Arc & arc)1246 void extendOnArc(const Arc& arc) { 1247 int base = _blossom_set->find(_graph.target(arc)); 1248 int tree = _tree_set->find(base); 1249 1250 int odd = _blossom_set->find(_graph.source(arc)); 1251 _tree_set->insert(odd, tree); 1252 (*_blossom_data)[odd].status = ODD; 1253 matchedToOdd(odd); 1254 (*_blossom_data)[odd].pred = arc; 1255 1256 int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next)); 1257 (*_blossom_data)[even].pred = (*_blossom_data)[even].next; 1258 _tree_set->insert(even, tree); 1259 (*_blossom_data)[even].status = EVEN; 1260 matchedToEven(even, tree); 1261 } 1262 shrinkOnEdge(const Edge & edge,int tree)1263 void shrinkOnEdge(const Edge& edge, int tree) { 1264 int nca = -1; 1265 std::vector<int> left_path, right_path; 1266 1267 { 1268 std::set<int> left_set, right_set; 1269 int left = _blossom_set->find(_graph.u(edge)); 1270 left_path.push_back(left); 1271 left_set.insert(left); 1272 1273 int right = _blossom_set->find(_graph.v(edge)); 1274 right_path.push_back(right); 1275 right_set.insert(right); 1276 1277 while (true) { 1278 1279 if ((*_blossom_data)[left].pred == INVALID) break; 1280 1281 left = 1282 _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); 1283 left_path.push_back(left); 1284 left = 1285 _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); 1286 left_path.push_back(left); 1287 1288 left_set.insert(left); 1289 1290 if (right_set.find(left) != right_set.end()) { 1291 nca = left; 1292 break; 1293 } 1294 1295 if ((*_blossom_data)[right].pred == INVALID) break; 1296 1297 right = 1298 _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); 1299 right_path.push_back(right); 1300 right = 1301 _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); 1302 right_path.push_back(right); 1303 1304 right_set.insert(right); 1305 1306 if (left_set.find(right) != left_set.end()) { 1307 nca = right; 1308 break; 1309 } 1310 1311 } 1312 1313 if (nca == -1) { 1314 if ((*_blossom_data)[left].pred == INVALID) { 1315 nca = right; 1316 while (left_set.find(nca) == left_set.end()) { 1317 nca = 1318 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); 1319 right_path.push_back(nca); 1320 nca = 1321 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); 1322 right_path.push_back(nca); 1323 } 1324 } else { 1325 nca = left; 1326 while (right_set.find(nca) == right_set.end()) { 1327 nca = 1328 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); 1329 left_path.push_back(nca); 1330 nca = 1331 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); 1332 left_path.push_back(nca); 1333 } 1334 } 1335 } 1336 } 1337 1338 std::vector<int> subblossoms; 1339 Arc prev; 1340 1341 prev = _graph.direct(edge, true); 1342 for (int i = 0; left_path[i] != nca; i += 2) { 1343 subblossoms.push_back(left_path[i]); 1344 (*_blossom_data)[left_path[i]].next = prev; 1345 _tree_set->erase(left_path[i]); 1346 1347 subblossoms.push_back(left_path[i + 1]); 1348 (*_blossom_data)[left_path[i + 1]].status = EVEN; 1349 oddToEven(left_path[i + 1], tree); 1350 _tree_set->erase(left_path[i + 1]); 1351 prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred); 1352 } 1353 1354 int k = 0; 1355 while (right_path[k] != nca) ++k; 1356 1357 subblossoms.push_back(nca); 1358 (*_blossom_data)[nca].next = prev; 1359 1360 for (int i = k - 2; i >= 0; i -= 2) { 1361 subblossoms.push_back(right_path[i + 1]); 1362 (*_blossom_data)[right_path[i + 1]].status = EVEN; 1363 oddToEven(right_path[i + 1], tree); 1364 _tree_set->erase(right_path[i + 1]); 1365 1366 (*_blossom_data)[right_path[i + 1]].next = 1367 (*_blossom_data)[right_path[i + 1]].pred; 1368 1369 subblossoms.push_back(right_path[i]); 1370 _tree_set->erase(right_path[i]); 1371 } 1372 1373 int surface = 1374 _blossom_set->join(subblossoms.begin(), subblossoms.end()); 1375 1376 for (int i = 0; i < int(subblossoms.size()); ++i) { 1377 if (!_blossom_set->trivial(subblossoms[i])) { 1378 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum; 1379 } 1380 (*_blossom_data)[subblossoms[i]].status = MATCHED; 1381 } 1382 1383 (*_blossom_data)[surface].pot = -2 * _delta_sum; 1384 (*_blossom_data)[surface].offset = 0; 1385 (*_blossom_data)[surface].status = EVEN; 1386 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred; 1387 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred; 1388 1389 _tree_set->insert(surface, tree); 1390 _tree_set->erase(nca); 1391 } 1392 splitBlossom(int blossom)1393 void splitBlossom(int blossom) { 1394 Arc next = (*_blossom_data)[blossom].next; 1395 Arc pred = (*_blossom_data)[blossom].pred; 1396 1397 int tree = _tree_set->find(blossom); 1398 1399 (*_blossom_data)[blossom].status = MATCHED; 1400 oddToMatched(blossom); 1401 if (_delta2->state(blossom) == _delta2->IN_HEAP) { 1402 _delta2->erase(blossom); 1403 } 1404 1405 std::vector<int> subblossoms; 1406 _blossom_set->split(blossom, std::back_inserter(subblossoms)); 1407 1408 Value offset = (*_blossom_data)[blossom].offset; 1409 int b = _blossom_set->find(_graph.source(pred)); 1410 int d = _blossom_set->find(_graph.source(next)); 1411 1412 int ib = -1, id = -1; 1413 for (int i = 0; i < int(subblossoms.size()); ++i) { 1414 if (subblossoms[i] == b) ib = i; 1415 if (subblossoms[i] == d) id = i; 1416 1417 (*_blossom_data)[subblossoms[i]].offset = offset; 1418 if (!_blossom_set->trivial(subblossoms[i])) { 1419 (*_blossom_data)[subblossoms[i]].pot -= 2 * offset; 1420 } 1421 if (_blossom_set->classPrio(subblossoms[i]) != 1422 std::numeric_limits<Value>::max()) { 1423 _delta2->push(subblossoms[i], 1424 _blossom_set->classPrio(subblossoms[i]) - 1425 (*_blossom_data)[subblossoms[i]].offset); 1426 } 1427 } 1428 1429 if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) { 1430 for (int i = (id + 1) % subblossoms.size(); 1431 i != ib; i = (i + 2) % subblossoms.size()) { 1432 int sb = subblossoms[i]; 1433 int tb = subblossoms[(i + 1) % subblossoms.size()]; 1434 (*_blossom_data)[sb].next = 1435 _graph.oppositeArc((*_blossom_data)[tb].next); 1436 } 1437 1438 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) { 1439 int sb = subblossoms[i]; 1440 int tb = subblossoms[(i + 1) % subblossoms.size()]; 1441 int ub = subblossoms[(i + 2) % subblossoms.size()]; 1442 1443 (*_blossom_data)[sb].status = ODD; 1444 matchedToOdd(sb); 1445 _tree_set->insert(sb, tree); 1446 (*_blossom_data)[sb].pred = pred; 1447 (*_blossom_data)[sb].next = 1448 _graph.oppositeArc((*_blossom_data)[tb].next); 1449 1450 pred = (*_blossom_data)[ub].next; 1451 1452 (*_blossom_data)[tb].status = EVEN; 1453 matchedToEven(tb, tree); 1454 _tree_set->insert(tb, tree); 1455 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next; 1456 } 1457 1458 (*_blossom_data)[subblossoms[id]].status = ODD; 1459 matchedToOdd(subblossoms[id]); 1460 _tree_set->insert(subblossoms[id], tree); 1461 (*_blossom_data)[subblossoms[id]].next = next; 1462 (*_blossom_data)[subblossoms[id]].pred = pred; 1463 1464 } else { 1465 1466 for (int i = (ib + 1) % subblossoms.size(); 1467 i != id; i = (i + 2) % subblossoms.size()) { 1468 int sb = subblossoms[i]; 1469 int tb = subblossoms[(i + 1) % subblossoms.size()]; 1470 (*_blossom_data)[sb].next = 1471 _graph.oppositeArc((*_blossom_data)[tb].next); 1472 } 1473 1474 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) { 1475 int sb = subblossoms[i]; 1476 int tb = subblossoms[(i + 1) % subblossoms.size()]; 1477 int ub = subblossoms[(i + 2) % subblossoms.size()]; 1478 1479 (*_blossom_data)[sb].status = ODD; 1480 matchedToOdd(sb); 1481 _tree_set->insert(sb, tree); 1482 (*_blossom_data)[sb].next = next; 1483 (*_blossom_data)[sb].pred = 1484 _graph.oppositeArc((*_blossom_data)[tb].next); 1485 1486 (*_blossom_data)[tb].status = EVEN; 1487 matchedToEven(tb, tree); 1488 _tree_set->insert(tb, tree); 1489 (*_blossom_data)[tb].pred = 1490 (*_blossom_data)[tb].next = 1491 _graph.oppositeArc((*_blossom_data)[ub].next); 1492 next = (*_blossom_data)[ub].next; 1493 } 1494 1495 (*_blossom_data)[subblossoms[ib]].status = ODD; 1496 matchedToOdd(subblossoms[ib]); 1497 _tree_set->insert(subblossoms[ib], tree); 1498 (*_blossom_data)[subblossoms[ib]].next = next; 1499 (*_blossom_data)[subblossoms[ib]].pred = pred; 1500 } 1501 _tree_set->erase(blossom); 1502 } 1503 extractBlossom(int blossom,const Node & base,const Arc & matching)1504 void extractBlossom(int blossom, const Node& base, const Arc& matching) { 1505 if (_blossom_set->trivial(blossom)) { 1506 int bi = (*_node_index)[base]; 1507 Value pot = (*_node_data)[bi].pot; 1508 1509 (*_matching)[base] = matching; 1510 _blossom_node_list.push_back(base); 1511 (*_node_potential)[base] = pot; 1512 } else { 1513 1514 Value pot = (*_blossom_data)[blossom].pot; 1515 int bn = _blossom_node_list.size(); 1516 1517 std::vector<int> subblossoms; 1518 _blossom_set->split(blossom, std::back_inserter(subblossoms)); 1519 int b = _blossom_set->find(base); 1520 int ib = -1; 1521 for (int i = 0; i < int(subblossoms.size()); ++i) { 1522 if (subblossoms[i] == b) { ib = i; break; } 1523 } 1524 1525 for (int i = 1; i < int(subblossoms.size()); i += 2) { 1526 int sb = subblossoms[(ib + i) % subblossoms.size()]; 1527 int tb = subblossoms[(ib + i + 1) % subblossoms.size()]; 1528 1529 Arc m = (*_blossom_data)[tb].next; 1530 extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m)); 1531 extractBlossom(tb, _graph.source(m), m); 1532 } 1533 extractBlossom(subblossoms[ib], base, matching); 1534 1535 int en = _blossom_node_list.size(); 1536 1537 _blossom_potential.push_back(BlossomVariable(bn, en, pot)); 1538 } 1539 } 1540 extractMatching()1541 void extractMatching() { 1542 std::vector<int> blossoms; 1543 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) { 1544 blossoms.push_back(c); 1545 } 1546 1547 for (int i = 0; i < int(blossoms.size()); ++i) { 1548 if ((*_blossom_data)[blossoms[i]].next != INVALID) { 1549 1550 Value offset = (*_blossom_data)[blossoms[i]].offset; 1551 (*_blossom_data)[blossoms[i]].pot += 2 * offset; 1552 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); 1553 n != INVALID; ++n) { 1554 (*_node_data)[(*_node_index)[n]].pot -= offset; 1555 } 1556 1557 Arc matching = (*_blossom_data)[blossoms[i]].next; 1558 Node base = _graph.source(matching); 1559 extractBlossom(blossoms[i], base, matching); 1560 } else { 1561 Node base = (*_blossom_data)[blossoms[i]].base; 1562 extractBlossom(blossoms[i], base, INVALID); 1563 } 1564 } 1565 } 1566 1567 public: 1568 1569 /// \brief Constructor 1570 /// 1571 /// Constructor. MaxWeightedMatching(const Graph & graph,const WeightMap & weight)1572 MaxWeightedMatching(const Graph& graph, const WeightMap& weight) 1573 : _graph(graph), _weight(weight), _matching(0), 1574 _node_potential(0), _blossom_potential(), _blossom_node_list(), 1575 _node_num(0), _blossom_num(0), 1576 1577 _blossom_index(0), _blossom_set(0), _blossom_data(0), 1578 _node_index(0), _node_heap_index(0), _node_data(0), 1579 _tree_set_index(0), _tree_set(0), 1580 1581 _delta1_index(0), _delta1(0), 1582 _delta2_index(0), _delta2(0), 1583 _delta3_index(0), _delta3(0), 1584 _delta4_index(0), _delta4(0), 1585 1586 _delta_sum(), _unmatched(0), 1587 1588 _fractional(0) 1589 {} 1590 ~MaxWeightedMatching()1591 ~MaxWeightedMatching() { 1592 destroyStructures(); 1593 if (_fractional) { 1594 delete _fractional; 1595 } 1596 } 1597 1598 /// \name Execution Control 1599 /// The simplest way to execute the algorithm is to use the 1600 /// \ref run() member function. 1601 1602 ///@{ 1603 1604 /// \brief Initialize the algorithm 1605 /// 1606 /// This function initializes the algorithm. init()1607 void init() { 1608 createStructures(); 1609 1610 _blossom_node_list.clear(); 1611 _blossom_potential.clear(); 1612 1613 for (ArcIt e(_graph); e != INVALID; ++e) { 1614 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP; 1615 } 1616 for (NodeIt n(_graph); n != INVALID; ++n) { 1617 (*_delta1_index)[n] = _delta1->PRE_HEAP; 1618 } 1619 for (EdgeIt e(_graph); e != INVALID; ++e) { 1620 (*_delta3_index)[e] = _delta3->PRE_HEAP; 1621 } 1622 for (int i = 0; i < _blossom_num; ++i) { 1623 (*_delta2_index)[i] = _delta2->PRE_HEAP; 1624 (*_delta4_index)[i] = _delta4->PRE_HEAP; 1625 } 1626 1627 _unmatched = _node_num; 1628 1629 _delta1->clear(); 1630 _delta2->clear(); 1631 _delta3->clear(); 1632 _delta4->clear(); 1633 _blossom_set->clear(); 1634 _tree_set->clear(); 1635 1636 int index = 0; 1637 for (NodeIt n(_graph); n != INVALID; ++n) { 1638 Value max = 0; 1639 for (OutArcIt e(_graph, n); e != INVALID; ++e) { 1640 if (_graph.target(e) == n) continue; 1641 if ((dualScale * _weight[e]) / 2 > max) { 1642 max = (dualScale * _weight[e]) / 2; 1643 } 1644 } 1645 (*_node_index)[n] = index; 1646 (*_node_data)[index].heap_index.clear(); 1647 (*_node_data)[index].heap.clear(); 1648 (*_node_data)[index].pot = max; 1649 _delta1->push(n, max); 1650 int blossom = 1651 _blossom_set->insert(n, std::numeric_limits<Value>::max()); 1652 1653 _tree_set->insert(blossom); 1654 1655 (*_blossom_data)[blossom].status = EVEN; 1656 (*_blossom_data)[blossom].pred = INVALID; 1657 (*_blossom_data)[blossom].next = INVALID; 1658 (*_blossom_data)[blossom].pot = 0; 1659 (*_blossom_data)[blossom].offset = 0; 1660 ++index; 1661 } 1662 for (EdgeIt e(_graph); e != INVALID; ++e) { 1663 int si = (*_node_index)[_graph.u(e)]; 1664 int ti = (*_node_index)[_graph.v(e)]; 1665 if (_graph.u(e) != _graph.v(e)) { 1666 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - 1667 dualScale * _weight[e]) / 2); 1668 } 1669 } 1670 } 1671 1672 /// \brief Initialize the algorithm with fractional matching 1673 /// 1674 /// This function initializes the algorithm with a fractional 1675 /// matching. This initialization is also called jumpstart heuristic. fractionalInit()1676 void fractionalInit() { 1677 createStructures(); 1678 1679 _blossom_node_list.clear(); 1680 _blossom_potential.clear(); 1681 1682 if (_fractional == 0) { 1683 _fractional = new FractionalMatching(_graph, _weight, false); 1684 } 1685 _fractional->run(); 1686 1687 for (ArcIt e(_graph); e != INVALID; ++e) { 1688 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP; 1689 } 1690 for (NodeIt n(_graph); n != INVALID; ++n) { 1691 (*_delta1_index)[n] = _delta1->PRE_HEAP; 1692 } 1693 for (EdgeIt e(_graph); e != INVALID; ++e) { 1694 (*_delta3_index)[e] = _delta3->PRE_HEAP; 1695 } 1696 for (int i = 0; i < _blossom_num; ++i) { 1697 (*_delta2_index)[i] = _delta2->PRE_HEAP; 1698 (*_delta4_index)[i] = _delta4->PRE_HEAP; 1699 } 1700 1701 _unmatched = 0; 1702 1703 _delta1->clear(); 1704 _delta2->clear(); 1705 _delta3->clear(); 1706 _delta4->clear(); 1707 _blossom_set->clear(); 1708 _tree_set->clear(); 1709 1710 int index = 0; 1711 for (NodeIt n(_graph); n != INVALID; ++n) { 1712 Value pot = _fractional->nodeValue(n); 1713 (*_node_index)[n] = index; 1714 (*_node_data)[index].pot = pot; 1715 (*_node_data)[index].heap_index.clear(); 1716 (*_node_data)[index].heap.clear(); 1717 int blossom = 1718 _blossom_set->insert(n, std::numeric_limits<Value>::max()); 1719 1720 (*_blossom_data)[blossom].status = MATCHED; 1721 (*_blossom_data)[blossom].pred = INVALID; 1722 (*_blossom_data)[blossom].next = _fractional->matching(n); 1723 if (_fractional->matching(n) == INVALID) { 1724 (*_blossom_data)[blossom].base = n; 1725 } 1726 (*_blossom_data)[blossom].pot = 0; 1727 (*_blossom_data)[blossom].offset = 0; 1728 ++index; 1729 } 1730 1731 typename Graph::template NodeMap<bool> processed(_graph, false); 1732 for (NodeIt n(_graph); n != INVALID; ++n) { 1733 if (processed[n]) continue; 1734 processed[n] = true; 1735 if (_fractional->matching(n) == INVALID) continue; 1736 int num = 1; 1737 Node v = _graph.target(_fractional->matching(n)); 1738 while (n != v) { 1739 processed[v] = true; 1740 v = _graph.target(_fractional->matching(v)); 1741 ++num; 1742 } 1743 1744 if (num % 2 == 1) { 1745 std::vector<int> subblossoms(num); 1746 1747 subblossoms[--num] = _blossom_set->find(n); 1748 _delta1->push(n, _fractional->nodeValue(n)); 1749 v = _graph.target(_fractional->matching(n)); 1750 while (n != v) { 1751 subblossoms[--num] = _blossom_set->find(v); 1752 _delta1->push(v, _fractional->nodeValue(v)); 1753 v = _graph.target(_fractional->matching(v)); 1754 } 1755 1756 int surface = 1757 _blossom_set->join(subblossoms.begin(), subblossoms.end()); 1758 (*_blossom_data)[surface].status = EVEN; 1759 (*_blossom_data)[surface].pred = INVALID; 1760 (*_blossom_data)[surface].next = INVALID; 1761 (*_blossom_data)[surface].pot = 0; 1762 (*_blossom_data)[surface].offset = 0; 1763 1764 _tree_set->insert(surface); 1765 ++_unmatched; 1766 } 1767 } 1768 1769 for (EdgeIt e(_graph); e != INVALID; ++e) { 1770 int si = (*_node_index)[_graph.u(e)]; 1771 int sb = _blossom_set->find(_graph.u(e)); 1772 int ti = (*_node_index)[_graph.v(e)]; 1773 int tb = _blossom_set->find(_graph.v(e)); 1774 if ((*_blossom_data)[sb].status == EVEN && 1775 (*_blossom_data)[tb].status == EVEN && sb != tb) { 1776 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - 1777 dualScale * _weight[e]) / 2); 1778 } 1779 } 1780 1781 for (NodeIt n(_graph); n != INVALID; ++n) { 1782 int nb = _blossom_set->find(n); 1783 if ((*_blossom_data)[nb].status != MATCHED) continue; 1784 int ni = (*_node_index)[n]; 1785 1786 for (OutArcIt e(_graph, n); e != INVALID; ++e) { 1787 Node v = _graph.target(e); 1788 int vb = _blossom_set->find(v); 1789 int vi = (*_node_index)[v]; 1790 1791 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 1792 dualScale * _weight[e]; 1793 1794 if ((*_blossom_data)[vb].status == EVEN) { 1795 1796 int vt = _tree_set->find(vb); 1797 1798 typename std::map<int, Arc>::iterator it = 1799 (*_node_data)[ni].heap_index.find(vt); 1800 1801 if (it != (*_node_data)[ni].heap_index.end()) { 1802 if ((*_node_data)[ni].heap[it->second] > rw) { 1803 (*_node_data)[ni].heap.replace(it->second, e); 1804 (*_node_data)[ni].heap.decrease(e, rw); 1805 it->second = e; 1806 } 1807 } else { 1808 (*_node_data)[ni].heap.push(e, rw); 1809 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e)); 1810 } 1811 } 1812 } 1813 1814 if (!(*_node_data)[ni].heap.empty()) { 1815 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); 1816 _delta2->push(nb, _blossom_set->classPrio(nb)); 1817 } 1818 } 1819 } 1820 1821 /// \brief Start the algorithm 1822 /// 1823 /// This function starts the algorithm. 1824 /// 1825 /// \pre \ref init() or \ref fractionalInit() must be called 1826 /// before using this function. start()1827 void start() { 1828 enum OpType { 1829 D1, D2, D3, D4 1830 }; 1831 1832 while (_unmatched > 0) { 1833 Value d1 = !_delta1->empty() ? 1834 _delta1->prio() : std::numeric_limits<Value>::max(); 1835 1836 Value d2 = !_delta2->empty() ? 1837 _delta2->prio() : std::numeric_limits<Value>::max(); 1838 1839 Value d3 = !_delta3->empty() ? 1840 _delta3->prio() : std::numeric_limits<Value>::max(); 1841 1842 Value d4 = !_delta4->empty() ? 1843 _delta4->prio() : std::numeric_limits<Value>::max(); 1844 1845 _delta_sum = d3; OpType ot = D3; 1846 if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; } 1847 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } 1848 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } 1849 1850 switch (ot) { 1851 case D1: 1852 { 1853 Node n = _delta1->top(); 1854 unmatchNode(n); 1855 --_unmatched; 1856 } 1857 break; 1858 case D2: 1859 { 1860 int blossom = _delta2->top(); 1861 Node n = _blossom_set->classTop(blossom); 1862 Arc a = (*_node_data)[(*_node_index)[n]].heap.top(); 1863 if ((*_blossom_data)[blossom].next == INVALID) { 1864 augmentOnArc(a); 1865 --_unmatched; 1866 } else { 1867 extendOnArc(a); 1868 } 1869 } 1870 break; 1871 case D3: 1872 { 1873 Edge e = _delta3->top(); 1874 1875 int left_blossom = _blossom_set->find(_graph.u(e)); 1876 int right_blossom = _blossom_set->find(_graph.v(e)); 1877 1878 if (left_blossom == right_blossom) { 1879 _delta3->pop(); 1880 } else { 1881 int left_tree = _tree_set->find(left_blossom); 1882 int right_tree = _tree_set->find(right_blossom); 1883 1884 if (left_tree == right_tree) { 1885 shrinkOnEdge(e, left_tree); 1886 } else { 1887 augmentOnEdge(e); 1888 _unmatched -= 2; 1889 } 1890 } 1891 } break; 1892 case D4: 1893 splitBlossom(_delta4->top()); 1894 break; 1895 } 1896 } 1897 extractMatching(); 1898 } 1899 1900 /// \brief Run the algorithm. 1901 /// 1902 /// This method runs the \c %MaxWeightedMatching algorithm. 1903 /// 1904 /// \note mwm.run() is just a shortcut of the following code. 1905 /// \code 1906 /// mwm.fractionalInit(); 1907 /// mwm.start(); 1908 /// \endcode run()1909 void run() { 1910 fractionalInit(); 1911 start(); 1912 } 1913 1914 /// @} 1915 1916 /// \name Primal Solution 1917 /// Functions to get the primal solution, i.e. the maximum weighted 1918 /// matching.\n 1919 /// Either \ref run() or \ref start() function should be called before 1920 /// using them. 1921 1922 /// @{ 1923 1924 /// \brief Return the weight of the matching. 1925 /// 1926 /// This function returns the weight of the found matching. 1927 /// 1928 /// \pre Either run() or start() must be called before using this function. matchingWeight()1929 Value matchingWeight() const { 1930 Value sum = 0; 1931 for (NodeIt n(_graph); n != INVALID; ++n) { 1932 if ((*_matching)[n] != INVALID) { 1933 sum += _weight[(*_matching)[n]]; 1934 } 1935 } 1936 return sum / 2; 1937 } 1938 1939 /// \brief Return the size (cardinality) of the matching. 1940 /// 1941 /// This function returns the size (cardinality) of the found matching. 1942 /// 1943 /// \pre Either run() or start() must be called before using this function. matchingSize()1944 int matchingSize() const { 1945 int num = 0; 1946 for (NodeIt n(_graph); n != INVALID; ++n) { 1947 if ((*_matching)[n] != INVALID) { 1948 ++num; 1949 } 1950 } 1951 return num /= 2; 1952 } 1953 1954 /// \brief Return \c true if the given edge is in the matching. 1955 /// 1956 /// This function returns \c true if the given edge is in the found 1957 /// matching. 1958 /// 1959 /// \pre Either run() or start() must be called before using this function. matching(const Edge & edge)1960 bool matching(const Edge& edge) const { 1961 return edge == (*_matching)[_graph.u(edge)]; 1962 } 1963 1964 /// \brief Return the matching arc (or edge) incident to the given node. 1965 /// 1966 /// This function returns the matching arc (or edge) incident to the 1967 /// given node in the found matching or \c INVALID if the node is 1968 /// not covered by the matching. 1969 /// 1970 /// \pre Either run() or start() must be called before using this function. matching(const Node & node)1971 Arc matching(const Node& node) const { 1972 return (*_matching)[node]; 1973 } 1974 1975 /// \brief Return a const reference to the matching map. 1976 /// 1977 /// This function returns a const reference to a node map that stores 1978 /// the matching arc (or edge) incident to each node. matchingMap()1979 const MatchingMap& matchingMap() const { 1980 return *_matching; 1981 } 1982 1983 /// \brief Return the mate of the given node. 1984 /// 1985 /// This function returns the mate of the given node in the found 1986 /// matching or \c INVALID if the node is not covered by the matching. 1987 /// 1988 /// \pre Either run() or start() must be called before using this function. mate(const Node & node)1989 Node mate(const Node& node) const { 1990 return (*_matching)[node] != INVALID ? 1991 _graph.target((*_matching)[node]) : INVALID; 1992 } 1993 1994 /// @} 1995 1996 /// \name Dual Solution 1997 /// Functions to get the dual solution.\n 1998 /// Either \ref run() or \ref start() function should be called before 1999 /// using them. 2000 2001 /// @{ 2002 2003 /// \brief Return the value of the dual solution. 2004 /// 2005 /// This function returns the value of the dual solution. 2006 /// It should be equal to the primal value scaled by \ref dualScale 2007 /// "dual scale". 2008 /// 2009 /// \pre Either run() or start() must be called before using this function. dualValue()2010 Value dualValue() const { 2011 Value sum = 0; 2012 for (NodeIt n(_graph); n != INVALID; ++n) { 2013 sum += nodeValue(n); 2014 } 2015 for (int i = 0; i < blossomNum(); ++i) { 2016 sum += blossomValue(i) * (blossomSize(i) / 2); 2017 } 2018 return sum; 2019 } 2020 2021 /// \brief Return the dual value (potential) of the given node. 2022 /// 2023 /// This function returns the dual value (potential) of the given node. 2024 /// 2025 /// \pre Either run() or start() must be called before using this function. nodeValue(const Node & n)2026 Value nodeValue(const Node& n) const { 2027 return (*_node_potential)[n]; 2028 } 2029 2030 /// \brief Return the number of the blossoms in the basis. 2031 /// 2032 /// This function returns the number of the blossoms in the basis. 2033 /// 2034 /// \pre Either run() or start() must be called before using this function. 2035 /// \see BlossomIt blossomNum()2036 int blossomNum() const { 2037 return _blossom_potential.size(); 2038 } 2039 2040 /// \brief Return the number of the nodes in the given blossom. 2041 /// 2042 /// This function returns the number of the nodes in the given blossom. 2043 /// 2044 /// \pre Either run() or start() must be called before using this function. 2045 /// \see BlossomIt blossomSize(int k)2046 int blossomSize(int k) const { 2047 return _blossom_potential[k].end - _blossom_potential[k].begin; 2048 } 2049 2050 /// \brief Return the dual value (ptential) of the given blossom. 2051 /// 2052 /// This function returns the dual value (ptential) of the given blossom. 2053 /// 2054 /// \pre Either run() or start() must be called before using this function. blossomValue(int k)2055 Value blossomValue(int k) const { 2056 return _blossom_potential[k].value; 2057 } 2058 2059 /// \brief Iterator for obtaining the nodes of a blossom. 2060 /// 2061 /// This class provides an iterator for obtaining the nodes of the 2062 /// given blossom. It lists a subset of the nodes. 2063 /// Before using this iterator, you must allocate a 2064 /// MaxWeightedMatching class and execute it. 2065 class BlossomIt { 2066 public: 2067 2068 /// \brief Constructor. 2069 /// 2070 /// Constructor to get the nodes of the given variable. 2071 /// 2072 /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or 2073 /// \ref MaxWeightedMatching::start() "algorithm.start()" must be 2074 /// called before initializing this iterator. BlossomIt(const MaxWeightedMatching & algorithm,int variable)2075 BlossomIt(const MaxWeightedMatching& algorithm, int variable) 2076 : _algorithm(&algorithm) 2077 { 2078 _index = _algorithm->_blossom_potential[variable].begin; 2079 _last = _algorithm->_blossom_potential[variable].end; 2080 } 2081 2082 /// \brief Conversion to \c Node. 2083 /// 2084 /// Conversion to \c Node. Node()2085 operator Node() const { 2086 return _algorithm->_blossom_node_list[_index]; 2087 } 2088 2089 /// \brief Increment operator. 2090 /// 2091 /// Increment operator. 2092 BlossomIt& operator++() { 2093 ++_index; 2094 return *this; 2095 } 2096 2097 /// \brief Validity checking 2098 /// 2099 /// Checks whether the iterator is invalid. 2100 bool operator==(Invalid) const { return _index == _last; } 2101 2102 /// \brief Validity checking 2103 /// 2104 /// Checks whether the iterator is valid. 2105 bool operator!=(Invalid) const { return _index != _last; } 2106 2107 private: 2108 const MaxWeightedMatching* _algorithm; 2109 int _last; 2110 int _index; 2111 }; 2112 2113 /// @} 2114 2115 }; 2116 2117 /// \ingroup matching 2118 /// 2119 /// \brief Weighted perfect matching in general graphs 2120 /// 2121 /// This class provides an efficient implementation of Edmond's 2122 /// maximum weighted perfect matching algorithm. The implementation 2123 /// is based on extensive use of priority queues and provides 2124 /// \f$O(nm\log n)\f$ time complexity. 2125 /// 2126 /// The maximum weighted perfect matching problem is to find a subset of 2127 /// the edges in an undirected graph with maximum overall weight for which 2128 /// each node has exactly one incident edge. 2129 /// It can be formulated with the following linear program. 2130 /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f] 2131 /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} 2132 \quad \forall B\in\mathcal{O}\f] */ 2133 /// \f[x_e \ge 0\quad \forall e\in E\f] 2134 /// \f[\max \sum_{e\in E}x_ew_e\f] 2135 /// where \f$\delta(X)\f$ is the set of edges incident to a node in 2136 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in 2137 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality 2138 /// subsets of the nodes. 2139 /// 2140 /// The algorithm calculates an optimal matching and a proof of the 2141 /// optimality. The solution of the dual problem can be used to check 2142 /// the result of the algorithm. The dual linear problem is the 2143 /// following. 2144 /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge 2145 w_{uv} \quad \forall uv\in E\f] */ 2146 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] 2147 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} 2148 \frac{\vert B \vert - 1}{2}z_B\f] */ 2149 /// 2150 /// The algorithm can be executed with the run() function. 2151 /// After it the matching (the primal solution) and the dual solution 2152 /// can be obtained using the query functions and the 2153 /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, 2154 /// which is able to iterate on the nodes of a blossom. 2155 /// If the value type is integer, then the dual solution is multiplied 2156 /// by \ref MaxWeightedMatching::dualScale "4". 2157 /// 2158 /// \tparam GR The undirected graph type the algorithm runs on. 2159 /// \tparam WM The type edge weight map. The default type is 2160 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>". 2161 #ifdef DOXYGEN 2162 template <typename GR, typename WM> 2163 #else 2164 template <typename GR, 2165 typename WM = typename GR::template EdgeMap<int> > 2166 #endif 2167 class MaxWeightedPerfectMatching { 2168 public: 2169 2170 /// The graph type of the algorithm 2171 typedef GR Graph; 2172 /// The type of the edge weight map 2173 typedef WM WeightMap; 2174 /// The value type of the edge weights 2175 typedef typename WeightMap::Value Value; 2176 2177 /// \brief Scaling factor for dual solution 2178 /// 2179 /// Scaling factor for dual solution, it is equal to 4 or 1 2180 /// according to the value type. 2181 static const int dualScale = 2182 std::numeric_limits<Value>::is_integer ? 4 : 1; 2183 2184 /// The type of the matching map 2185 typedef typename Graph::template NodeMap<typename Graph::Arc> 2186 MatchingMap; 2187 2188 private: 2189 2190 TEMPLATE_GRAPH_TYPEDEFS(Graph); 2191 2192 typedef typename Graph::template NodeMap<Value> NodePotential; 2193 typedef std::vector<Node> BlossomNodeList; 2194 2195 struct BlossomVariable { 2196 int begin, end; 2197 Value value; 2198 BlossomVariableBlossomVariable2199 BlossomVariable(int _begin, int _end, Value _value) 2200 : begin(_begin), end(_end), value(_value) {} 2201 2202 }; 2203 2204 typedef std::vector<BlossomVariable> BlossomPotential; 2205 2206 const Graph& _graph; 2207 const WeightMap& _weight; 2208 2209 MatchingMap* _matching; 2210 2211 NodePotential* _node_potential; 2212 2213 BlossomPotential _blossom_potential; 2214 BlossomNodeList _blossom_node_list; 2215 2216 int _node_num; 2217 int _blossom_num; 2218 2219 typedef RangeMap<int> IntIntMap; 2220 2221 enum Status { 2222 EVEN = -1, MATCHED = 0, ODD = 1 2223 }; 2224 2225 typedef HeapUnionFind<Value, IntNodeMap> BlossomSet; 2226 struct BlossomData { 2227 int tree; 2228 Status status; 2229 Arc pred, next; 2230 Value pot, offset; 2231 }; 2232 2233 IntNodeMap *_blossom_index; 2234 BlossomSet *_blossom_set; 2235 RangeMap<BlossomData>* _blossom_data; 2236 2237 IntNodeMap *_node_index; 2238 IntArcMap *_node_heap_index; 2239 2240 struct NodeData { 2241 NodeDataNodeData2242 NodeData(IntArcMap& node_heap_index) 2243 : heap(node_heap_index) {} 2244 2245 int blossom; 2246 Value pot; 2247 BinHeap<Value, IntArcMap> heap; 2248 std::map<int, Arc> heap_index; 2249 2250 int tree; 2251 }; 2252 2253 RangeMap<NodeData>* _node_data; 2254 2255 typedef ExtendFindEnum<IntIntMap> TreeSet; 2256 2257 IntIntMap *_tree_set_index; 2258 TreeSet *_tree_set; 2259 2260 IntIntMap *_delta2_index; 2261 BinHeap<Value, IntIntMap> *_delta2; 2262 2263 IntEdgeMap *_delta3_index; 2264 BinHeap<Value, IntEdgeMap> *_delta3; 2265 2266 IntIntMap *_delta4_index; 2267 BinHeap<Value, IntIntMap> *_delta4; 2268 2269 Value _delta_sum; 2270 int _unmatched; 2271 2272 typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap> 2273 FractionalMatching; 2274 FractionalMatching *_fractional; 2275 createStructures()2276 void createStructures() { 2277 _node_num = countNodes(_graph); 2278 _blossom_num = _node_num * 3 / 2; 2279 2280 if (!_matching) { 2281 _matching = new MatchingMap(_graph); 2282 } 2283 2284 if (!_node_potential) { 2285 _node_potential = new NodePotential(_graph); 2286 } 2287 2288 if (!_blossom_set) { 2289 _blossom_index = new IntNodeMap(_graph); 2290 _blossom_set = new BlossomSet(*_blossom_index); 2291 _blossom_data = new RangeMap<BlossomData>(_blossom_num); 2292 } else if (_blossom_data->size() != _blossom_num) { 2293 delete _blossom_data; 2294 _blossom_data = new RangeMap<BlossomData>(_blossom_num); 2295 } 2296 2297 if (!_node_index) { 2298 _node_index = new IntNodeMap(_graph); 2299 _node_heap_index = new IntArcMap(_graph); 2300 _node_data = new RangeMap<NodeData>(_node_num, 2301 NodeData(*_node_heap_index)); 2302 } else if (_node_data->size() != _node_num) { 2303 delete _node_data; 2304 _node_data = new RangeMap<NodeData>(_node_num, 2305 NodeData(*_node_heap_index)); 2306 } 2307 2308 if (!_tree_set) { 2309 _tree_set_index = new IntIntMap(_blossom_num); 2310 _tree_set = new TreeSet(*_tree_set_index); 2311 } else { 2312 _tree_set_index->resize(_blossom_num); 2313 } 2314 2315 if (!_delta2) { 2316 _delta2_index = new IntIntMap(_blossom_num); 2317 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index); 2318 } else { 2319 _delta2_index->resize(_blossom_num); 2320 } 2321 2322 if (!_delta3) { 2323 _delta3_index = new IntEdgeMap(_graph); 2324 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index); 2325 } 2326 2327 if (!_delta4) { 2328 _delta4_index = new IntIntMap(_blossom_num); 2329 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index); 2330 } else { 2331 _delta4_index->resize(_blossom_num); 2332 } 2333 } 2334 destroyStructures()2335 void destroyStructures() { 2336 if (_matching) { 2337 delete _matching; 2338 } 2339 if (_node_potential) { 2340 delete _node_potential; 2341 } 2342 if (_blossom_set) { 2343 delete _blossom_index; 2344 delete _blossom_set; 2345 delete _blossom_data; 2346 } 2347 2348 if (_node_index) { 2349 delete _node_index; 2350 delete _node_heap_index; 2351 delete _node_data; 2352 } 2353 2354 if (_tree_set) { 2355 delete _tree_set_index; 2356 delete _tree_set; 2357 } 2358 if (_delta2) { 2359 delete _delta2_index; 2360 delete _delta2; 2361 } 2362 if (_delta3) { 2363 delete _delta3_index; 2364 delete _delta3; 2365 } 2366 if (_delta4) { 2367 delete _delta4_index; 2368 delete _delta4; 2369 } 2370 } 2371 matchedToEven(int blossom,int tree)2372 void matchedToEven(int blossom, int tree) { 2373 if (_delta2->state(blossom) == _delta2->IN_HEAP) { 2374 _delta2->erase(blossom); 2375 } 2376 2377 if (!_blossom_set->trivial(blossom)) { 2378 (*_blossom_data)[blossom].pot -= 2379 2 * (_delta_sum - (*_blossom_data)[blossom].offset); 2380 } 2381 2382 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 2383 n != INVALID; ++n) { 2384 2385 _blossom_set->increase(n, std::numeric_limits<Value>::max()); 2386 int ni = (*_node_index)[n]; 2387 2388 (*_node_data)[ni].heap.clear(); 2389 (*_node_data)[ni].heap_index.clear(); 2390 2391 (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset; 2392 2393 for (InArcIt e(_graph, n); e != INVALID; ++e) { 2394 Node v = _graph.source(e); 2395 int vb = _blossom_set->find(v); 2396 int vi = (*_node_index)[v]; 2397 2398 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 2399 dualScale * _weight[e]; 2400 2401 if ((*_blossom_data)[vb].status == EVEN) { 2402 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { 2403 _delta3->push(e, rw / 2); 2404 } 2405 } else { 2406 typename std::map<int, Arc>::iterator it = 2407 (*_node_data)[vi].heap_index.find(tree); 2408 2409 if (it != (*_node_data)[vi].heap_index.end()) { 2410 if ((*_node_data)[vi].heap[it->second] > rw) { 2411 (*_node_data)[vi].heap.replace(it->second, e); 2412 (*_node_data)[vi].heap.decrease(e, rw); 2413 it->second = e; 2414 } 2415 } else { 2416 (*_node_data)[vi].heap.push(e, rw); 2417 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); 2418 } 2419 2420 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { 2421 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); 2422 2423 if ((*_blossom_data)[vb].status == MATCHED) { 2424 if (_delta2->state(vb) != _delta2->IN_HEAP) { 2425 _delta2->push(vb, _blossom_set->classPrio(vb) - 2426 (*_blossom_data)[vb].offset); 2427 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - 2428 (*_blossom_data)[vb].offset){ 2429 _delta2->decrease(vb, _blossom_set->classPrio(vb) - 2430 (*_blossom_data)[vb].offset); 2431 } 2432 } 2433 } 2434 } 2435 } 2436 } 2437 (*_blossom_data)[blossom].offset = 0; 2438 } 2439 matchedToOdd(int blossom)2440 void matchedToOdd(int blossom) { 2441 if (_delta2->state(blossom) == _delta2->IN_HEAP) { 2442 _delta2->erase(blossom); 2443 } 2444 (*_blossom_data)[blossom].offset += _delta_sum; 2445 if (!_blossom_set->trivial(blossom)) { 2446 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + 2447 (*_blossom_data)[blossom].offset); 2448 } 2449 } 2450 evenToMatched(int blossom,int tree)2451 void evenToMatched(int blossom, int tree) { 2452 if (!_blossom_set->trivial(blossom)) { 2453 (*_blossom_data)[blossom].pot += 2 * _delta_sum; 2454 } 2455 2456 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 2457 n != INVALID; ++n) { 2458 int ni = (*_node_index)[n]; 2459 (*_node_data)[ni].pot -= _delta_sum; 2460 2461 for (InArcIt e(_graph, n); e != INVALID; ++e) { 2462 Node v = _graph.source(e); 2463 int vb = _blossom_set->find(v); 2464 int vi = (*_node_index)[v]; 2465 2466 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 2467 dualScale * _weight[e]; 2468 2469 if (vb == blossom) { 2470 if (_delta3->state(e) == _delta3->IN_HEAP) { 2471 _delta3->erase(e); 2472 } 2473 } else if ((*_blossom_data)[vb].status == EVEN) { 2474 2475 if (_delta3->state(e) == _delta3->IN_HEAP) { 2476 _delta3->erase(e); 2477 } 2478 2479 int vt = _tree_set->find(vb); 2480 2481 if (vt != tree) { 2482 2483 Arc r = _graph.oppositeArc(e); 2484 2485 typename std::map<int, Arc>::iterator it = 2486 (*_node_data)[ni].heap_index.find(vt); 2487 2488 if (it != (*_node_data)[ni].heap_index.end()) { 2489 if ((*_node_data)[ni].heap[it->second] > rw) { 2490 (*_node_data)[ni].heap.replace(it->second, r); 2491 (*_node_data)[ni].heap.decrease(r, rw); 2492 it->second = r; 2493 } 2494 } else { 2495 (*_node_data)[ni].heap.push(r, rw); 2496 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); 2497 } 2498 2499 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { 2500 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); 2501 2502 if (_delta2->state(blossom) != _delta2->IN_HEAP) { 2503 _delta2->push(blossom, _blossom_set->classPrio(blossom) - 2504 (*_blossom_data)[blossom].offset); 2505 } else if ((*_delta2)[blossom] > 2506 _blossom_set->classPrio(blossom) - 2507 (*_blossom_data)[blossom].offset){ 2508 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - 2509 (*_blossom_data)[blossom].offset); 2510 } 2511 } 2512 } 2513 } else { 2514 2515 typename std::map<int, Arc>::iterator it = 2516 (*_node_data)[vi].heap_index.find(tree); 2517 2518 if (it != (*_node_data)[vi].heap_index.end()) { 2519 (*_node_data)[vi].heap.erase(it->second); 2520 (*_node_data)[vi].heap_index.erase(it); 2521 if ((*_node_data)[vi].heap.empty()) { 2522 _blossom_set->increase(v, std::numeric_limits<Value>::max()); 2523 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { 2524 _blossom_set->increase(v, (*_node_data)[vi].heap.prio()); 2525 } 2526 2527 if ((*_blossom_data)[vb].status == MATCHED) { 2528 if (_blossom_set->classPrio(vb) == 2529 std::numeric_limits<Value>::max()) { 2530 _delta2->erase(vb); 2531 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - 2532 (*_blossom_data)[vb].offset) { 2533 _delta2->increase(vb, _blossom_set->classPrio(vb) - 2534 (*_blossom_data)[vb].offset); 2535 } 2536 } 2537 } 2538 } 2539 } 2540 } 2541 } 2542 oddToMatched(int blossom)2543 void oddToMatched(int blossom) { 2544 (*_blossom_data)[blossom].offset -= _delta_sum; 2545 2546 if (_blossom_set->classPrio(blossom) != 2547 std::numeric_limits<Value>::max()) { 2548 _delta2->push(blossom, _blossom_set->classPrio(blossom) - 2549 (*_blossom_data)[blossom].offset); 2550 } 2551 2552 if (!_blossom_set->trivial(blossom)) { 2553 _delta4->erase(blossom); 2554 } 2555 } 2556 oddToEven(int blossom,int tree)2557 void oddToEven(int blossom, int tree) { 2558 if (!_blossom_set->trivial(blossom)) { 2559 _delta4->erase(blossom); 2560 (*_blossom_data)[blossom].pot -= 2561 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); 2562 } 2563 2564 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 2565 n != INVALID; ++n) { 2566 int ni = (*_node_index)[n]; 2567 2568 _blossom_set->increase(n, std::numeric_limits<Value>::max()); 2569 2570 (*_node_data)[ni].heap.clear(); 2571 (*_node_data)[ni].heap_index.clear(); 2572 (*_node_data)[ni].pot += 2573 2 * _delta_sum - (*_blossom_data)[blossom].offset; 2574 2575 for (InArcIt e(_graph, n); e != INVALID; ++e) { 2576 Node v = _graph.source(e); 2577 int vb = _blossom_set->find(v); 2578 int vi = (*_node_index)[v]; 2579 2580 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 2581 dualScale * _weight[e]; 2582 2583 if ((*_blossom_data)[vb].status == EVEN) { 2584 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { 2585 _delta3->push(e, rw / 2); 2586 } 2587 } else { 2588 2589 typename std::map<int, Arc>::iterator it = 2590 (*_node_data)[vi].heap_index.find(tree); 2591 2592 if (it != (*_node_data)[vi].heap_index.end()) { 2593 if ((*_node_data)[vi].heap[it->second] > rw) { 2594 (*_node_data)[vi].heap.replace(it->second, e); 2595 (*_node_data)[vi].heap.decrease(e, rw); 2596 it->second = e; 2597 } 2598 } else { 2599 (*_node_data)[vi].heap.push(e, rw); 2600 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); 2601 } 2602 2603 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { 2604 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); 2605 2606 if ((*_blossom_data)[vb].status == MATCHED) { 2607 if (_delta2->state(vb) != _delta2->IN_HEAP) { 2608 _delta2->push(vb, _blossom_set->classPrio(vb) - 2609 (*_blossom_data)[vb].offset); 2610 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - 2611 (*_blossom_data)[vb].offset) { 2612 _delta2->decrease(vb, _blossom_set->classPrio(vb) - 2613 (*_blossom_data)[vb].offset); 2614 } 2615 } 2616 } 2617 } 2618 } 2619 } 2620 (*_blossom_data)[blossom].offset = 0; 2621 } 2622 alternatePath(int even,int tree)2623 void alternatePath(int even, int tree) { 2624 int odd; 2625 2626 evenToMatched(even, tree); 2627 (*_blossom_data)[even].status = MATCHED; 2628 2629 while ((*_blossom_data)[even].pred != INVALID) { 2630 odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred)); 2631 (*_blossom_data)[odd].status = MATCHED; 2632 oddToMatched(odd); 2633 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred; 2634 2635 even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred)); 2636 (*_blossom_data)[even].status = MATCHED; 2637 evenToMatched(even, tree); 2638 (*_blossom_data)[even].next = 2639 _graph.oppositeArc((*_blossom_data)[odd].pred); 2640 } 2641 2642 } 2643 destroyTree(int tree)2644 void destroyTree(int tree) { 2645 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) { 2646 if ((*_blossom_data)[b].status == EVEN) { 2647 (*_blossom_data)[b].status = MATCHED; 2648 evenToMatched(b, tree); 2649 } else if ((*_blossom_data)[b].status == ODD) { 2650 (*_blossom_data)[b].status = MATCHED; 2651 oddToMatched(b); 2652 } 2653 } 2654 _tree_set->eraseClass(tree); 2655 } 2656 augmentOnEdge(const Edge & edge)2657 void augmentOnEdge(const Edge& edge) { 2658 2659 int left = _blossom_set->find(_graph.u(edge)); 2660 int right = _blossom_set->find(_graph.v(edge)); 2661 2662 int left_tree = _tree_set->find(left); 2663 alternatePath(left, left_tree); 2664 destroyTree(left_tree); 2665 2666 int right_tree = _tree_set->find(right); 2667 alternatePath(right, right_tree); 2668 destroyTree(right_tree); 2669 2670 (*_blossom_data)[left].next = _graph.direct(edge, true); 2671 (*_blossom_data)[right].next = _graph.direct(edge, false); 2672 } 2673 extendOnArc(const Arc & arc)2674 void extendOnArc(const Arc& arc) { 2675 int base = _blossom_set->find(_graph.target(arc)); 2676 int tree = _tree_set->find(base); 2677 2678 int odd = _blossom_set->find(_graph.source(arc)); 2679 _tree_set->insert(odd, tree); 2680 (*_blossom_data)[odd].status = ODD; 2681 matchedToOdd(odd); 2682 (*_blossom_data)[odd].pred = arc; 2683 2684 int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next)); 2685 (*_blossom_data)[even].pred = (*_blossom_data)[even].next; 2686 _tree_set->insert(even, tree); 2687 (*_blossom_data)[even].status = EVEN; 2688 matchedToEven(even, tree); 2689 } 2690 shrinkOnEdge(const Edge & edge,int tree)2691 void shrinkOnEdge(const Edge& edge, int tree) { 2692 int nca = -1; 2693 std::vector<int> left_path, right_path; 2694 2695 { 2696 std::set<int> left_set, right_set; 2697 int left = _blossom_set->find(_graph.u(edge)); 2698 left_path.push_back(left); 2699 left_set.insert(left); 2700 2701 int right = _blossom_set->find(_graph.v(edge)); 2702 right_path.push_back(right); 2703 right_set.insert(right); 2704 2705 while (true) { 2706 2707 if ((*_blossom_data)[left].pred == INVALID) break; 2708 2709 left = 2710 _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); 2711 left_path.push_back(left); 2712 left = 2713 _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); 2714 left_path.push_back(left); 2715 2716 left_set.insert(left); 2717 2718 if (right_set.find(left) != right_set.end()) { 2719 nca = left; 2720 break; 2721 } 2722 2723 if ((*_blossom_data)[right].pred == INVALID) break; 2724 2725 right = 2726 _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); 2727 right_path.push_back(right); 2728 right = 2729 _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); 2730 right_path.push_back(right); 2731 2732 right_set.insert(right); 2733 2734 if (left_set.find(right) != left_set.end()) { 2735 nca = right; 2736 break; 2737 } 2738 2739 } 2740 2741 if (nca == -1) { 2742 if ((*_blossom_data)[left].pred == INVALID) { 2743 nca = right; 2744 while (left_set.find(nca) == left_set.end()) { 2745 nca = 2746 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); 2747 right_path.push_back(nca); 2748 nca = 2749 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); 2750 right_path.push_back(nca); 2751 } 2752 } else { 2753 nca = left; 2754 while (right_set.find(nca) == right_set.end()) { 2755 nca = 2756 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); 2757 left_path.push_back(nca); 2758 nca = 2759 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); 2760 left_path.push_back(nca); 2761 } 2762 } 2763 } 2764 } 2765 2766 std::vector<int> subblossoms; 2767 Arc prev; 2768 2769 prev = _graph.direct(edge, true); 2770 for (int i = 0; left_path[i] != nca; i += 2) { 2771 subblossoms.push_back(left_path[i]); 2772 (*_blossom_data)[left_path[i]].next = prev; 2773 _tree_set->erase(left_path[i]); 2774 2775 subblossoms.push_back(left_path[i + 1]); 2776 (*_blossom_data)[left_path[i + 1]].status = EVEN; 2777 oddToEven(left_path[i + 1], tree); 2778 _tree_set->erase(left_path[i + 1]); 2779 prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred); 2780 } 2781 2782 int k = 0; 2783 while (right_path[k] != nca) ++k; 2784 2785 subblossoms.push_back(nca); 2786 (*_blossom_data)[nca].next = prev; 2787 2788 for (int i = k - 2; i >= 0; i -= 2) { 2789 subblossoms.push_back(right_path[i + 1]); 2790 (*_blossom_data)[right_path[i + 1]].status = EVEN; 2791 oddToEven(right_path[i + 1], tree); 2792 _tree_set->erase(right_path[i + 1]); 2793 2794 (*_blossom_data)[right_path[i + 1]].next = 2795 (*_blossom_data)[right_path[i + 1]].pred; 2796 2797 subblossoms.push_back(right_path[i]); 2798 _tree_set->erase(right_path[i]); 2799 } 2800 2801 int surface = 2802 _blossom_set->join(subblossoms.begin(), subblossoms.end()); 2803 2804 for (int i = 0; i < int(subblossoms.size()); ++i) { 2805 if (!_blossom_set->trivial(subblossoms[i])) { 2806 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum; 2807 } 2808 (*_blossom_data)[subblossoms[i]].status = MATCHED; 2809 } 2810 2811 (*_blossom_data)[surface].pot = -2 * _delta_sum; 2812 (*_blossom_data)[surface].offset = 0; 2813 (*_blossom_data)[surface].status = EVEN; 2814 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred; 2815 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred; 2816 2817 _tree_set->insert(surface, tree); 2818 _tree_set->erase(nca); 2819 } 2820 splitBlossom(int blossom)2821 void splitBlossom(int blossom) { 2822 Arc next = (*_blossom_data)[blossom].next; 2823 Arc pred = (*_blossom_data)[blossom].pred; 2824 2825 int tree = _tree_set->find(blossom); 2826 2827 (*_blossom_data)[blossom].status = MATCHED; 2828 oddToMatched(blossom); 2829 if (_delta2->state(blossom) == _delta2->IN_HEAP) { 2830 _delta2->erase(blossom); 2831 } 2832 2833 std::vector<int> subblossoms; 2834 _blossom_set->split(blossom, std::back_inserter(subblossoms)); 2835 2836 Value offset = (*_blossom_data)[blossom].offset; 2837 int b = _blossom_set->find(_graph.source(pred)); 2838 int d = _blossom_set->find(_graph.source(next)); 2839 2840 int ib = -1, id = -1; 2841 for (int i = 0; i < int(subblossoms.size()); ++i) { 2842 if (subblossoms[i] == b) ib = i; 2843 if (subblossoms[i] == d) id = i; 2844 2845 (*_blossom_data)[subblossoms[i]].offset = offset; 2846 if (!_blossom_set->trivial(subblossoms[i])) { 2847 (*_blossom_data)[subblossoms[i]].pot -= 2 * offset; 2848 } 2849 if (_blossom_set->classPrio(subblossoms[i]) != 2850 std::numeric_limits<Value>::max()) { 2851 _delta2->push(subblossoms[i], 2852 _blossom_set->classPrio(subblossoms[i]) - 2853 (*_blossom_data)[subblossoms[i]].offset); 2854 } 2855 } 2856 2857 if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) { 2858 for (int i = (id + 1) % subblossoms.size(); 2859 i != ib; i = (i + 2) % subblossoms.size()) { 2860 int sb = subblossoms[i]; 2861 int tb = subblossoms[(i + 1) % subblossoms.size()]; 2862 (*_blossom_data)[sb].next = 2863 _graph.oppositeArc((*_blossom_data)[tb].next); 2864 } 2865 2866 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) { 2867 int sb = subblossoms[i]; 2868 int tb = subblossoms[(i + 1) % subblossoms.size()]; 2869 int ub = subblossoms[(i + 2) % subblossoms.size()]; 2870 2871 (*_blossom_data)[sb].status = ODD; 2872 matchedToOdd(sb); 2873 _tree_set->insert(sb, tree); 2874 (*_blossom_data)[sb].pred = pred; 2875 (*_blossom_data)[sb].next = 2876 _graph.oppositeArc((*_blossom_data)[tb].next); 2877 2878 pred = (*_blossom_data)[ub].next; 2879 2880 (*_blossom_data)[tb].status = EVEN; 2881 matchedToEven(tb, tree); 2882 _tree_set->insert(tb, tree); 2883 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next; 2884 } 2885 2886 (*_blossom_data)[subblossoms[id]].status = ODD; 2887 matchedToOdd(subblossoms[id]); 2888 _tree_set->insert(subblossoms[id], tree); 2889 (*_blossom_data)[subblossoms[id]].next = next; 2890 (*_blossom_data)[subblossoms[id]].pred = pred; 2891 2892 } else { 2893 2894 for (int i = (ib + 1) % subblossoms.size(); 2895 i != id; i = (i + 2) % subblossoms.size()) { 2896 int sb = subblossoms[i]; 2897 int tb = subblossoms[(i + 1) % subblossoms.size()]; 2898 (*_blossom_data)[sb].next = 2899 _graph.oppositeArc((*_blossom_data)[tb].next); 2900 } 2901 2902 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) { 2903 int sb = subblossoms[i]; 2904 int tb = subblossoms[(i + 1) % subblossoms.size()]; 2905 int ub = subblossoms[(i + 2) % subblossoms.size()]; 2906 2907 (*_blossom_data)[sb].status = ODD; 2908 matchedToOdd(sb); 2909 _tree_set->insert(sb, tree); 2910 (*_blossom_data)[sb].next = next; 2911 (*_blossom_data)[sb].pred = 2912 _graph.oppositeArc((*_blossom_data)[tb].next); 2913 2914 (*_blossom_data)[tb].status = EVEN; 2915 matchedToEven(tb, tree); 2916 _tree_set->insert(tb, tree); 2917 (*_blossom_data)[tb].pred = 2918 (*_blossom_data)[tb].next = 2919 _graph.oppositeArc((*_blossom_data)[ub].next); 2920 next = (*_blossom_data)[ub].next; 2921 } 2922 2923 (*_blossom_data)[subblossoms[ib]].status = ODD; 2924 matchedToOdd(subblossoms[ib]); 2925 _tree_set->insert(subblossoms[ib], tree); 2926 (*_blossom_data)[subblossoms[ib]].next = next; 2927 (*_blossom_data)[subblossoms[ib]].pred = pred; 2928 } 2929 _tree_set->erase(blossom); 2930 } 2931 extractBlossom(int blossom,const Node & base,const Arc & matching)2932 void extractBlossom(int blossom, const Node& base, const Arc& matching) { 2933 if (_blossom_set->trivial(blossom)) { 2934 int bi = (*_node_index)[base]; 2935 Value pot = (*_node_data)[bi].pot; 2936 2937 (*_matching)[base] = matching; 2938 _blossom_node_list.push_back(base); 2939 (*_node_potential)[base] = pot; 2940 } else { 2941 2942 Value pot = (*_blossom_data)[blossom].pot; 2943 int bn = _blossom_node_list.size(); 2944 2945 std::vector<int> subblossoms; 2946 _blossom_set->split(blossom, std::back_inserter(subblossoms)); 2947 int b = _blossom_set->find(base); 2948 int ib = -1; 2949 for (int i = 0; i < int(subblossoms.size()); ++i) { 2950 if (subblossoms[i] == b) { ib = i; break; } 2951 } 2952 2953 for (int i = 1; i < int(subblossoms.size()); i += 2) { 2954 int sb = subblossoms[(ib + i) % subblossoms.size()]; 2955 int tb = subblossoms[(ib + i + 1) % subblossoms.size()]; 2956 2957 Arc m = (*_blossom_data)[tb].next; 2958 extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m)); 2959 extractBlossom(tb, _graph.source(m), m); 2960 } 2961 extractBlossom(subblossoms[ib], base, matching); 2962 2963 int en = _blossom_node_list.size(); 2964 2965 _blossom_potential.push_back(BlossomVariable(bn, en, pot)); 2966 } 2967 } 2968 extractMatching()2969 void extractMatching() { 2970 std::vector<int> blossoms; 2971 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) { 2972 blossoms.push_back(c); 2973 } 2974 2975 for (int i = 0; i < int(blossoms.size()); ++i) { 2976 2977 Value offset = (*_blossom_data)[blossoms[i]].offset; 2978 (*_blossom_data)[blossoms[i]].pot += 2 * offset; 2979 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); 2980 n != INVALID; ++n) { 2981 (*_node_data)[(*_node_index)[n]].pot -= offset; 2982 } 2983 2984 Arc matching = (*_blossom_data)[blossoms[i]].next; 2985 Node base = _graph.source(matching); 2986 extractBlossom(blossoms[i], base, matching); 2987 } 2988 } 2989 2990 public: 2991 2992 /// \brief Constructor 2993 /// 2994 /// Constructor. MaxWeightedPerfectMatching(const Graph & graph,const WeightMap & weight)2995 MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight) 2996 : _graph(graph), _weight(weight), _matching(0), 2997 _node_potential(0), _blossom_potential(), _blossom_node_list(), 2998 _node_num(0), _blossom_num(0), 2999 3000 _blossom_index(0), _blossom_set(0), _blossom_data(0), 3001 _node_index(0), _node_heap_index(0), _node_data(0), 3002 _tree_set_index(0), _tree_set(0), 3003 3004 _delta2_index(0), _delta2(0), 3005 _delta3_index(0), _delta3(0), 3006 _delta4_index(0), _delta4(0), 3007 3008 _delta_sum(), _unmatched(0), 3009 3010 _fractional(0) 3011 {} 3012 ~MaxWeightedPerfectMatching()3013 ~MaxWeightedPerfectMatching() { 3014 destroyStructures(); 3015 if (_fractional) { 3016 delete _fractional; 3017 } 3018 } 3019 3020 /// \name Execution Control 3021 /// The simplest way to execute the algorithm is to use the 3022 /// \ref run() member function. 3023 3024 ///@{ 3025 3026 /// \brief Initialize the algorithm 3027 /// 3028 /// This function initializes the algorithm. init()3029 void init() { 3030 createStructures(); 3031 3032 _blossom_node_list.clear(); 3033 _blossom_potential.clear(); 3034 3035 for (ArcIt e(_graph); e != INVALID; ++e) { 3036 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP; 3037 } 3038 for (EdgeIt e(_graph); e != INVALID; ++e) { 3039 (*_delta3_index)[e] = _delta3->PRE_HEAP; 3040 } 3041 for (int i = 0; i < _blossom_num; ++i) { 3042 (*_delta2_index)[i] = _delta2->PRE_HEAP; 3043 (*_delta4_index)[i] = _delta4->PRE_HEAP; 3044 } 3045 3046 _unmatched = _node_num; 3047 3048 _delta2->clear(); 3049 _delta3->clear(); 3050 _delta4->clear(); 3051 _blossom_set->clear(); 3052 _tree_set->clear(); 3053 3054 int index = 0; 3055 for (NodeIt n(_graph); n != INVALID; ++n) { 3056 Value max = - std::numeric_limits<Value>::max(); 3057 for (OutArcIt e(_graph, n); e != INVALID; ++e) { 3058 if (_graph.target(e) == n) continue; 3059 if ((dualScale * _weight[e]) / 2 > max) { 3060 max = (dualScale * _weight[e]) / 2; 3061 } 3062 } 3063 (*_node_index)[n] = index; 3064 (*_node_data)[index].heap_index.clear(); 3065 (*_node_data)[index].heap.clear(); 3066 (*_node_data)[index].pot = max; 3067 int blossom = 3068 _blossom_set->insert(n, std::numeric_limits<Value>::max()); 3069 3070 _tree_set->insert(blossom); 3071 3072 (*_blossom_data)[blossom].status = EVEN; 3073 (*_blossom_data)[blossom].pred = INVALID; 3074 (*_blossom_data)[blossom].next = INVALID; 3075 (*_blossom_data)[blossom].pot = 0; 3076 (*_blossom_data)[blossom].offset = 0; 3077 ++index; 3078 } 3079 for (EdgeIt e(_graph); e != INVALID; ++e) { 3080 int si = (*_node_index)[_graph.u(e)]; 3081 int ti = (*_node_index)[_graph.v(e)]; 3082 if (_graph.u(e) != _graph.v(e)) { 3083 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - 3084 dualScale * _weight[e]) / 2); 3085 } 3086 } 3087 } 3088 3089 /// \brief Initialize the algorithm with fractional matching 3090 /// 3091 /// This function initializes the algorithm with a fractional 3092 /// matching. This initialization is also called jumpstart heuristic. fractionalInit()3093 void fractionalInit() { 3094 createStructures(); 3095 3096 _blossom_node_list.clear(); 3097 _blossom_potential.clear(); 3098 3099 if (_fractional == 0) { 3100 _fractional = new FractionalMatching(_graph, _weight, false); 3101 } 3102 if (!_fractional->run()) { 3103 _unmatched = -1; 3104 return; 3105 } 3106 3107 for (ArcIt e(_graph); e != INVALID; ++e) { 3108 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP; 3109 } 3110 for (EdgeIt e(_graph); e != INVALID; ++e) { 3111 (*_delta3_index)[e] = _delta3->PRE_HEAP; 3112 } 3113 for (int i = 0; i < _blossom_num; ++i) { 3114 (*_delta2_index)[i] = _delta2->PRE_HEAP; 3115 (*_delta4_index)[i] = _delta4->PRE_HEAP; 3116 } 3117 3118 _unmatched = 0; 3119 3120 _delta2->clear(); 3121 _delta3->clear(); 3122 _delta4->clear(); 3123 _blossom_set->clear(); 3124 _tree_set->clear(); 3125 3126 int index = 0; 3127 for (NodeIt n(_graph); n != INVALID; ++n) { 3128 Value pot = _fractional->nodeValue(n); 3129 (*_node_index)[n] = index; 3130 (*_node_data)[index].pot = pot; 3131 (*_node_data)[index].heap_index.clear(); 3132 (*_node_data)[index].heap.clear(); 3133 int blossom = 3134 _blossom_set->insert(n, std::numeric_limits<Value>::max()); 3135 3136 (*_blossom_data)[blossom].status = MATCHED; 3137 (*_blossom_data)[blossom].pred = INVALID; 3138 (*_blossom_data)[blossom].next = _fractional->matching(n); 3139 (*_blossom_data)[blossom].pot = 0; 3140 (*_blossom_data)[blossom].offset = 0; 3141 ++index; 3142 } 3143 3144 typename Graph::template NodeMap<bool> processed(_graph, false); 3145 for (NodeIt n(_graph); n != INVALID; ++n) { 3146 if (processed[n]) continue; 3147 processed[n] = true; 3148 if (_fractional->matching(n) == INVALID) continue; 3149 int num = 1; 3150 Node v = _graph.target(_fractional->matching(n)); 3151 while (n != v) { 3152 processed[v] = true; 3153 v = _graph.target(_fractional->matching(v)); 3154 ++num; 3155 } 3156 3157 if (num % 2 == 1) { 3158 std::vector<int> subblossoms(num); 3159 3160 subblossoms[--num] = _blossom_set->find(n); 3161 v = _graph.target(_fractional->matching(n)); 3162 while (n != v) { 3163 subblossoms[--num] = _blossom_set->find(v); 3164 v = _graph.target(_fractional->matching(v)); 3165 } 3166 3167 int surface = 3168 _blossom_set->join(subblossoms.begin(), subblossoms.end()); 3169 (*_blossom_data)[surface].status = EVEN; 3170 (*_blossom_data)[surface].pred = INVALID; 3171 (*_blossom_data)[surface].next = INVALID; 3172 (*_blossom_data)[surface].pot = 0; 3173 (*_blossom_data)[surface].offset = 0; 3174 3175 _tree_set->insert(surface); 3176 ++_unmatched; 3177 } 3178 } 3179 3180 for (EdgeIt e(_graph); e != INVALID; ++e) { 3181 int si = (*_node_index)[_graph.u(e)]; 3182 int sb = _blossom_set->find(_graph.u(e)); 3183 int ti = (*_node_index)[_graph.v(e)]; 3184 int tb = _blossom_set->find(_graph.v(e)); 3185 if ((*_blossom_data)[sb].status == EVEN && 3186 (*_blossom_data)[tb].status == EVEN && sb != tb) { 3187 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - 3188 dualScale * _weight[e]) / 2); 3189 } 3190 } 3191 3192 for (NodeIt n(_graph); n != INVALID; ++n) { 3193 int nb = _blossom_set->find(n); 3194 if ((*_blossom_data)[nb].status != MATCHED) continue; 3195 int ni = (*_node_index)[n]; 3196 3197 for (OutArcIt e(_graph, n); e != INVALID; ++e) { 3198 Node v = _graph.target(e); 3199 int vb = _blossom_set->find(v); 3200 int vi = (*_node_index)[v]; 3201 3202 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 3203 dualScale * _weight[e]; 3204 3205 if ((*_blossom_data)[vb].status == EVEN) { 3206 3207 int vt = _tree_set->find(vb); 3208 3209 typename std::map<int, Arc>::iterator it = 3210 (*_node_data)[ni].heap_index.find(vt); 3211 3212 if (it != (*_node_data)[ni].heap_index.end()) { 3213 if ((*_node_data)[ni].heap[it->second] > rw) { 3214 (*_node_data)[ni].heap.replace(it->second, e); 3215 (*_node_data)[ni].heap.decrease(e, rw); 3216 it->second = e; 3217 } 3218 } else { 3219 (*_node_data)[ni].heap.push(e, rw); 3220 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e)); 3221 } 3222 } 3223 } 3224 3225 if (!(*_node_data)[ni].heap.empty()) { 3226 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); 3227 _delta2->push(nb, _blossom_set->classPrio(nb)); 3228 } 3229 } 3230 } 3231 3232 /// \brief Start the algorithm 3233 /// 3234 /// This function starts the algorithm. 3235 /// 3236 /// \pre \ref init() or \ref fractionalInit() must be called before 3237 /// using this function. start()3238 bool start() { 3239 enum OpType { 3240 D2, D3, D4 3241 }; 3242 3243 if (_unmatched == -1) return false; 3244 3245 while (_unmatched > 0) { 3246 Value d2 = !_delta2->empty() ? 3247 _delta2->prio() : std::numeric_limits<Value>::max(); 3248 3249 Value d3 = !_delta3->empty() ? 3250 _delta3->prio() : std::numeric_limits<Value>::max(); 3251 3252 Value d4 = !_delta4->empty() ? 3253 _delta4->prio() : std::numeric_limits<Value>::max(); 3254 3255 _delta_sum = d3; OpType ot = D3; 3256 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } 3257 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } 3258 3259 if (_delta_sum == std::numeric_limits<Value>::max()) { 3260 return false; 3261 } 3262 3263 switch (ot) { 3264 case D2: 3265 { 3266 int blossom = _delta2->top(); 3267 Node n = _blossom_set->classTop(blossom); 3268 Arc e = (*_node_data)[(*_node_index)[n]].heap.top(); 3269 extendOnArc(e); 3270 } 3271 break; 3272 case D3: 3273 { 3274 Edge e = _delta3->top(); 3275 3276 int left_blossom = _blossom_set->find(_graph.u(e)); 3277 int right_blossom = _blossom_set->find(_graph.v(e)); 3278 3279 if (left_blossom == right_blossom) { 3280 _delta3->pop(); 3281 } else { 3282 int left_tree = _tree_set->find(left_blossom); 3283 int right_tree = _tree_set->find(right_blossom); 3284 3285 if (left_tree == right_tree) { 3286 shrinkOnEdge(e, left_tree); 3287 } else { 3288 augmentOnEdge(e); 3289 _unmatched -= 2; 3290 } 3291 } 3292 } break; 3293 case D4: 3294 splitBlossom(_delta4->top()); 3295 break; 3296 } 3297 } 3298 extractMatching(); 3299 return true; 3300 } 3301 3302 /// \brief Run the algorithm. 3303 /// 3304 /// This method runs the \c %MaxWeightedPerfectMatching algorithm. 3305 /// 3306 /// \note mwpm.run() is just a shortcut of the following code. 3307 /// \code 3308 /// mwpm.fractionalInit(); 3309 /// mwpm.start(); 3310 /// \endcode run()3311 bool run() { 3312 fractionalInit(); 3313 return start(); 3314 } 3315 3316 /// @} 3317 3318 /// \name Primal Solution 3319 /// Functions to get the primal solution, i.e. the maximum weighted 3320 /// perfect matching.\n 3321 /// Either \ref run() or \ref start() function should be called before 3322 /// using them. 3323 3324 /// @{ 3325 3326 /// \brief Return the weight of the matching. 3327 /// 3328 /// This function returns the weight of the found matching. 3329 /// 3330 /// \pre Either run() or start() must be called before using this function. matchingWeight()3331 Value matchingWeight() const { 3332 Value sum = 0; 3333 for (NodeIt n(_graph); n != INVALID; ++n) { 3334 if ((*_matching)[n] != INVALID) { 3335 sum += _weight[(*_matching)[n]]; 3336 } 3337 } 3338 return sum / 2; 3339 } 3340 3341 /// \brief Return \c true if the given edge is in the matching. 3342 /// 3343 /// This function returns \c true if the given edge is in the found 3344 /// matching. 3345 /// 3346 /// \pre Either run() or start() must be called before using this function. matching(const Edge & edge)3347 bool matching(const Edge& edge) const { 3348 return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge; 3349 } 3350 3351 /// \brief Return the matching arc (or edge) incident to the given node. 3352 /// 3353 /// This function returns the matching arc (or edge) incident to the 3354 /// given node in the found matching or \c INVALID if the node is 3355 /// not covered by the matching. 3356 /// 3357 /// \pre Either run() or start() must be called before using this function. matching(const Node & node)3358 Arc matching(const Node& node) const { 3359 return (*_matching)[node]; 3360 } 3361 3362 /// \brief Return a const reference to the matching map. 3363 /// 3364 /// This function returns a const reference to a node map that stores 3365 /// the matching arc (or edge) incident to each node. matchingMap()3366 const MatchingMap& matchingMap() const { 3367 return *_matching; 3368 } 3369 3370 /// \brief Return the mate of the given node. 3371 /// 3372 /// This function returns the mate of the given node in the found 3373 /// matching or \c INVALID if the node is not covered by the matching. 3374 /// 3375 /// \pre Either run() or start() must be called before using this function. mate(const Node & node)3376 Node mate(const Node& node) const { 3377 return _graph.target((*_matching)[node]); 3378 } 3379 3380 /// @} 3381 3382 /// \name Dual Solution 3383 /// Functions to get the dual solution.\n 3384 /// Either \ref run() or \ref start() function should be called before 3385 /// using them. 3386 3387 /// @{ 3388 3389 /// \brief Return the value of the dual solution. 3390 /// 3391 /// This function returns the value of the dual solution. 3392 /// It should be equal to the primal value scaled by \ref dualScale 3393 /// "dual scale". 3394 /// 3395 /// \pre Either run() or start() must be called before using this function. dualValue()3396 Value dualValue() const { 3397 Value sum = 0; 3398 for (NodeIt n(_graph); n != INVALID; ++n) { 3399 sum += nodeValue(n); 3400 } 3401 for (int i = 0; i < blossomNum(); ++i) { 3402 sum += blossomValue(i) * (blossomSize(i) / 2); 3403 } 3404 return sum; 3405 } 3406 3407 /// \brief Return the dual value (potential) of the given node. 3408 /// 3409 /// This function returns the dual value (potential) of the given node. 3410 /// 3411 /// \pre Either run() or start() must be called before using this function. nodeValue(const Node & n)3412 Value nodeValue(const Node& n) const { 3413 return (*_node_potential)[n]; 3414 } 3415 3416 /// \brief Return the number of the blossoms in the basis. 3417 /// 3418 /// This function returns the number of the blossoms in the basis. 3419 /// 3420 /// \pre Either run() or start() must be called before using this function. 3421 /// \see BlossomIt blossomNum()3422 int blossomNum() const { 3423 return _blossom_potential.size(); 3424 } 3425 3426 /// \brief Return the number of the nodes in the given blossom. 3427 /// 3428 /// This function returns the number of the nodes in the given blossom. 3429 /// 3430 /// \pre Either run() or start() must be called before using this function. 3431 /// \see BlossomIt blossomSize(int k)3432 int blossomSize(int k) const { 3433 return _blossom_potential[k].end - _blossom_potential[k].begin; 3434 } 3435 3436 /// \brief Return the dual value (ptential) of the given blossom. 3437 /// 3438 /// This function returns the dual value (ptential) of the given blossom. 3439 /// 3440 /// \pre Either run() or start() must be called before using this function. blossomValue(int k)3441 Value blossomValue(int k) const { 3442 return _blossom_potential[k].value; 3443 } 3444 3445 /// \brief Iterator for obtaining the nodes of a blossom. 3446 /// 3447 /// This class provides an iterator for obtaining the nodes of the 3448 /// given blossom. It lists a subset of the nodes. 3449 /// Before using this iterator, you must allocate a 3450 /// MaxWeightedPerfectMatching class and execute it. 3451 class BlossomIt { 3452 public: 3453 3454 /// \brief Constructor. 3455 /// 3456 /// Constructor to get the nodes of the given variable. 3457 /// 3458 /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" 3459 /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" 3460 /// must be called before initializing this iterator. BlossomIt(const MaxWeightedPerfectMatching & algorithm,int variable)3461 BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable) 3462 : _algorithm(&algorithm) 3463 { 3464 _index = _algorithm->_blossom_potential[variable].begin; 3465 _last = _algorithm->_blossom_potential[variable].end; 3466 } 3467 3468 /// \brief Conversion to \c Node. 3469 /// 3470 /// Conversion to \c Node. Node()3471 operator Node() const { 3472 return _algorithm->_blossom_node_list[_index]; 3473 } 3474 3475 /// \brief Increment operator. 3476 /// 3477 /// Increment operator. 3478 BlossomIt& operator++() { 3479 ++_index; 3480 return *this; 3481 } 3482 3483 /// \brief Validity checking 3484 /// 3485 /// This function checks whether the iterator is invalid. 3486 bool operator==(Invalid) const { return _index == _last; } 3487 3488 /// \brief Validity checking 3489 /// 3490 /// This function checks whether the iterator is valid. 3491 bool operator!=(Invalid) const { return _index != _last; } 3492 3493 private: 3494 const MaxWeightedPerfectMatching* _algorithm; 3495 int _last; 3496 int _index; 3497 }; 3498 3499 /// @} 3500 3501 }; 3502 3503 } //END OF NAMESPACE LEMON 3504 3505 #endif //LEMON_MATCHING_H 3506