Changeset 339:91d63b8b1a4c in lemon for lemon
 Timestamp:
 10/13/08 14:00:11 (13 years ago)
 Branch:
 default
 Phase:
 public
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/max_matching.h
r338 r339 32 32 ///\ingroup matching 33 33 ///\file 34 ///\brief Maximum matching algorithms in g raph.34 ///\brief Maximum matching algorithms in general graphs. 35 35 36 36 namespace lemon { 37 37 38 /// \ingroup matching38 /// \ingroup matching 39 39 /// 40 /// \brief Edmonds' alternating forest maximum matching algorithm.40 /// \brief Edmonds' alternating forest maximum matching algorithm. 41 41 /// 42 /// This class provides Edmonds' alternating forest matching43 /// algorithm. The starting matching (if any) can be passed to the44 /// algorithm using some of init functions.42 /// This class provides Edmonds' alternating forest matching 43 /// algorithm. The starting matching (if any) can be passed to the 44 /// algorithm using some of init functions. 45 45 /// 46 ///The dual side of a matching is a map of the nodes to 47 ///MaxMatching::DecompType, having values \c D, \c A and \c C 48 ///showing the GallaiEdmonds decomposition of the digraph. The nodes 49 ///in \c D induce a digraph with factorcritical components, the nodes 50 ///in \c A form the barrier, and the nodes in \c C induce a digraph 51 ///having a perfect matching. This decomposition can be attained by 52 ///calling \c decomposition() after running the algorithm. 46 /// The dual side of a matching is a map of the nodes to 47 /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c 48 /// MATCHED/C showing the GallaiEdmonds decomposition of the 49 /// graph. The nodes in \c EVEN/D induce a graph with 50 /// factorcritical components, the nodes in \c ODD/A form the 51 /// barrier, and the nodes in \c MATCHED/C induce a graph having a 52 /// perfect matching. The number of the fractor critical components 53 /// minus the number of barrier nodes is a lower bound on the 54 /// unmatched nodes, and if the matching is optimal this bound is 55 /// tight. This decomposition can be attained by calling \c 56 /// decomposition() after running the algorithm. 53 57 /// 54 /// \param Digraph The graph type the algorithm runs on.55 template <typename Graph>58 /// \param _Graph The graph type the algorithm runs on. 59 template <typename _Graph> 56 60 class MaxMatching { 57 58 protected: 61 public: 62 63 typedef _Graph Graph; 64 typedef typename Graph::template NodeMap<typename Graph::Arc> 65 MatchingMap; 66 67 ///\brief Indicates the GallaiEdmonds decomposition of the graph. 68 /// 69 ///Indicates the GallaiEdmonds decomposition of the graph, which 70 ///shows an upper bound on the size of a maximum matching. The 71 ///nodes with Status \c EVEN/D induce a graph with factorcritical 72 ///components, the nodes in \c ODD/A form the canonical barrier, 73 ///and the nodes in \c MATCHED/C induce a graph having a perfect 74 ///matching. 75 enum Status { 76 EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = 1, A = 1, UNMATCHED = 2 77 }; 78 79 typedef typename Graph::template NodeMap<Status> StatusMap; 80 81 private: 59 82 60 83 TEMPLATE_GRAPH_TYPEDEFS(Graph); 61 84 62 typedef typename Graph::template NodeMap<int> UFECrossRef; 63 typedef UnionFindEnum<UFECrossRef> UFE; 64 typedef std::vector<Node> NV; 65 66 typedef typename Graph::template NodeMap<int> EFECrossRef; 67 typedef ExtendFindEnum<EFECrossRef> EFE; 85 typedef UnionFindEnum<IntNodeMap> BlossomSet; 86 typedef ExtendFindEnum<IntNodeMap> TreeSet; 87 typedef RangeMap<Node> NodeIntMap; 88 typedef MatchingMap EarMap; 89 typedef std::vector<Node> NodeQueue; 90 91 const Graph& _graph; 92 MatchingMap* _matching; 93 StatusMap* _status; 94 95 EarMap* _ear; 96 97 IntNodeMap* _blossom_set_index; 98 BlossomSet* _blossom_set; 99 NodeIntMap* _blossom_rep; 100 101 IntNodeMap* _tree_set_index; 102 TreeSet* _tree_set; 103 104 NodeQueue _node_queue; 105 int _process, _postpone, _last; 106 107 int _node_num; 108 109 private: 110 111 void createStructures() { 112 _node_num = countNodes(_graph); 113 if (!_matching) { 114 _matching = new MatchingMap(_graph); 115 } 116 if (!_status) { 117 _status = new StatusMap(_graph); 118 } 119 if (!_ear) { 120 _ear = new EarMap(_graph); 121 } 122 if (!_blossom_set) { 123 _blossom_set_index = new IntNodeMap(_graph); 124 _blossom_set = new BlossomSet(*_blossom_set_index); 125 } 126 if (!_blossom_rep) { 127 _blossom_rep = new NodeIntMap(_node_num); 128 } 129 if (!_tree_set) { 130 _tree_set_index = new IntNodeMap(_graph); 131 _tree_set = new TreeSet(*_tree_set_index); 132 } 133 _node_queue.resize(_node_num); 134 } 135 136 void destroyStructures() { 137 if (_matching) { 138 delete _matching; 139 } 140 if (_status) { 141 delete _status; 142 } 143 if (_ear) { 144 delete _ear; 145 } 146 if (_blossom_set) { 147 delete _blossom_set; 148 delete _blossom_set_index; 149 } 150 if (_blossom_rep) { 151 delete _blossom_rep; 152 } 153 if (_tree_set) { 154 delete _tree_set_index; 155 delete _tree_set; 156 } 157 } 158 159 void processDense(const Node& n) { 160 _process = _postpone = _last = 0; 161 _node_queue[_last++] = n; 162 163 while (_process != _last) { 164 Node u = _node_queue[_process++]; 165 for (OutArcIt a(_graph, u); a != INVALID; ++a) { 166 Node v = _graph.target(a); 167 if ((*_status)[v] == MATCHED) { 168 extendOnArc(a); 169 } else if ((*_status)[v] == UNMATCHED) { 170 augmentOnArc(a); 171 return; 172 } 173 } 174 } 175 176 while (_postpone != _last) { 177 Node u = _node_queue[_postpone++]; 178 179 for (OutArcIt a(_graph, u); a != INVALID ; ++a) { 180 Node v = _graph.target(a); 181 182 if ((*_status)[v] == EVEN) { 183 if (_blossom_set>find(u) != _blossom_set>find(v)) { 184 shrinkOnEdge(a); 185 } 186 } 187 188 while (_process != _last) { 189 Node w = _node_queue[_process++]; 190 for (OutArcIt b(_graph, w); b != INVALID; ++b) { 191 Node x = _graph.target(b); 192 if ((*_status)[x] == MATCHED) { 193 extendOnArc(b); 194 } else if ((*_status)[x] == UNMATCHED) { 195 augmentOnArc(b); 196 return; 197 } 198 } 199 } 200 } 201 } 202 } 203 204 void processSparse(const Node& n) { 205 _process = _last = 0; 206 _node_queue[_last++] = n; 207 while (_process != _last) { 208 Node u = _node_queue[_process++]; 209 for (OutArcIt a(_graph, u); a != INVALID; ++a) { 210 Node v = _graph.target(a); 211 212 if ((*_status)[v] == EVEN) { 213 if (_blossom_set>find(u) != _blossom_set>find(v)) { 214 shrinkOnEdge(a); 215 } 216 } else if ((*_status)[v] == MATCHED) { 217 extendOnArc(a); 218 } else if ((*_status)[v] == UNMATCHED) { 219 augmentOnArc(a); 220 return; 221 } 222 } 223 } 224 } 225 226 void shrinkOnEdge(const Edge& e) { 227 Node nca = INVALID; 228 229 { 230 std::set<Node> left_set, right_set; 231 232 Node left = (*_blossom_rep)[_blossom_set>find(_graph.u(e))]; 233 left_set.insert(left); 234 235 Node right = (*_blossom_rep)[_blossom_set>find(_graph.v(e))]; 236 right_set.insert(right); 237 238 while (true) { 239 if ((*_matching)[left] == INVALID) break; 240 left = _graph.target((*_matching)[left]); 241 left = (*_blossom_rep)[_blossom_set> 242 find(_graph.target((*_ear)[left]))]; 243 if (right_set.find(left) != right_set.end()) { 244 nca = left; 245 break; 246 } 247 left_set.insert(left); 248 249 if ((*_matching)[right] == INVALID) break; 250 right = _graph.target((*_matching)[right]); 251 right = (*_blossom_rep)[_blossom_set> 252 find(_graph.target((*_ear)[right]))]; 253 if (left_set.find(right) != left_set.end()) { 254 nca = right; 255 break; 256 } 257 right_set.insert(right); 258 } 259 260 if (nca == INVALID) { 261 if ((*_matching)[left] == INVALID) { 262 nca = right; 263 while (left_set.find(nca) == left_set.end()) { 264 nca = _graph.target((*_matching)[nca]); 265 nca =(*_blossom_rep)[_blossom_set> 266 find(_graph.target((*_ear)[nca]))]; 267 } 268 } else { 269 nca = left; 270 while (right_set.find(nca) == right_set.end()) { 271 nca = _graph.target((*_matching)[nca]); 272 nca = (*_blossom_rep)[_blossom_set> 273 find(_graph.target((*_ear)[nca]))]; 274 } 275 } 276 } 277 } 278 279 { 280 281 Node node = _graph.u(e); 282 Arc arc = _graph.direct(e, true); 283 Node base = (*_blossom_rep)[_blossom_set>find(node)]; 284 285 while (base != nca) { 286 _ear>set(node, arc); 287 288 Node n = node; 289 while (n != base) { 290 n = _graph.target((*_matching)[n]); 291 Arc a = (*_ear)[n]; 292 n = _graph.target(a); 293 _ear>set(n, _graph.oppositeArc(a)); 294 } 295 node = _graph.target((*_matching)[base]); 296 _tree_set>erase(base); 297 _tree_set>erase(node); 298 _blossom_set>insert(node, _blossom_set>find(base)); 299 _status>set(node, EVEN); 300 _node_queue[_last++] = node; 301 arc = _graph.oppositeArc((*_ear)[node]); 302 node = _graph.target((*_ear)[node]); 303 base = (*_blossom_rep)[_blossom_set>find(node)]; 304 _blossom_set>join(_graph.target(arc), base); 305 } 306 } 307 308 _blossom_rep>set(_blossom_set>find(nca), nca); 309 310 { 311 312 Node node = _graph.v(e); 313 Arc arc = _graph.direct(e, false); 314 Node base = (*_blossom_rep)[_blossom_set>find(node)]; 315 316 while (base != nca) { 317 _ear>set(node, arc); 318 319 Node n = node; 320 while (n != base) { 321 n = _graph.target((*_matching)[n]); 322 Arc a = (*_ear)[n]; 323 n = _graph.target(a); 324 _ear>set(n, _graph.oppositeArc(a)); 325 } 326 node = _graph.target((*_matching)[base]); 327 _tree_set>erase(base); 328 _tree_set>erase(node); 329 _blossom_set>insert(node, _blossom_set>find(base)); 330 _status>set(node, EVEN); 331 _node_queue[_last++] = node; 332 arc = _graph.oppositeArc((*_ear)[node]); 333 node = _graph.target((*_ear)[node]); 334 base = (*_blossom_rep)[_blossom_set>find(node)]; 335 _blossom_set>join(_graph.target(arc), base); 336 } 337 } 338 339 _blossom_rep>set(_blossom_set>find(nca), nca); 340 } 341 342 343 344 void extendOnArc(const Arc& a) { 345 Node base = _graph.source(a); 346 Node odd = _graph.target(a); 347 348 _ear>set(odd, _graph.oppositeArc(a)); 349 Node even = _graph.target((*_matching)[odd]); 350 _blossom_rep>set(_blossom_set>insert(even), even); 351 _status>set(odd, ODD); 352 _status>set(even, EVEN); 353 int tree = _tree_set>find((*_blossom_rep)[_blossom_set>find(base)]); 354 _tree_set>insert(odd, tree); 355 _tree_set>insert(even, tree); 356 _node_queue[_last++] = even; 357 358 } 359 360 void augmentOnArc(const Arc& a) { 361 Node even = _graph.source(a); 362 Node odd = _graph.target(a); 363 364 int tree = _tree_set>find((*_blossom_rep)[_blossom_set>find(even)]); 365 366 _matching>set(odd, _graph.oppositeArc(a)); 367 _status>set(odd, MATCHED); 368 369 Arc arc = (*_matching)[even]; 370 _matching>set(even, a); 371 372 while (arc != INVALID) { 373 odd = _graph.target(arc); 374 arc = (*_ear)[odd]; 375 even = _graph.target(arc); 376 _matching>set(odd, arc); 377 arc = (*_matching)[even]; 378 _matching>set(even, _graph.oppositeArc((*_matching)[odd])); 379 } 380 381 for (typename TreeSet::ItemIt it(*_tree_set, tree); 382 it != INVALID; ++it) { 383 if ((*_status)[it] == ODD) { 384 _status>set(it, MATCHED); 385 } else { 386 int blossom = _blossom_set>find(it); 387 for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom); 388 jt != INVALID; ++jt) { 389 _status>set(jt, MATCHED); 390 } 391 _blossom_set>eraseClass(blossom); 392 } 393 } 394 _tree_set>eraseClass(tree); 395 396 } 68 397 69 398 public: 70 399 71 ///\brief Indicates the GallaiEdmonds decomposition of the digraph. 72 /// 73 ///Indicates the GallaiEdmonds decomposition of the digraph, which 74 ///shows an upper bound on the size of a maximum matching. The 75 ///nodes with DecompType \c D induce a digraph with factorcritical 76 ///components, the nodes in \c A form the canonical barrier, and the 77 ///nodes in \c C induce a digraph having a perfect matching. 78 enum DecompType { 79 D=0, 80 A=1, 81 C=2 82 }; 83 84 protected: 85 86 static const int HEUR_density=2; 87 const Graph& g; 88 typename Graph::template NodeMap<Node> _mate; 89 typename Graph::template NodeMap<DecompType> position; 90 91 public: 92 93 MaxMatching(const Graph& _g) 94 : g(_g), _mate(_g), position(_g) {} 95 96 ///\brief Sets the actual matching to the empty matching. 97 /// 98 ///Sets the actual matching to the empty matching. 400 /// \brief Constructor 401 /// 402 /// Constructor. 403 MaxMatching(const Graph& graph) 404 : _graph(graph), _matching(0), _status(0), _ear(0), 405 _blossom_set_index(0), _blossom_set(0), _blossom_rep(0), 406 _tree_set_index(0), _tree_set(0) {} 407 408 ~MaxMatching() { 409 destroyStructures(); 410 } 411 412 /// \name Execution control 413 /// The simplest way to execute the algorithm is to use the member 414 /// \c run() member function. 415 /// \n 416 417 /// If you need more control on the execution, first you must call 418 /// \ref init(), \ref greedyInit() or \ref matchingInit() 419 /// functions, then you can start the algorithm with the \ref 420 /// startParse() or startDense() functions. 421 422 ///@{ 423 424 /// \brief Sets the actual matching to the empty matching. 425 /// 426 /// Sets the actual matching to the empty matching. 99 427 /// 100 428 void init() { 101 for(NodeIt v(g); v!=INVALID; ++v) { 102 _mate.set(v,INVALID); 103 position.set(v,C); 429 createStructures(); 430 for(NodeIt n(_graph); n != INVALID; ++n) { 431 _matching>set(n, INVALID); 432 _status>set(n, UNMATCHED); 104 433 } 105 434 } … … 109 438 ///For initial matchig it finds a maximal greedy matching. 110 439 void greedyInit() { 111 for(NodeIt v(g); v!=INVALID; ++v) { 112 _mate.set(v,INVALID); 113 position.set(v,C); 114 } 115 for(NodeIt v(g); v!=INVALID; ++v) 116 if ( _mate[v]==INVALID ) { 117 for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) { 118 Node y=g.runningNode(e); 119 if ( _mate[y]==INVALID && y!=v ) { 120 _mate.set(v,y); 121 _mate.set(y,v); 440 createStructures(); 441 for (NodeIt n(_graph); n != INVALID; ++n) { 442 _matching>set(n, INVALID); 443 _status>set(n, UNMATCHED); 444 } 445 for (NodeIt n(_graph); n != INVALID; ++n) { 446 if ((*_matching)[n] == INVALID) { 447 for (OutArcIt a(_graph, n); a != INVALID ; ++a) { 448 Node v = _graph.target(a); 449 if ((*_matching)[v] == INVALID && v != n) { 450 _matching>set(n, a); 451 _status>set(n, MATCHED); 452 _matching>set(v, _graph.oppositeArc(a)); 453 _status>set(v, MATCHED); 122 454 break; 123 455 } 124 456 } 125 457 } 126 } 127 128 ///\brief Initialize the matching from each nodes' mate. 129 /// 130 ///Initialize the matching from a \c Node valued \c Node map. This 131 ///map must be \e symmetric, i.e. if \c map[u]==v then \c 132 ///map[v]==u must hold, and \c uv will be an arc of the initial 133 ///matching. 134 template <typename MateMap> 135 void mateMapInit(MateMap& map) { 136 for(NodeIt v(g); v!=INVALID; ++v) { 137 _mate.set(v,map[v]); 138 position.set(v,C); 139 } 140 } 141 142 ///\brief Initialize the matching from a node map with the 143 ///incident matching arcs. 144 /// 145 ///Initialize the matching from an \c Edge valued \c Node map. \c 146 ///map[v] must be an \c Edge incident to \c v. This map must have 147 ///the property that if \c g.oppositeNode(u,map[u])==v then \c \c 148 ///g.oppositeNode(v,map[v])==u holds, and now some arc joining \c 149 ///u to \c v will be an arc of the matching. 150 template<typename MatchingMap> 151 void matchingMapInit(MatchingMap& map) { 152 for(NodeIt v(g); v!=INVALID; ++v) { 153 position.set(v,C); 154 Edge e=map[v]; 155 if ( e!=INVALID ) 156 _mate.set(v,g.oppositeNode(v,e)); 157 else 158 _mate.set(v,INVALID); 159 } 160 } 161 162 ///\brief Initialize the matching from the map containing the 163 ///undirected matching arcs. 164 /// 165 ///Initialize the matching from a \c bool valued \c Edge map. This 166 ///map must have the property that there are no two incident arcs 167 ///\c e, \c f with \c map[e]==map[f]==true. The arcs \c e with \c 168 ///map[e]==true form the matching. 458 } 459 } 460 461 462 /// \brief Initialize the matching from the map containing a matching. 463 /// 464 /// Initialize the matching from a \c bool valued \c Edge map. This 465 /// map must have the property that there are no two incident edges 466 /// with true value, ie. it contains a matching. 467 /// \return %True if the map contains a matching. 169 468 template <typename MatchingMap> 170 void matchingInit(MatchingMap& map) { 171 for(NodeIt v(g); v!=INVALID; ++v) { 172 _mate.set(v,INVALID); 173 position.set(v,C); 174 } 175 for(EdgeIt e(g); e!=INVALID; ++e) { 176 if ( map[e] ) { 177 Node u=g.u(e); 178 Node v=g.v(e); 179 _mate.set(u,v); 180 _mate.set(v,u); 181 } 182 } 183 } 184 185 186 ///\brief Runs Edmonds' algorithm. 187 /// 188 ///Runs Edmonds' algorithm for sparse digraphs (number of arcs < 189 ///2*number of nodes), and a heuristical Edmonds' algorithm with a 190 ///heuristic of postponing shrinks for dense digraphs. 469 bool matchingInit(const MatchingMap& matching) { 470 createStructures(); 471 472 for (NodeIt n(_graph); n != INVALID; ++n) { 473 _matching>set(n, INVALID); 474 _status>set(n, UNMATCHED); 475 } 476 for(EdgeIt e(_graph); e!=INVALID; ++e) { 477 if (matching[e]) { 478 479 Node u = _graph.u(e); 480 if ((*_matching)[u] != INVALID) return false; 481 _matching>set(u, _graph.direct(e, true)); 482 _status>set(u, MATCHED); 483 484 Node v = _graph.v(e); 485 if ((*_matching)[v] != INVALID) return false; 486 _matching>set(v, _graph.direct(e, false)); 487 _status>set(v, MATCHED); 488 } 489 } 490 return true; 491 } 492 493 /// \brief Starts Edmonds' algorithm 494 /// 495 /// If runs the original Edmonds' algorithm. 496 void startSparse() { 497 for(NodeIt n(_graph); n != INVALID; ++n) { 498 if ((*_status)[n] == UNMATCHED) { 499 (*_blossom_rep)[_blossom_set>insert(n)] = n; 500 _tree_set>insert(n); 501 _status>set(n, EVEN); 502 processSparse(n); 503 } 504 } 505 } 506 507 /// \brief Starts Edmonds' algorithm. 508 /// 509 /// It runs Edmonds' algorithm with a heuristic of postponing 510 /// shrinks, giving a faster algorithm for dense graphs. 511 void startDense() { 512 for(NodeIt n(_graph); n != INVALID; ++n) { 513 if ((*_status)[n] == UNMATCHED) { 514 (*_blossom_rep)[_blossom_set>insert(n)] = n; 515 _tree_set>insert(n); 516 _status>set(n, EVEN); 517 processDense(n); 518 } 519 } 520 } 521 522 523 /// \brief Runs Edmonds' algorithm 524 /// 525 /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>) 526 /// or Edmonds' algorithm with a heuristic of 527 /// postponing shrinks for dense graphs. 191 528 void run() { 192 if (countEdges( g) < HEUR_density * countNodes(g)) {529 if (countEdges(_graph) < 2 * countNodes(_graph)) { 193 530 greedyInit(); 194 531 startSparse(); … … 199 536 } 200 537 201 202 ///\brief Starts Edmonds' algorithm. 203 /// 204 ///If runs the original Edmonds' algorithm. 205 void startSparse() { 206 207 typename Graph::template NodeMap<Node> ear(g,INVALID); 208 //undefined for the base nodes of the blossoms (i.e. for the 209 //representative elements of UFE blossom) and for the nodes in C 210 211 UFECrossRef blossom_base(g); 212 UFE blossom(blossom_base); 213 NV rep(countNodes(g)); 214 215 EFECrossRef tree_base(g); 216 EFE tree(tree_base); 217 218 //If these UFE's would be members of the class then also 219 //blossom_base and tree_base should be a member. 220 221 //We build only one tree and the other vertices uncovered by the 222 //matching belong to C. (They can be considered as singleton 223 //trees.) If this tree can be augmented or no more 224 //grow/augmentation/shrink is possible then we return to this 225 //"for" cycle. 226 for(NodeIt v(g); v!=INVALID; ++v) { 227 if (position[v]==C && _mate[v]==INVALID) { 228 rep[blossom.insert(v)] = v; 229 tree.insert(v); 230 position.set(v,D); 231 normShrink(v, ear, blossom, rep, tree); 232 } 233 } 234 } 235 236 ///\brief Starts Edmonds' algorithm. 237 /// 238 ///It runs Edmonds' algorithm with a heuristic of postponing 239 ///shrinks, giving a faster algorithm for dense digraphs. 240 void startDense() { 241 242 typename Graph::template NodeMap<Node> ear(g,INVALID); 243 //undefined for the base nodes of the blossoms (i.e. for the 244 //representative elements of UFE blossom) and for the nodes in C 245 246 UFECrossRef blossom_base(g); 247 UFE blossom(blossom_base); 248 NV rep(countNodes(g)); 249 250 EFECrossRef tree_base(g); 251 EFE tree(tree_base); 252 253 //If these UFE's would be members of the class then also 254 //blossom_base and tree_base should be a member. 255 256 //We build only one tree and the other vertices uncovered by the 257 //matching belong to C. (They can be considered as singleton 258 //trees.) If this tree can be augmented or no more 259 //grow/augmentation/shrink is possible then we return to this 260 //"for" cycle. 261 for(NodeIt v(g); v!=INVALID; ++v) { 262 if ( position[v]==C && _mate[v]==INVALID ) { 263 rep[blossom.insert(v)] = v; 264 tree.insert(v); 265 position.set(v,D); 266 lateShrink(v, ear, blossom, rep, tree); 267 } 268 } 269 } 270 271 538 /// @} 539 540 /// \name Primal solution 541 /// Functions for get the primal solution, ie. the matching. 542 543 /// @{ 272 544 273 545 ///\brief Returns the size of the actual matching stored. 274 546 /// 275 547 ///Returns the size of the actual matching stored. After \ref 276 ///run() it returns the size of a maximum matching in the digraph. 277 int size() const { 278 int s=0; 279 for(NodeIt v(g); v!=INVALID; ++v) { 280 if ( _mate[v]!=INVALID ) { 281 ++s; 282 } 283 } 284 return s/2; 285 } 286 548 ///run() it returns the size of the maximum matching in the graph. 549 int matchingSize() const { 550 int size = 0; 551 for (NodeIt n(_graph); n != INVALID; ++n) { 552 if ((*_matching)[n] != INVALID) { 553 ++size; 554 } 555 } 556 return size / 2; 557 } 558 559 /// \brief Returns true when the edge is in the matching. 560 /// 561 /// Returns true when the edge is in the matching. 562 bool matching(const Edge& edge) const { 563 return edge == (*_matching)[_graph.u(edge)]; 564 } 565 566 /// \brief Returns the matching edge incident to the given node. 567 /// 568 /// Returns the matching edge of a \c node in the actual matching or 569 /// INVALID if the \c node is not covered by the actual matching. 570 Arc matching(const Node& n) const { 571 return (*_matching)[n]; 572 } 287 573 288 574 ///\brief Returns the mate of a node in the actual matching. 289 575 /// 290 ///Returns the mate of a \c node in the actual matching. 291 ///Returns INVALID if the \c node is not covered by the actual matching. 292 Node mate(const Node& node) const { 293 return _mate[node]; 294 } 295 296 ///\brief Returns the matching arc incident to the given node. 297 /// 298 ///Returns the matching arc of a \c node in the actual matching. 299 ///Returns INVALID if the \c node is not covered by the actual matching. 300 Edge matchingArc(const Node& node) const { 301 if (_mate[node] == INVALID) return INVALID; 302 Node n = node < _mate[node] ? node : _mate[node]; 303 for (IncEdgeIt e(g, n); e != INVALID; ++e) { 304 if (g.oppositeNode(n, e) == _mate[n]) { 305 return e; 306 } 307 } 308 return INVALID; 309 } 576 ///Returns the mate of a \c node in the actual matching or 577 ///INVALID if the \c node is not covered by the actual matching. 578 Node mate(const Node& n) const { 579 return (*_matching)[n] != INVALID ? 580 _graph.target((*_matching)[n]) : INVALID; 581 } 582 583 /// @} 584 585 /// \name Dual solution 586 /// Functions for get the dual solution, ie. the decomposition. 587 588 /// @{ 310 589 311 590 /// \brief Returns the class of the node in the EdmondsGallai … … 314 593 /// Returns the class of the node in the EdmondsGallai 315 594 /// decomposition. 316 DecompType decomposition(const Node& n){317 return position[n] == A;595 Status decomposition(const Node& n) const { 596 return (*_status)[n]; 318 597 } 319 598 … … 321 600 /// 322 601 /// Returns true when the node is in the barrier. 323 bool barrier(const Node& n) { 324 return position[n] == A; 325 } 326 327 ///\brief Gives back the matching in a \c Node of mates. 328 /// 329 ///Writes the stored matching to a \c Node valued \c Node map. The 330 ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c 331 ///map[v]==u will hold, and now \c uv is an arc of the matching. 332 template <typename MateMap> 333 void mateMap(MateMap& map) const { 334 for(NodeIt v(g); v!=INVALID; ++v) { 335 map.set(v,_mate[v]); 336 } 337 } 338 339 ///\brief Gives back the matching in an \c Edge valued \c Node 340 ///map. 341 /// 342 ///Writes the stored matching to an \c Edge valued \c Node 343 ///map. \c map[v] will be an \c Edge incident to \c v. This 344 ///map will have the property that if \c g.oppositeNode(u,map[u]) 345 ///== v then \c map[u]==map[v] holds, and now this arc is an arc 346 ///of the matching. 347 template <typename MatchingMap> 348 void matchingMap(MatchingMap& map) const { 349 typename Graph::template NodeMap<bool> todo(g,true); 350 for(NodeIt v(g); v!=INVALID; ++v) { 351 if (_mate[v]!=INVALID && v < _mate[v]) { 352 Node u=_mate[v]; 353 for(IncEdgeIt e(g,v); e!=INVALID; ++e) { 354 if ( g.runningNode(e) == u ) { 355 map.set(u,e); 356 map.set(v,e); 357 todo.set(u,false); 358 todo.set(v,false); 359 break; 360 } 361 } 362 } 363 } 364 } 365 366 367 ///\brief Gives back the matching in a \c bool valued \c Edge 368 ///map. 369 /// 370 ///Writes the matching stored to a \c bool valued \c Arc 371 ///map. This map will have the property that there are no two 372 ///incident arcs \c e, \c f with \c map[e]==map[f]==true. The 373 ///arcs \c e with \c map[e]==true form the matching. 374 template<typename MatchingMap> 375 void matching(MatchingMap& map) const { 376 for(EdgeIt e(g); e!=INVALID; ++e) map.set(e,false); 377 378 typename Graph::template NodeMap<bool> todo(g,true); 379 for(NodeIt v(g); v!=INVALID; ++v) { 380 if ( todo[v] && _mate[v]!=INVALID ) { 381 Node u=_mate[v]; 382 for(IncEdgeIt e(g,v); e!=INVALID; ++e) { 383 if ( g.runningNode(e) == u ) { 384 map.set(e,true); 385 todo.set(u,false); 386 todo.set(v,false); 387 break; 388 } 389 } 390 } 391 } 392 } 393 394 395 ///\brief Returns the canonical decomposition of the digraph after running 396 ///the algorithm. 397 /// 398 ///After calling any run methods of the class, it writes the 399 ///GallaiEdmonds canonical decomposition of the digraph. \c map 400 ///must be a node map of \ref DecompType 's. 401 template <typename DecompositionMap> 402 void decomposition(DecompositionMap& map) const { 403 for(NodeIt v(g); v!=INVALID; ++v) map.set(v,position[v]); 404 } 405 406 ///\brief Returns a barrier on the nodes. 407 /// 408 ///After calling any run methods of the class, it writes a 409 ///canonical barrier on the nodes. The odd component number of the 410 ///remaining digraph minus the barrier size is a lower bound for the 411 ///uncovered nodes in the digraph. The \c map must be a node map of 412 ///bools. 413 template <typename BarrierMap> 414 void barrier(BarrierMap& barrier) { 415 for(NodeIt v(g); v!=INVALID; ++v) barrier.set(v,position[v] == A); 416 } 417 418 private: 419 420 421 void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear, 422 UFE& blossom, NV& rep, EFE& tree) { 423 //We have one tree which we grow, and also shrink but only if it 424 //cannot be postponed. If we augment then we return to the "for" 425 //cycle of runEdmonds(). 426 427 std::queue<Node> Q; //queue of the totally unscanned nodes 428 Q.push(v); 429 std::queue<Node> R; 430 //queue of the nodes which must be scanned for a possible shrink 431 432 while ( !Q.empty() ) { 433 Node x=Q.front(); 434 Q.pop(); 435 for( IncEdgeIt e(g,x); e!= INVALID; ++e ) { 436 Node y=g.runningNode(e); 437 //growOrAugment grows if y is covered by the matching and 438 //augments if not. In this latter case it returns 1. 439 if (position[y]==C && 440 growOrAugment(y, x, ear, blossom, rep, tree, Q)) return; 441 } 442 R.push(x); 443 } 444 445 while ( !R.empty() ) { 446 Node x=R.front(); 447 R.pop(); 448 449 for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) { 450 Node y=g.runningNode(e); 451 452 if ( position[y] == D && blossom.find(x) != blossom.find(y) ) 453 //Recall that we have only one tree. 454 shrink( x, y, ear, blossom, rep, tree, Q); 455 456 while ( !Q.empty() ) { 457 Node z=Q.front(); 458 Q.pop(); 459 for( IncEdgeIt f(g,z); f!= INVALID; ++f ) { 460 Node w=g.runningNode(f); 461 //growOrAugment grows if y is covered by the matching and 462 //augments if not. In this latter case it returns 1. 463 if (position[w]==C && 464 growOrAugment(w, z, ear, blossom, rep, tree, Q)) return; 465 } 466 R.push(z); 467 } 468 } //for e 469 } // while ( !R.empty() ) 470 } 471 472 void normShrink(Node v, typename Graph::template NodeMap<Node>& ear, 473 UFE& blossom, NV& rep, EFE& tree) { 474 //We have one tree, which we grow and shrink. If we augment then we 475 //return to the "for" cycle of runEdmonds(). 476 477 std::queue<Node> Q; //queue of the unscanned nodes 478 Q.push(v); 479 while ( !Q.empty() ) { 480 481 Node x=Q.front(); 482 Q.pop(); 483 484 for( IncEdgeIt e(g,x); e!=INVALID; ++e ) { 485 Node y=g.runningNode(e); 486 487 switch ( position[y] ) { 488 case D: //x and y must be in the same tree 489 if ( blossom.find(x) != blossom.find(y)) 490 //x and y are in the same tree 491 shrink(x, y, ear, blossom, rep, tree, Q); 492 break; 493 case C: 494 //growOrAugment grows if y is covered by the matching and 495 //augments if not. In this latter case it returns 1. 496 if (growOrAugment(y, x, ear, blossom, rep, tree, Q)) return; 497 break; 498 default: break; 499 } 500 } 501 } 502 } 503 504 void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear, 505 UFE& blossom, NV& rep, EFE& tree,std::queue<Node>& Q) { 506 //x and y are the two adjacent vertices in two blossoms. 507 508 typename Graph::template NodeMap<bool> path(g,false); 509 510 Node b=rep[blossom.find(x)]; 511 path.set(b,true); 512 b=_mate[b]; 513 while ( b!=INVALID ) { 514 b=rep[blossom.find(ear[b])]; 515 path.set(b,true); 516 b=_mate[b]; 517 } //we go until the root through bases of blossoms and odd vertices 518 519 Node top=y; 520 Node middle=rep[blossom.find(top)]; 521 Node bottom=x; 522 while ( !path[middle] ) 523 shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q); 524 //Until we arrive to a node on the path, we update blossom, tree 525 //and the positions of the odd nodes. 526 527 Node base=middle; 528 top=x; 529 middle=rep[blossom.find(top)]; 530 bottom=y; 531 Node blossom_base=rep[blossom.find(base)]; 532 while ( middle!=blossom_base ) 533 shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q); 534 //Until we arrive to a node on the path, we update blossom, tree 535 //and the positions of the odd nodes. 536 537 rep[blossom.find(base)] = base; 538 } 539 540 void shrinkStep(Node& top, Node& middle, Node& bottom, 541 typename Graph::template NodeMap<Node>& ear, 542 UFE& blossom, NV& rep, EFE& tree, std::queue<Node>& Q) { 543 //We traverse a blossom and update everything. 544 545 ear.set(top,bottom); 546 Node t=top; 547 while ( t!=middle ) { 548 Node u=_mate[t]; 549 t=ear[u]; 550 ear.set(t,u); 551 } 552 bottom=_mate[middle]; 553 position.set(bottom,D); 554 Q.push(bottom); 555 top=ear[bottom]; 556 Node oldmiddle=middle; 557 middle=rep[blossom.find(top)]; 558 tree.erase(bottom); 559 tree.erase(oldmiddle); 560 blossom.insert(bottom); 561 blossom.join(bottom, oldmiddle); 562 blossom.join(top, oldmiddle); 563 } 564 565 566 567 bool growOrAugment(Node& y, Node& x, typename Graph::template 568 NodeMap<Node>& ear, UFE& blossom, NV& rep, EFE& tree, 569 std::queue<Node>& Q) { 570 //x is in a blossom in the tree, y is outside. If y is covered by 571 //the matching we grow, otherwise we augment. In this case we 572 //return 1. 573 574 if ( _mate[y]!=INVALID ) { //grow 575 ear.set(y,x); 576 Node w=_mate[y]; 577 rep[blossom.insert(w)] = w; 578 position.set(y,A); 579 position.set(w,D); 580 int t = tree.find(rep[blossom.find(x)]); 581 tree.insert(y,t); 582 tree.insert(w,t); 583 Q.push(w); 584 } else { //augment 585 augment(x, ear, blossom, rep, tree); 586 _mate.set(x,y); 587 _mate.set(y,x); 588 return true; 589 } 590 return false; 591 } 592 593 void augment(Node x, typename Graph::template NodeMap<Node>& ear, 594 UFE& blossom, NV& rep, EFE& tree) { 595 Node v=_mate[x]; 596 while ( v!=INVALID ) { 597 598 Node u=ear[v]; 599 _mate.set(v,u); 600 Node tmp=v; 601 v=_mate[u]; 602 _mate.set(u,tmp); 603 } 604 int y = tree.find(rep[blossom.find(x)]); 605 for (typename EFE::ItemIt tit(tree, y); tit != INVALID; ++tit) { 606 if ( position[tit] == D ) { 607 int b = blossom.find(tit); 608 for (typename UFE::ItemIt bit(blossom, b); bit != INVALID; ++bit) { 609 position.set(bit, C); 610 } 611 blossom.eraseClass(b); 612 } else position.set(tit, C); 613 } 614 tree.eraseClass(y); 615 616 } 602 bool barrier(const Node& n) const { 603 return (*_status)[n] == ODD; 604 } 605 606 /// @} 617 607 618 608 }; … … 628 618 /// 629 619 /// The maximum weighted matching problem is to find undirected 630 /// arcs in the digraph with maximum overall weight and no two of631 /// them shares their end points. The problem can be formulated with632 /// the next linear program:620 /// edges in the graph with maximum overall weight and no two of 621 /// them shares their ends. The problem can be formulated with the 622 /// following linear program. 633 623 /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f] 634 ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert  1}{2} \quad \forall B\in\mathcal{O}\f] 624 /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert  1}{2} 625 \quad \forall B\in\mathcal{O}\f] */ 635 626 /// \f[x_e \ge 0\quad \forall e\in E\f] 636 627 /// \f[\max \sum_{e\in E}x_ew_e\f] 637 /// where \f$\delta(X)\f$ is the set of arcs incident to a node in638 /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in639 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of640 /// the nodes.628 /// where \f$\delta(X)\f$ is the set of edges incident to a node in 629 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in 630 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality 631 /// subsets of the nodes. 641 632 /// 642 633 /// The algorithm calculates an optimal matching and a proof of the 643 634 /// optimality. The solution of the dual problem can be used to check 644 /// the result of the algorithm. The dual linear problem is the next: 645 /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f] 635 /// the result of the algorithm. The dual linear problem is the 636 /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)} 637 z_B \ge w_{uv} \quad \forall uv\in E\f] */ 646 638 /// \f[y_u \ge 0 \quad \forall u \in V\f] 647 639 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] 648 /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert  1}{2}z_B\f] 640 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} 641 \frac{\vert B \vert  1}{2}z_B\f] */ 649 642 /// 650 643 /// The algorithm can be executed with \c run() or the \c init() and … … 706 699 int _blossom_num; 707 700 708 typedef typename Graph::template NodeMap<int> NodeIntMap;709 typedef typename Graph::template ArcMap<int> ArcIntMap;710 typedef typename Graph::template EdgeMap<int> EdgeIntMap;711 701 typedef RangeMap<int> IntIntMap; 712 702 … … 715 705 }; 716 706 717 typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;707 typedef HeapUnionFind<Value, IntNodeMap> BlossomSet; 718 708 struct BlossomData { 719 709 int tree; … … 724 714 }; 725 715 726 NodeIntMap *_blossom_index;716 IntNodeMap *_blossom_index; 727 717 BlossomSet *_blossom_set; 728 718 RangeMap<BlossomData>* _blossom_data; 729 719 730 NodeIntMap *_node_index;731 ArcIntMap *_node_heap_index;720 IntNodeMap *_node_index; 721 IntArcMap *_node_heap_index; 732 722 733 723 struct NodeData { 734 724 735 NodeData( ArcIntMap& node_heap_index)725 NodeData(IntArcMap& node_heap_index) 736 726 : heap(node_heap_index) {} 737 727 738 728 int blossom; 739 729 Value pot; 740 BinHeap<Value, ArcIntMap> heap;730 BinHeap<Value, IntArcMap> heap; 741 731 std::map<int, Arc> heap_index; 742 732 … … 751 741 TreeSet *_tree_set; 752 742 753 NodeIntMap *_delta1_index;754 BinHeap<Value, NodeIntMap> *_delta1;743 IntNodeMap *_delta1_index; 744 BinHeap<Value, IntNodeMap> *_delta1; 755 745 756 746 IntIntMap *_delta2_index; 757 747 BinHeap<Value, IntIntMap> *_delta2; 758 748 759 EdgeIntMap *_delta3_index;760 BinHeap<Value, EdgeIntMap> *_delta3;749 IntEdgeMap *_delta3_index; 750 BinHeap<Value, IntEdgeMap> *_delta3; 761 751 762 752 IntIntMap *_delta4_index; … … 776 766 } 777 767 if (!_blossom_set) { 778 _blossom_index = new NodeIntMap(_graph);768 _blossom_index = new IntNodeMap(_graph); 779 769 _blossom_set = new BlossomSet(*_blossom_index); 780 770 _blossom_data = new RangeMap<BlossomData>(_blossom_num); … … 782 772 783 773 if (!_node_index) { 784 _node_index = new NodeIntMap(_graph);785 _node_heap_index = new ArcIntMap(_graph);774 _node_index = new IntNodeMap(_graph); 775 _node_heap_index = new IntArcMap(_graph); 786 776 _node_data = new RangeMap<NodeData>(_node_num, 787 777 NodeData(*_node_heap_index)); … … 793 783 } 794 784 if (!_delta1) { 795 _delta1_index = new NodeIntMap(_graph);796 _delta1 = new BinHeap<Value, NodeIntMap>(*_delta1_index);785 _delta1_index = new IntNodeMap(_graph); 786 _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index); 797 787 } 798 788 if (!_delta2) { … … 801 791 } 802 792 if (!_delta3) { 803 _delta3_index = new EdgeIntMap(_graph);804 _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index);793 _delta3_index = new IntEdgeMap(_graph); 794 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index); 805 795 } 806 796 if (!_delta4) { … … 1267 1257 1268 1258 1269 void augmentOn Arc(const Edge& arc) {1270 1271 int left = _blossom_set>find(_graph.u( arc));1272 int right = _blossom_set>find(_graph.v( arc));1259 void augmentOnEdge(const Edge& edge) { 1260 1261 int left = _blossom_set>find(_graph.u(edge)); 1262 int right = _blossom_set>find(_graph.v(edge)); 1273 1263 1274 1264 if ((*_blossom_data)[left].status == EVEN) { … … 1290 1280 } 1291 1281 1292 (*_blossom_data)[left].next = _graph.direct( arc, true);1293 (*_blossom_data)[right].next = _graph.direct( arc, false);1282 (*_blossom_data)[left].next = _graph.direct(edge, true); 1283 (*_blossom_data)[right].next = _graph.direct(edge, false); 1294 1284 } 1295 1285 … … 1311 1301 } 1312 1302 1313 void shrinkOn Arc(const Edge& edge, int tree) {1303 void shrinkOnEdge(const Edge& edge, int tree) { 1314 1304 int nca = 1; 1315 1305 std::vector<int> left_path, right_path; … … 1653 1643 1654 1644 for (ArcIt e(_graph); e != INVALID; ++e) { 1655 _node_heap_index>set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP);1645 _node_heap_index>set(e, BinHeap<Value, IntArcMap>::PRE_HEAP); 1656 1646 } 1657 1647 for (NodeIt n(_graph); n != INVALID; ++n) { … … 1770 1760 1771 1761 if (left_tree == right_tree) { 1772 shrinkOn Arc(e, left_tree);1762 shrinkOnEdge(e, left_tree); 1773 1763 } else { 1774 augmentOn Arc(e);1764 augmentOnEdge(e); 1775 1765 unmatched = 2; 1776 1766 } … … 1819 1809 } 1820 1810 1821 /// \brief Returns true when the arc is in the matching. 1822 /// 1823 /// Returns true when the arc is in the matching. 1824 bool matching(const Edge& arc) const { 1825 return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true); 1811 /// \brief Returns the cardinality of the matching. 1812 /// 1813 /// Returns the cardinality of the matching. 1814 int matchingSize() const { 1815 int num = 0; 1816 for (NodeIt n(_graph); n != INVALID; ++n) { 1817 if ((*_matching)[n] != INVALID) { 1818 ++num; 1819 } 1820 } 1821 return num /= 2; 1822 } 1823 1824 /// \brief Returns true when the edge is in the matching. 1825 /// 1826 /// Returns true when the edge is in the matching. 1827 bool matching(const Edge& edge) const { 1828 return edge == (*_matching)[_graph.u(edge)]; 1826 1829 } 1827 1830 … … 1914 1917 } 1915 1918 1916 /// \brief Invalid constructor.1917 ///1918 /// Invalid constructor.1919 BlossomIt(Invalid) : _index(1) {}1920 1921 1919 /// \brief Conversion to node. 1922 1920 /// 1923 1921 /// Conversion to node. 1924 1922 operator Node() const { 1925 return _algorithm ? _algorithm>_blossom_node_list[_index] : INVALID;1923 return _algorithm>_blossom_node_list[_index]; 1926 1924 } 1927 1925 … … 1931 1929 BlossomIt& operator++() { 1932 1930 ++_index; 1933 if (_index == _last) {1934 _index = 1;1935 }1936 1931 return *this; 1937 1932 } 1938 1933 1939 bool operator==(const BlossomIt& it) const { 1940 return _index == it._index; 1941 } 1942 bool operator!=(const BlossomIt& it) const { 1943 return _index != it._index; 1944 } 1934 /// \brief Validity checking 1935 /// 1936 /// Checks whether the iterator is invalid. 1937 bool operator==(Invalid) const { return _index == _last; } 1938 1939 /// \brief Validity checking 1940 /// 1941 /// Checks whether the iterator is valid. 1942 bool operator!=(Invalid) const { return _index != _last; } 1945 1943 1946 1944 private: … … 1959 1957 /// 1960 1958 /// This class provides an efficient implementation of Edmond's 1961 /// maximum weighted perfec rmatching algorithm. The implementation1959 /// maximum weighted perfect matching algorithm. The implementation 1962 1960 /// is based on extensive use of priority queues and provides 1963 1961 /// \f$O(nm\log(n))\f$ time complexity. 1964 1962 /// 1965 1963 /// The maximum weighted matching problem is to find undirected 1966 /// arcs in the digraph with maximum overall weight and no two of1967 /// them shares their end points and covers all nodes. The problem1968 /// can be formulated with the next linear program:1964 /// edges in the graph with maximum overall weight and no two of 1965 /// them shares their ends and covers all nodes. The problem can be 1966 /// formulated with the following linear program. 1969 1967 /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f] 1970 ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert  1}{2} \quad \forall B\in\mathcal{O}\f] 1968 /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert  1}{2} 1969 \quad \forall B\in\mathcal{O}\f] */ 1971 1970 /// \f[x_e \ge 0\quad \forall e\in E\f] 1972 1971 /// \f[\max \sum_{e\in E}x_ew_e\f] 1973 /// where \f$\delta(X)\f$ is the set of arcs incident to a node in1974 /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in1975 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of1976 /// the nodes.1972 /// where \f$\delta(X)\f$ is the set of edges incident to a node in 1973 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in 1974 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality 1975 /// subsets of the nodes. 1977 1976 /// 1978 1977 /// The algorithm calculates an optimal matching and a proof of the 1979 1978 /// optimality. The solution of the dual problem can be used to check 1980 /// the result of the algorithm. The dual linear problem is the next: 1981 /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f] 1979 /// the result of the algorithm. The dual linear problem is the 1980 /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge 1981 w_{uv} \quad \forall uv\in E\f] */ 1982 1982 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] 1983 /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert  1}{2}z_B\f] 1983 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} 1984 \frac{\vert B \vert  1}{2}z_B\f] */ 1984 1985 /// 1985 1986 /// The algorithm can be executed with \c run() or the \c init() and … … 2041 2042 int _blossom_num; 2042 2043 2043 typedef typename Graph::template NodeMap<int> NodeIntMap;2044 typedef typename Graph::template ArcMap<int> ArcIntMap;2045 typedef typename Graph::template EdgeMap<int> EdgeIntMap;2046 2044 typedef RangeMap<int> IntIntMap; 2047 2045 … … 2050 2048 }; 2051 2049 2052 typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;2050 typedef HeapUnionFind<Value, IntNodeMap> BlossomSet; 2053 2051 struct BlossomData { 2054 2052 int tree; … … 2058 2056 }; 2059 2057 2060 NodeIntMap *_blossom_index;2058 IntNodeMap *_blossom_index; 2061 2059 BlossomSet *_blossom_set; 2062 2060 RangeMap<BlossomData>* _blossom_data; 2063 2061 2064 NodeIntMap *_node_index;2065 ArcIntMap *_node_heap_index;2062 IntNodeMap *_node_index; 2063 IntArcMap *_node_heap_index; 2066 2064 2067 2065 struct NodeData { 2068 2066 2069 NodeData( ArcIntMap& node_heap_index)2067 NodeData(IntArcMap& node_heap_index) 2070 2068 : heap(node_heap_index) {} 2071 2069 2072 2070 int blossom; 2073 2071 Value pot; 2074 BinHeap<Value, ArcIntMap> heap;2072 BinHeap<Value, IntArcMap> heap; 2075 2073 std::map<int, Arc> heap_index; 2076 2074 … … 2088 2086 BinHeap<Value, IntIntMap> *_delta2; 2089 2087 2090 EdgeIntMap *_delta3_index;2091 BinHeap<Value, EdgeIntMap> *_delta3;2088 IntEdgeMap *_delta3_index; 2089 BinHeap<Value, IntEdgeMap> *_delta3; 2092 2090 2093 2091 IntIntMap *_delta4_index; … … 2107 2105 } 2108 2106 if (!_blossom_set) { 2109 _blossom_index = new NodeIntMap(_graph);2107 _blossom_index = new IntNodeMap(_graph); 2110 2108 _blossom_set = new BlossomSet(*_blossom_index); 2111 2109 _blossom_data = new RangeMap<BlossomData>(_blossom_num); … … 2113 2111 2114 2112 if (!_node_index) { 2115 _node_index = new NodeIntMap(_graph);2116 _node_heap_index = new ArcIntMap(_graph);2113 _node_index = new IntNodeMap(_graph); 2114 _node_heap_index = new IntArcMap(_graph); 2117 2115 _node_data = new RangeMap<NodeData>(_node_num, 2118 2116 NodeData(*_node_heap_index)); 2119 2117 } 2120 2118 … … 2128 2126 } 2129 2127 if (!_delta3) { 2130 _delta3_index = new EdgeIntMap(_graph);2131 _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index);2128 _delta3_index = new IntEdgeMap(_graph); 2129 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index); 2132 2130 } 2133 2131 if (!_delta4) { … … 2462 2460 } 2463 2461 2464 void augmentOn Arc(const Edge& arc) {2465 2466 int left = _blossom_set>find(_graph.u( arc));2467 int right = _blossom_set>find(_graph.v( arc));2462 void augmentOnEdge(const Edge& edge) { 2463 2464 int left = _blossom_set>find(_graph.u(edge)); 2465 int right = _blossom_set>find(_graph.v(edge)); 2468 2466 2469 2467 int left_tree = _tree_set>find(left); … … 2475 2473 destroyTree(right_tree); 2476 2474 2477 (*_blossom_data)[left].next = _graph.direct( arc, true);2478 (*_blossom_data)[right].next = _graph.direct( arc, false);2475 (*_blossom_data)[left].next = _graph.direct(edge, true); 2476 (*_blossom_data)[right].next = _graph.direct(edge, false); 2479 2477 } 2480 2478 … … 2496 2494 } 2497 2495 2498 void shrinkOn Arc(const Edge& edge, int tree) {2496 void shrinkOnEdge(const Edge& edge, int tree) { 2499 2497 int nca = 1; 2500 2498 std::vector<int> left_path, right_path; … … 2832 2830 2833 2831 for (ArcIt e(_graph); e != INVALID; ++e) { 2834 _node_heap_index>set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP);2832 _node_heap_index>set(e, BinHeap<Value, IntArcMap>::PRE_HEAP); 2835 2833 } 2836 2834 for (EdgeIt e(_graph); e != INVALID; ++e) { … … 2925 2923 2926 2924 if (left_tree == right_tree) { 2927 shrinkOn Arc(e, left_tree);2925 shrinkOnEdge(e, left_tree); 2928 2926 } else { 2929 augmentOn Arc(e);2927 augmentOnEdge(e); 2930 2928 unmatched = 2; 2931 2929 } … … 2975 2973 } 2976 2974 2977 /// \brief Returns true when the arcis in the matching.2978 /// 2979 /// Returns true when the arcis in the matching.2980 bool matching(const Edge& arc) const {2981 return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true);2982 } 2983 2984 /// \brief Returns the incident matching arc.2985 /// 2986 /// Returns the incident matching arc from given node.2975 /// \brief Returns true when the edge is in the matching. 2976 /// 2977 /// Returns true when the edge is in the matching. 2978 bool matching(const Edge& edge) const { 2979 return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge; 2980 } 2981 2982 /// \brief Returns the incident matching edge. 2983 /// 2984 /// Returns the incident matching arc from given edge. 2987 2985 Arc matching(const Node& node) const { 2988 2986 return (*_matching)[node]; … … 3067 3065 } 3068 3066 3069 /// \brief Invalid constructor.3070 ///3071 /// Invalid constructor.3072 BlossomIt(Invalid) : _index(1) {}3073 3074 3067 /// \brief Conversion to node. 3075 3068 /// 3076 3069 /// Conversion to node. 3077 3070 operator Node() const { 3078 return _algorithm ? _algorithm>_blossom_node_list[_index] : INVALID;3071 return _algorithm>_blossom_node_list[_index]; 3079 3072 } 3080 3073 … … 3084 3077 BlossomIt& operator++() { 3085 3078 ++_index; 3086 if (_index == _last) {3087 _index = 1;3088 }3089 3079 return *this; 3090 3080 } 3091 3081 3092 bool operator==(const BlossomIt& it) const { 3093 return _index == it._index; 3094 } 3095 bool operator!=(const BlossomIt& it) const { 3096 return _index != it._index; 3097 } 3082 /// \brief Validity checking 3083 /// 3084 /// Checks whether the iterator is invalid. 3085 bool operator==(Invalid) const { return _index == _last; } 3086 3087 /// \brief Validity checking 3088 /// 3089 /// Checks whether the iterator is valid. 3090 bool operator!=(Invalid) const { return _index != _last; } 3098 3091 3099 3092 private:
Note: See TracChangeset
for help on using the changeset viewer.