|
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-2008 |
|
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_MAX_MATCHING_H |
|
20 #define LEMON_MAX_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 |
|
32 ///\ingroup matching |
|
33 ///\file |
|
34 ///\brief Maximum matching algorithms in graph. |
|
35 |
|
36 namespace lemon { |
|
37 |
|
38 ///\ingroup matching |
|
39 /// |
|
40 ///\brief Edmonds' alternating forest maximum matching algorithm. |
|
41 /// |
|
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 /// |
|
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 Gallai-Edmonds decomposition of the digraph. The nodes |
|
49 ///in \c D induce a digraph with factor-critical 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. |
|
53 /// |
|
54 ///\param Digraph The graph type the algorithm runs on. |
|
55 template <typename Graph> |
|
56 class MaxMatching { |
|
57 |
|
58 protected: |
|
59 |
|
60 TEMPLATE_GRAPH_TYPEDEFS(Graph); |
|
61 |
|
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; |
|
68 |
|
69 public: |
|
70 |
|
71 ///\brief Indicates the Gallai-Edmonds decomposition of the digraph. |
|
72 /// |
|
73 ///Indicates the Gallai-Edmonds 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 factor-critical |
|
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. |
|
99 /// |
|
100 void init() { |
|
101 for(NodeIt v(g); v!=INVALID; ++v) { |
|
102 _mate.set(v,INVALID); |
|
103 position.set(v,C); |
|
104 } |
|
105 } |
|
106 |
|
107 ///\brief Finds a greedy matching for initial matching. |
|
108 /// |
|
109 ///For initial matchig it finds a maximal greedy matching. |
|
110 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); |
|
122 break; |
|
123 } |
|
124 } |
|
125 } |
|
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. |
|
169 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. |
|
191 void run() { |
|
192 if (countEdges(g) < HEUR_density * countNodes(g)) { |
|
193 greedyInit(); |
|
194 startSparse(); |
|
195 } else { |
|
196 init(); |
|
197 startDense(); |
|
198 } |
|
199 } |
|
200 |
|
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 |
|
272 |
|
273 ///\brief Returns the size of the actual matching stored. |
|
274 /// |
|
275 ///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 |
|
287 |
|
288 ///\brief Returns the mate of a node in the actual matching. |
|
289 /// |
|
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 } |
|
310 |
|
311 /// \brief Returns the class of the node in the Edmonds-Gallai |
|
312 /// decomposition. |
|
313 /// |
|
314 /// Returns the class of the node in the Edmonds-Gallai |
|
315 /// decomposition. |
|
316 DecompType decomposition(const Node& n) { |
|
317 return position[n] == A; |
|
318 } |
|
319 |
|
320 /// \brief Returns true when the node is in the barrier. |
|
321 /// |
|
322 /// 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 ///Gallai-Edmonds 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 } |
|
617 |
|
618 }; |
|
619 |
|
620 /// \ingroup matching |
|
621 /// |
|
622 /// \brief Weighted matching in general graphs |
|
623 /// |
|
624 /// This class provides an efficient implementation of Edmond's |
|
625 /// maximum weighted matching algorithm. The implementation is based |
|
626 /// on extensive use of priority queues and provides |
|
627 /// \f$O(nm\log(n))\f$ time complexity. |
|
628 /// |
|
629 /// The maximum weighted matching problem is to find undirected |
|
630 /// arcs in the digraph with maximum overall weight and no two of |
|
631 /// them shares their endpoints. The problem can be formulated with |
|
632 /// the next linear program: |
|
633 /// \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] |
|
635 /// \f[x_e \ge 0\quad \forall e\in E\f] |
|
636 /// \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 in |
|
638 /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in |
|
639 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of |
|
640 /// the nodes. |
|
641 /// |
|
642 /// The algorithm calculates an optimal matching and a proof of the |
|
643 /// 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] |
|
646 /// \f[y_u \ge 0 \quad \forall u \in V\f] |
|
647 /// \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] |
|
649 /// |
|
650 /// The algorithm can be executed with \c run() or the \c init() and |
|
651 /// then the \c start() member functions. After it the matching can |
|
652 /// be asked with \c matching() or mate() functions. The dual |
|
653 /// solution can be get with \c nodeValue(), \c blossomNum() and \c |
|
654 /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt |
|
655 /// "BlossomIt" nested class which is able to iterate on the nodes |
|
656 /// of a blossom. If the value type is integral then the dual |
|
657 /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4". |
|
658 template <typename _Graph, |
|
659 typename _WeightMap = typename _Graph::template EdgeMap<int> > |
|
660 class MaxWeightedMatching { |
|
661 public: |
|
662 |
|
663 typedef _Graph Graph; |
|
664 typedef _WeightMap WeightMap; |
|
665 typedef typename WeightMap::Value Value; |
|
666 |
|
667 /// \brief Scaling factor for dual solution |
|
668 /// |
|
669 /// Scaling factor for dual solution, it is equal to 4 or 1 |
|
670 /// according to the value type. |
|
671 static const int dualScale = |
|
672 std::numeric_limits<Value>::is_integer ? 4 : 1; |
|
673 |
|
674 typedef typename Graph::template NodeMap<typename Graph::Arc> |
|
675 MatchingMap; |
|
676 |
|
677 private: |
|
678 |
|
679 TEMPLATE_GRAPH_TYPEDEFS(Graph); |
|
680 |
|
681 typedef typename Graph::template NodeMap<Value> NodePotential; |
|
682 typedef std::vector<Node> BlossomNodeList; |
|
683 |
|
684 struct BlossomVariable { |
|
685 int begin, end; |
|
686 Value value; |
|
687 |
|
688 BlossomVariable(int _begin, int _end, Value _value) |
|
689 : begin(_begin), end(_end), value(_value) {} |
|
690 |
|
691 }; |
|
692 |
|
693 typedef std::vector<BlossomVariable> BlossomPotential; |
|
694 |
|
695 const Graph& _graph; |
|
696 const WeightMap& _weight; |
|
697 |
|
698 MatchingMap* _matching; |
|
699 |
|
700 NodePotential* _node_potential; |
|
701 |
|
702 BlossomPotential _blossom_potential; |
|
703 BlossomNodeList _blossom_node_list; |
|
704 |
|
705 int _node_num; |
|
706 int _blossom_num; |
|
707 |
|
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 typedef RangeMap<int> IntIntMap; |
|
712 |
|
713 enum Status { |
|
714 EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2 |
|
715 }; |
|
716 |
|
717 typedef HeapUnionFind<Value, NodeIntMap> BlossomSet; |
|
718 struct BlossomData { |
|
719 int tree; |
|
720 Status status; |
|
721 Arc pred, next; |
|
722 Value pot, offset; |
|
723 Node base; |
|
724 }; |
|
725 |
|
726 NodeIntMap *_blossom_index; |
|
727 BlossomSet *_blossom_set; |
|
728 RangeMap<BlossomData>* _blossom_data; |
|
729 |
|
730 NodeIntMap *_node_index; |
|
731 ArcIntMap *_node_heap_index; |
|
732 |
|
733 struct NodeData { |
|
734 |
|
735 NodeData(ArcIntMap& node_heap_index) |
|
736 : heap(node_heap_index) {} |
|
737 |
|
738 int blossom; |
|
739 Value pot; |
|
740 BinHeap<Value, ArcIntMap> heap; |
|
741 std::map<int, Arc> heap_index; |
|
742 |
|
743 int tree; |
|
744 }; |
|
745 |
|
746 RangeMap<NodeData>* _node_data; |
|
747 |
|
748 typedef ExtendFindEnum<IntIntMap> TreeSet; |
|
749 |
|
750 IntIntMap *_tree_set_index; |
|
751 TreeSet *_tree_set; |
|
752 |
|
753 NodeIntMap *_delta1_index; |
|
754 BinHeap<Value, NodeIntMap> *_delta1; |
|
755 |
|
756 IntIntMap *_delta2_index; |
|
757 BinHeap<Value, IntIntMap> *_delta2; |
|
758 |
|
759 EdgeIntMap *_delta3_index; |
|
760 BinHeap<Value, EdgeIntMap> *_delta3; |
|
761 |
|
762 IntIntMap *_delta4_index; |
|
763 BinHeap<Value, IntIntMap> *_delta4; |
|
764 |
|
765 Value _delta_sum; |
|
766 |
|
767 void createStructures() { |
|
768 _node_num = countNodes(_graph); |
|
769 _blossom_num = _node_num * 3 / 2; |
|
770 |
|
771 if (!_matching) { |
|
772 _matching = new MatchingMap(_graph); |
|
773 } |
|
774 if (!_node_potential) { |
|
775 _node_potential = new NodePotential(_graph); |
|
776 } |
|
777 if (!_blossom_set) { |
|
778 _blossom_index = new NodeIntMap(_graph); |
|
779 _blossom_set = new BlossomSet(*_blossom_index); |
|
780 _blossom_data = new RangeMap<BlossomData>(_blossom_num); |
|
781 } |
|
782 |
|
783 if (!_node_index) { |
|
784 _node_index = new NodeIntMap(_graph); |
|
785 _node_heap_index = new ArcIntMap(_graph); |
|
786 _node_data = new RangeMap<NodeData>(_node_num, |
|
787 NodeData(*_node_heap_index)); |
|
788 } |
|
789 |
|
790 if (!_tree_set) { |
|
791 _tree_set_index = new IntIntMap(_blossom_num); |
|
792 _tree_set = new TreeSet(*_tree_set_index); |
|
793 } |
|
794 if (!_delta1) { |
|
795 _delta1_index = new NodeIntMap(_graph); |
|
796 _delta1 = new BinHeap<Value, NodeIntMap>(*_delta1_index); |
|
797 } |
|
798 if (!_delta2) { |
|
799 _delta2_index = new IntIntMap(_blossom_num); |
|
800 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index); |
|
801 } |
|
802 if (!_delta3) { |
|
803 _delta3_index = new EdgeIntMap(_graph); |
|
804 _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index); |
|
805 } |
|
806 if (!_delta4) { |
|
807 _delta4_index = new IntIntMap(_blossom_num); |
|
808 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index); |
|
809 } |
|
810 } |
|
811 |
|
812 void destroyStructures() { |
|
813 _node_num = countNodes(_graph); |
|
814 _blossom_num = _node_num * 3 / 2; |
|
815 |
|
816 if (_matching) { |
|
817 delete _matching; |
|
818 } |
|
819 if (_node_potential) { |
|
820 delete _node_potential; |
|
821 } |
|
822 if (_blossom_set) { |
|
823 delete _blossom_index; |
|
824 delete _blossom_set; |
|
825 delete _blossom_data; |
|
826 } |
|
827 |
|
828 if (_node_index) { |
|
829 delete _node_index; |
|
830 delete _node_heap_index; |
|
831 delete _node_data; |
|
832 } |
|
833 |
|
834 if (_tree_set) { |
|
835 delete _tree_set_index; |
|
836 delete _tree_set; |
|
837 } |
|
838 if (_delta1) { |
|
839 delete _delta1_index; |
|
840 delete _delta1; |
|
841 } |
|
842 if (_delta2) { |
|
843 delete _delta2_index; |
|
844 delete _delta2; |
|
845 } |
|
846 if (_delta3) { |
|
847 delete _delta3_index; |
|
848 delete _delta3; |
|
849 } |
|
850 if (_delta4) { |
|
851 delete _delta4_index; |
|
852 delete _delta4; |
|
853 } |
|
854 } |
|
855 |
|
856 void matchedToEven(int blossom, int tree) { |
|
857 if (_delta2->state(blossom) == _delta2->IN_HEAP) { |
|
858 _delta2->erase(blossom); |
|
859 } |
|
860 |
|
861 if (!_blossom_set->trivial(blossom)) { |
|
862 (*_blossom_data)[blossom].pot -= |
|
863 2 * (_delta_sum - (*_blossom_data)[blossom].offset); |
|
864 } |
|
865 |
|
866 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
|
867 n != INVALID; ++n) { |
|
868 |
|
869 _blossom_set->increase(n, std::numeric_limits<Value>::max()); |
|
870 int ni = (*_node_index)[n]; |
|
871 |
|
872 (*_node_data)[ni].heap.clear(); |
|
873 (*_node_data)[ni].heap_index.clear(); |
|
874 |
|
875 (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset; |
|
876 |
|
877 _delta1->push(n, (*_node_data)[ni].pot); |
|
878 |
|
879 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
880 Node v = _graph.source(e); |
|
881 int vb = _blossom_set->find(v); |
|
882 int vi = (*_node_index)[v]; |
|
883 |
|
884 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - |
|
885 dualScale * _weight[e]; |
|
886 |
|
887 if ((*_blossom_data)[vb].status == EVEN) { |
|
888 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { |
|
889 _delta3->push(e, rw / 2); |
|
890 } |
|
891 } else if ((*_blossom_data)[vb].status == UNMATCHED) { |
|
892 if (_delta3->state(e) != _delta3->IN_HEAP) { |
|
893 _delta3->push(e, rw); |
|
894 } |
|
895 } else { |
|
896 typename std::map<int, Arc>::iterator it = |
|
897 (*_node_data)[vi].heap_index.find(tree); |
|
898 |
|
899 if (it != (*_node_data)[vi].heap_index.end()) { |
|
900 if ((*_node_data)[vi].heap[it->second] > rw) { |
|
901 (*_node_data)[vi].heap.replace(it->second, e); |
|
902 (*_node_data)[vi].heap.decrease(e, rw); |
|
903 it->second = e; |
|
904 } |
|
905 } else { |
|
906 (*_node_data)[vi].heap.push(e, rw); |
|
907 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); |
|
908 } |
|
909 |
|
910 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { |
|
911 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); |
|
912 |
|
913 if ((*_blossom_data)[vb].status == MATCHED) { |
|
914 if (_delta2->state(vb) != _delta2->IN_HEAP) { |
|
915 _delta2->push(vb, _blossom_set->classPrio(vb) - |
|
916 (*_blossom_data)[vb].offset); |
|
917 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - |
|
918 (*_blossom_data)[vb].offset){ |
|
919 _delta2->decrease(vb, _blossom_set->classPrio(vb) - |
|
920 (*_blossom_data)[vb].offset); |
|
921 } |
|
922 } |
|
923 } |
|
924 } |
|
925 } |
|
926 } |
|
927 (*_blossom_data)[blossom].offset = 0; |
|
928 } |
|
929 |
|
930 void matchedToOdd(int blossom) { |
|
931 if (_delta2->state(blossom) == _delta2->IN_HEAP) { |
|
932 _delta2->erase(blossom); |
|
933 } |
|
934 (*_blossom_data)[blossom].offset += _delta_sum; |
|
935 if (!_blossom_set->trivial(blossom)) { |
|
936 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + |
|
937 (*_blossom_data)[blossom].offset); |
|
938 } |
|
939 } |
|
940 |
|
941 void evenToMatched(int blossom, int tree) { |
|
942 if (!_blossom_set->trivial(blossom)) { |
|
943 (*_blossom_data)[blossom].pot += 2 * _delta_sum; |
|
944 } |
|
945 |
|
946 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
|
947 n != INVALID; ++n) { |
|
948 int ni = (*_node_index)[n]; |
|
949 (*_node_data)[ni].pot -= _delta_sum; |
|
950 |
|
951 _delta1->erase(n); |
|
952 |
|
953 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
954 Node v = _graph.source(e); |
|
955 int vb = _blossom_set->find(v); |
|
956 int vi = (*_node_index)[v]; |
|
957 |
|
958 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - |
|
959 dualScale * _weight[e]; |
|
960 |
|
961 if (vb == blossom) { |
|
962 if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
963 _delta3->erase(e); |
|
964 } |
|
965 } else if ((*_blossom_data)[vb].status == EVEN) { |
|
966 |
|
967 if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
968 _delta3->erase(e); |
|
969 } |
|
970 |
|
971 int vt = _tree_set->find(vb); |
|
972 |
|
973 if (vt != tree) { |
|
974 |
|
975 Arc r = _graph.oppositeArc(e); |
|
976 |
|
977 typename std::map<int, Arc>::iterator it = |
|
978 (*_node_data)[ni].heap_index.find(vt); |
|
979 |
|
980 if (it != (*_node_data)[ni].heap_index.end()) { |
|
981 if ((*_node_data)[ni].heap[it->second] > rw) { |
|
982 (*_node_data)[ni].heap.replace(it->second, r); |
|
983 (*_node_data)[ni].heap.decrease(r, rw); |
|
984 it->second = r; |
|
985 } |
|
986 } else { |
|
987 (*_node_data)[ni].heap.push(r, rw); |
|
988 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); |
|
989 } |
|
990 |
|
991 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { |
|
992 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); |
|
993 |
|
994 if (_delta2->state(blossom) != _delta2->IN_HEAP) { |
|
995 _delta2->push(blossom, _blossom_set->classPrio(blossom) - |
|
996 (*_blossom_data)[blossom].offset); |
|
997 } else if ((*_delta2)[blossom] > |
|
998 _blossom_set->classPrio(blossom) - |
|
999 (*_blossom_data)[blossom].offset){ |
|
1000 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - |
|
1001 (*_blossom_data)[blossom].offset); |
|
1002 } |
|
1003 } |
|
1004 } |
|
1005 |
|
1006 } else if ((*_blossom_data)[vb].status == UNMATCHED) { |
|
1007 if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
1008 _delta3->erase(e); |
|
1009 } |
|
1010 } else { |
|
1011 |
|
1012 typename std::map<int, Arc>::iterator it = |
|
1013 (*_node_data)[vi].heap_index.find(tree); |
|
1014 |
|
1015 if (it != (*_node_data)[vi].heap_index.end()) { |
|
1016 (*_node_data)[vi].heap.erase(it->second); |
|
1017 (*_node_data)[vi].heap_index.erase(it); |
|
1018 if ((*_node_data)[vi].heap.empty()) { |
|
1019 _blossom_set->increase(v, std::numeric_limits<Value>::max()); |
|
1020 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { |
|
1021 _blossom_set->increase(v, (*_node_data)[vi].heap.prio()); |
|
1022 } |
|
1023 |
|
1024 if ((*_blossom_data)[vb].status == MATCHED) { |
|
1025 if (_blossom_set->classPrio(vb) == |
|
1026 std::numeric_limits<Value>::max()) { |
|
1027 _delta2->erase(vb); |
|
1028 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - |
|
1029 (*_blossom_data)[vb].offset) { |
|
1030 _delta2->increase(vb, _blossom_set->classPrio(vb) - |
|
1031 (*_blossom_data)[vb].offset); |
|
1032 } |
|
1033 } |
|
1034 } |
|
1035 } |
|
1036 } |
|
1037 } |
|
1038 } |
|
1039 |
|
1040 void oddToMatched(int blossom) { |
|
1041 (*_blossom_data)[blossom].offset -= _delta_sum; |
|
1042 |
|
1043 if (_blossom_set->classPrio(blossom) != |
|
1044 std::numeric_limits<Value>::max()) { |
|
1045 _delta2->push(blossom, _blossom_set->classPrio(blossom) - |
|
1046 (*_blossom_data)[blossom].offset); |
|
1047 } |
|
1048 |
|
1049 if (!_blossom_set->trivial(blossom)) { |
|
1050 _delta4->erase(blossom); |
|
1051 } |
|
1052 } |
|
1053 |
|
1054 void oddToEven(int blossom, int tree) { |
|
1055 if (!_blossom_set->trivial(blossom)) { |
|
1056 _delta4->erase(blossom); |
|
1057 (*_blossom_data)[blossom].pot -= |
|
1058 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); |
|
1059 } |
|
1060 |
|
1061 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
|
1062 n != INVALID; ++n) { |
|
1063 int ni = (*_node_index)[n]; |
|
1064 |
|
1065 _blossom_set->increase(n, std::numeric_limits<Value>::max()); |
|
1066 |
|
1067 (*_node_data)[ni].heap.clear(); |
|
1068 (*_node_data)[ni].heap_index.clear(); |
|
1069 (*_node_data)[ni].pot += |
|
1070 2 * _delta_sum - (*_blossom_data)[blossom].offset; |
|
1071 |
|
1072 _delta1->push(n, (*_node_data)[ni].pot); |
|
1073 |
|
1074 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
1075 Node v = _graph.source(e); |
|
1076 int vb = _blossom_set->find(v); |
|
1077 int vi = (*_node_index)[v]; |
|
1078 |
|
1079 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - |
|
1080 dualScale * _weight[e]; |
|
1081 |
|
1082 if ((*_blossom_data)[vb].status == EVEN) { |
|
1083 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { |
|
1084 _delta3->push(e, rw / 2); |
|
1085 } |
|
1086 } else if ((*_blossom_data)[vb].status == UNMATCHED) { |
|
1087 if (_delta3->state(e) != _delta3->IN_HEAP) { |
|
1088 _delta3->push(e, rw); |
|
1089 } |
|
1090 } else { |
|
1091 |
|
1092 typename std::map<int, Arc>::iterator it = |
|
1093 (*_node_data)[vi].heap_index.find(tree); |
|
1094 |
|
1095 if (it != (*_node_data)[vi].heap_index.end()) { |
|
1096 if ((*_node_data)[vi].heap[it->second] > rw) { |
|
1097 (*_node_data)[vi].heap.replace(it->second, e); |
|
1098 (*_node_data)[vi].heap.decrease(e, rw); |
|
1099 it->second = e; |
|
1100 } |
|
1101 } else { |
|
1102 (*_node_data)[vi].heap.push(e, rw); |
|
1103 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); |
|
1104 } |
|
1105 |
|
1106 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { |
|
1107 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); |
|
1108 |
|
1109 if ((*_blossom_data)[vb].status == MATCHED) { |
|
1110 if (_delta2->state(vb) != _delta2->IN_HEAP) { |
|
1111 _delta2->push(vb, _blossom_set->classPrio(vb) - |
|
1112 (*_blossom_data)[vb].offset); |
|
1113 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - |
|
1114 (*_blossom_data)[vb].offset) { |
|
1115 _delta2->decrease(vb, _blossom_set->classPrio(vb) - |
|
1116 (*_blossom_data)[vb].offset); |
|
1117 } |
|
1118 } |
|
1119 } |
|
1120 } |
|
1121 } |
|
1122 } |
|
1123 (*_blossom_data)[blossom].offset = 0; |
|
1124 } |
|
1125 |
|
1126 |
|
1127 void matchedToUnmatched(int blossom) { |
|
1128 if (_delta2->state(blossom) == _delta2->IN_HEAP) { |
|
1129 _delta2->erase(blossom); |
|
1130 } |
|
1131 |
|
1132 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
|
1133 n != INVALID; ++n) { |
|
1134 int ni = (*_node_index)[n]; |
|
1135 |
|
1136 _blossom_set->increase(n, std::numeric_limits<Value>::max()); |
|
1137 |
|
1138 (*_node_data)[ni].heap.clear(); |
|
1139 (*_node_data)[ni].heap_index.clear(); |
|
1140 |
|
1141 for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
|
1142 Node v = _graph.target(e); |
|
1143 int vb = _blossom_set->find(v); |
|
1144 int vi = (*_node_index)[v]; |
|
1145 |
|
1146 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - |
|
1147 dualScale * _weight[e]; |
|
1148 |
|
1149 if ((*_blossom_data)[vb].status == EVEN) { |
|
1150 if (_delta3->state(e) != _delta3->IN_HEAP) { |
|
1151 _delta3->push(e, rw); |
|
1152 } |
|
1153 } |
|
1154 } |
|
1155 } |
|
1156 } |
|
1157 |
|
1158 void unmatchedToMatched(int blossom) { |
|
1159 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
|
1160 n != INVALID; ++n) { |
|
1161 int ni = (*_node_index)[n]; |
|
1162 |
|
1163 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
1164 Node v = _graph.source(e); |
|
1165 int vb = _blossom_set->find(v); |
|
1166 int vi = (*_node_index)[v]; |
|
1167 |
|
1168 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - |
|
1169 dualScale * _weight[e]; |
|
1170 |
|
1171 if (vb == blossom) { |
|
1172 if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
1173 _delta3->erase(e); |
|
1174 } |
|
1175 } else if ((*_blossom_data)[vb].status == EVEN) { |
|
1176 |
|
1177 if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
1178 _delta3->erase(e); |
|
1179 } |
|
1180 |
|
1181 int vt = _tree_set->find(vb); |
|
1182 |
|
1183 Arc r = _graph.oppositeArc(e); |
|
1184 |
|
1185 typename std::map<int, Arc>::iterator it = |
|
1186 (*_node_data)[ni].heap_index.find(vt); |
|
1187 |
|
1188 if (it != (*_node_data)[ni].heap_index.end()) { |
|
1189 if ((*_node_data)[ni].heap[it->second] > rw) { |
|
1190 (*_node_data)[ni].heap.replace(it->second, r); |
|
1191 (*_node_data)[ni].heap.decrease(r, rw); |
|
1192 it->second = r; |
|
1193 } |
|
1194 } else { |
|
1195 (*_node_data)[ni].heap.push(r, rw); |
|
1196 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); |
|
1197 } |
|
1198 |
|
1199 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { |
|
1200 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); |
|
1201 |
|
1202 if (_delta2->state(blossom) != _delta2->IN_HEAP) { |
|
1203 _delta2->push(blossom, _blossom_set->classPrio(blossom) - |
|
1204 (*_blossom_data)[blossom].offset); |
|
1205 } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)- |
|
1206 (*_blossom_data)[blossom].offset){ |
|
1207 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - |
|
1208 (*_blossom_data)[blossom].offset); |
|
1209 } |
|
1210 } |
|
1211 |
|
1212 } else if ((*_blossom_data)[vb].status == UNMATCHED) { |
|
1213 if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
1214 _delta3->erase(e); |
|
1215 } |
|
1216 } |
|
1217 } |
|
1218 } |
|
1219 } |
|
1220 |
|
1221 void alternatePath(int even, int tree) { |
|
1222 int odd; |
|
1223 |
|
1224 evenToMatched(even, tree); |
|
1225 (*_blossom_data)[even].status = MATCHED; |
|
1226 |
|
1227 while ((*_blossom_data)[even].pred != INVALID) { |
|
1228 odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred)); |
|
1229 (*_blossom_data)[odd].status = MATCHED; |
|
1230 oddToMatched(odd); |
|
1231 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred; |
|
1232 |
|
1233 even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred)); |
|
1234 (*_blossom_data)[even].status = MATCHED; |
|
1235 evenToMatched(even, tree); |
|
1236 (*_blossom_data)[even].next = |
|
1237 _graph.oppositeArc((*_blossom_data)[odd].pred); |
|
1238 } |
|
1239 |
|
1240 } |
|
1241 |
|
1242 void destroyTree(int tree) { |
|
1243 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) { |
|
1244 if ((*_blossom_data)[b].status == EVEN) { |
|
1245 (*_blossom_data)[b].status = MATCHED; |
|
1246 evenToMatched(b, tree); |
|
1247 } else if ((*_blossom_data)[b].status == ODD) { |
|
1248 (*_blossom_data)[b].status = MATCHED; |
|
1249 oddToMatched(b); |
|
1250 } |
|
1251 } |
|
1252 _tree_set->eraseClass(tree); |
|
1253 } |
|
1254 |
|
1255 |
|
1256 void unmatchNode(const Node& node) { |
|
1257 int blossom = _blossom_set->find(node); |
|
1258 int tree = _tree_set->find(blossom); |
|
1259 |
|
1260 alternatePath(blossom, tree); |
|
1261 destroyTree(tree); |
|
1262 |
|
1263 (*_blossom_data)[blossom].status = UNMATCHED; |
|
1264 (*_blossom_data)[blossom].base = node; |
|
1265 matchedToUnmatched(blossom); |
|
1266 } |
|
1267 |
|
1268 |
|
1269 void augmentOnArc(const Edge& arc) { |
|
1270 |
|
1271 int left = _blossom_set->find(_graph.u(arc)); |
|
1272 int right = _blossom_set->find(_graph.v(arc)); |
|
1273 |
|
1274 if ((*_blossom_data)[left].status == EVEN) { |
|
1275 int left_tree = _tree_set->find(left); |
|
1276 alternatePath(left, left_tree); |
|
1277 destroyTree(left_tree); |
|
1278 } else { |
|
1279 (*_blossom_data)[left].status = MATCHED; |
|
1280 unmatchedToMatched(left); |
|
1281 } |
|
1282 |
|
1283 if ((*_blossom_data)[right].status == EVEN) { |
|
1284 int right_tree = _tree_set->find(right); |
|
1285 alternatePath(right, right_tree); |
|
1286 destroyTree(right_tree); |
|
1287 } else { |
|
1288 (*_blossom_data)[right].status = MATCHED; |
|
1289 unmatchedToMatched(right); |
|
1290 } |
|
1291 |
|
1292 (*_blossom_data)[left].next = _graph.direct(arc, true); |
|
1293 (*_blossom_data)[right].next = _graph.direct(arc, false); |
|
1294 } |
|
1295 |
|
1296 void extendOnArc(const Arc& arc) { |
|
1297 int base = _blossom_set->find(_graph.target(arc)); |
|
1298 int tree = _tree_set->find(base); |
|
1299 |
|
1300 int odd = _blossom_set->find(_graph.source(arc)); |
|
1301 _tree_set->insert(odd, tree); |
|
1302 (*_blossom_data)[odd].status = ODD; |
|
1303 matchedToOdd(odd); |
|
1304 (*_blossom_data)[odd].pred = arc; |
|
1305 |
|
1306 int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next)); |
|
1307 (*_blossom_data)[even].pred = (*_blossom_data)[even].next; |
|
1308 _tree_set->insert(even, tree); |
|
1309 (*_blossom_data)[even].status = EVEN; |
|
1310 matchedToEven(even, tree); |
|
1311 } |
|
1312 |
|
1313 void shrinkOnArc(const Edge& edge, int tree) { |
|
1314 int nca = -1; |
|
1315 std::vector<int> left_path, right_path; |
|
1316 |
|
1317 { |
|
1318 std::set<int> left_set, right_set; |
|
1319 int left = _blossom_set->find(_graph.u(edge)); |
|
1320 left_path.push_back(left); |
|
1321 left_set.insert(left); |
|
1322 |
|
1323 int right = _blossom_set->find(_graph.v(edge)); |
|
1324 right_path.push_back(right); |
|
1325 right_set.insert(right); |
|
1326 |
|
1327 while (true) { |
|
1328 |
|
1329 if ((*_blossom_data)[left].pred == INVALID) break; |
|
1330 |
|
1331 left = |
|
1332 _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); |
|
1333 left_path.push_back(left); |
|
1334 left = |
|
1335 _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); |
|
1336 left_path.push_back(left); |
|
1337 |
|
1338 left_set.insert(left); |
|
1339 |
|
1340 if (right_set.find(left) != right_set.end()) { |
|
1341 nca = left; |
|
1342 break; |
|
1343 } |
|
1344 |
|
1345 if ((*_blossom_data)[right].pred == INVALID) break; |
|
1346 |
|
1347 right = |
|
1348 _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); |
|
1349 right_path.push_back(right); |
|
1350 right = |
|
1351 _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); |
|
1352 right_path.push_back(right); |
|
1353 |
|
1354 right_set.insert(right); |
|
1355 |
|
1356 if (left_set.find(right) != left_set.end()) { |
|
1357 nca = right; |
|
1358 break; |
|
1359 } |
|
1360 |
|
1361 } |
|
1362 |
|
1363 if (nca == -1) { |
|
1364 if ((*_blossom_data)[left].pred == INVALID) { |
|
1365 nca = right; |
|
1366 while (left_set.find(nca) == left_set.end()) { |
|
1367 nca = |
|
1368 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); |
|
1369 right_path.push_back(nca); |
|
1370 nca = |
|
1371 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); |
|
1372 right_path.push_back(nca); |
|
1373 } |
|
1374 } else { |
|
1375 nca = left; |
|
1376 while (right_set.find(nca) == right_set.end()) { |
|
1377 nca = |
|
1378 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); |
|
1379 left_path.push_back(nca); |
|
1380 nca = |
|
1381 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); |
|
1382 left_path.push_back(nca); |
|
1383 } |
|
1384 } |
|
1385 } |
|
1386 } |
|
1387 |
|
1388 std::vector<int> subblossoms; |
|
1389 Arc prev; |
|
1390 |
|
1391 prev = _graph.direct(edge, true); |
|
1392 for (int i = 0; left_path[i] != nca; i += 2) { |
|
1393 subblossoms.push_back(left_path[i]); |
|
1394 (*_blossom_data)[left_path[i]].next = prev; |
|
1395 _tree_set->erase(left_path[i]); |
|
1396 |
|
1397 subblossoms.push_back(left_path[i + 1]); |
|
1398 (*_blossom_data)[left_path[i + 1]].status = EVEN; |
|
1399 oddToEven(left_path[i + 1], tree); |
|
1400 _tree_set->erase(left_path[i + 1]); |
|
1401 prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred); |
|
1402 } |
|
1403 |
|
1404 int k = 0; |
|
1405 while (right_path[k] != nca) ++k; |
|
1406 |
|
1407 subblossoms.push_back(nca); |
|
1408 (*_blossom_data)[nca].next = prev; |
|
1409 |
|
1410 for (int i = k - 2; i >= 0; i -= 2) { |
|
1411 subblossoms.push_back(right_path[i + 1]); |
|
1412 (*_blossom_data)[right_path[i + 1]].status = EVEN; |
|
1413 oddToEven(right_path[i + 1], tree); |
|
1414 _tree_set->erase(right_path[i + 1]); |
|
1415 |
|
1416 (*_blossom_data)[right_path[i + 1]].next = |
|
1417 (*_blossom_data)[right_path[i + 1]].pred; |
|
1418 |
|
1419 subblossoms.push_back(right_path[i]); |
|
1420 _tree_set->erase(right_path[i]); |
|
1421 } |
|
1422 |
|
1423 int surface = |
|
1424 _blossom_set->join(subblossoms.begin(), subblossoms.end()); |
|
1425 |
|
1426 for (int i = 0; i < int(subblossoms.size()); ++i) { |
|
1427 if (!_blossom_set->trivial(subblossoms[i])) { |
|
1428 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum; |
|
1429 } |
|
1430 (*_blossom_data)[subblossoms[i]].status = MATCHED; |
|
1431 } |
|
1432 |
|
1433 (*_blossom_data)[surface].pot = -2 * _delta_sum; |
|
1434 (*_blossom_data)[surface].offset = 0; |
|
1435 (*_blossom_data)[surface].status = EVEN; |
|
1436 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred; |
|
1437 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred; |
|
1438 |
|
1439 _tree_set->insert(surface, tree); |
|
1440 _tree_set->erase(nca); |
|
1441 } |
|
1442 |
|
1443 void splitBlossom(int blossom) { |
|
1444 Arc next = (*_blossom_data)[blossom].next; |
|
1445 Arc pred = (*_blossom_data)[blossom].pred; |
|
1446 |
|
1447 int tree = _tree_set->find(blossom); |
|
1448 |
|
1449 (*_blossom_data)[blossom].status = MATCHED; |
|
1450 oddToMatched(blossom); |
|
1451 if (_delta2->state(blossom) == _delta2->IN_HEAP) { |
|
1452 _delta2->erase(blossom); |
|
1453 } |
|
1454 |
|
1455 std::vector<int> subblossoms; |
|
1456 _blossom_set->split(blossom, std::back_inserter(subblossoms)); |
|
1457 |
|
1458 Value offset = (*_blossom_data)[blossom].offset; |
|
1459 int b = _blossom_set->find(_graph.source(pred)); |
|
1460 int d = _blossom_set->find(_graph.source(next)); |
|
1461 |
|
1462 int ib = -1, id = -1; |
|
1463 for (int i = 0; i < int(subblossoms.size()); ++i) { |
|
1464 if (subblossoms[i] == b) ib = i; |
|
1465 if (subblossoms[i] == d) id = i; |
|
1466 |
|
1467 (*_blossom_data)[subblossoms[i]].offset = offset; |
|
1468 if (!_blossom_set->trivial(subblossoms[i])) { |
|
1469 (*_blossom_data)[subblossoms[i]].pot -= 2 * offset; |
|
1470 } |
|
1471 if (_blossom_set->classPrio(subblossoms[i]) != |
|
1472 std::numeric_limits<Value>::max()) { |
|
1473 _delta2->push(subblossoms[i], |
|
1474 _blossom_set->classPrio(subblossoms[i]) - |
|
1475 (*_blossom_data)[subblossoms[i]].offset); |
|
1476 } |
|
1477 } |
|
1478 |
|
1479 if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) { |
|
1480 for (int i = (id + 1) % subblossoms.size(); |
|
1481 i != ib; i = (i + 2) % subblossoms.size()) { |
|
1482 int sb = subblossoms[i]; |
|
1483 int tb = subblossoms[(i + 1) % subblossoms.size()]; |
|
1484 (*_blossom_data)[sb].next = |
|
1485 _graph.oppositeArc((*_blossom_data)[tb].next); |
|
1486 } |
|
1487 |
|
1488 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) { |
|
1489 int sb = subblossoms[i]; |
|
1490 int tb = subblossoms[(i + 1) % subblossoms.size()]; |
|
1491 int ub = subblossoms[(i + 2) % subblossoms.size()]; |
|
1492 |
|
1493 (*_blossom_data)[sb].status = ODD; |
|
1494 matchedToOdd(sb); |
|
1495 _tree_set->insert(sb, tree); |
|
1496 (*_blossom_data)[sb].pred = pred; |
|
1497 (*_blossom_data)[sb].next = |
|
1498 _graph.oppositeArc((*_blossom_data)[tb].next); |
|
1499 |
|
1500 pred = (*_blossom_data)[ub].next; |
|
1501 |
|
1502 (*_blossom_data)[tb].status = EVEN; |
|
1503 matchedToEven(tb, tree); |
|
1504 _tree_set->insert(tb, tree); |
|
1505 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next; |
|
1506 } |
|
1507 |
|
1508 (*_blossom_data)[subblossoms[id]].status = ODD; |
|
1509 matchedToOdd(subblossoms[id]); |
|
1510 _tree_set->insert(subblossoms[id], tree); |
|
1511 (*_blossom_data)[subblossoms[id]].next = next; |
|
1512 (*_blossom_data)[subblossoms[id]].pred = pred; |
|
1513 |
|
1514 } else { |
|
1515 |
|
1516 for (int i = (ib + 1) % subblossoms.size(); |
|
1517 i != id; i = (i + 2) % subblossoms.size()) { |
|
1518 int sb = subblossoms[i]; |
|
1519 int tb = subblossoms[(i + 1) % subblossoms.size()]; |
|
1520 (*_blossom_data)[sb].next = |
|
1521 _graph.oppositeArc((*_blossom_data)[tb].next); |
|
1522 } |
|
1523 |
|
1524 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) { |
|
1525 int sb = subblossoms[i]; |
|
1526 int tb = subblossoms[(i + 1) % subblossoms.size()]; |
|
1527 int ub = subblossoms[(i + 2) % subblossoms.size()]; |
|
1528 |
|
1529 (*_blossom_data)[sb].status = ODD; |
|
1530 matchedToOdd(sb); |
|
1531 _tree_set->insert(sb, tree); |
|
1532 (*_blossom_data)[sb].next = next; |
|
1533 (*_blossom_data)[sb].pred = |
|
1534 _graph.oppositeArc((*_blossom_data)[tb].next); |
|
1535 |
|
1536 (*_blossom_data)[tb].status = EVEN; |
|
1537 matchedToEven(tb, tree); |
|
1538 _tree_set->insert(tb, tree); |
|
1539 (*_blossom_data)[tb].pred = |
|
1540 (*_blossom_data)[tb].next = |
|
1541 _graph.oppositeArc((*_blossom_data)[ub].next); |
|
1542 next = (*_blossom_data)[ub].next; |
|
1543 } |
|
1544 |
|
1545 (*_blossom_data)[subblossoms[ib]].status = ODD; |
|
1546 matchedToOdd(subblossoms[ib]); |
|
1547 _tree_set->insert(subblossoms[ib], tree); |
|
1548 (*_blossom_data)[subblossoms[ib]].next = next; |
|
1549 (*_blossom_data)[subblossoms[ib]].pred = pred; |
|
1550 } |
|
1551 _tree_set->erase(blossom); |
|
1552 } |
|
1553 |
|
1554 void extractBlossom(int blossom, const Node& base, const Arc& matching) { |
|
1555 if (_blossom_set->trivial(blossom)) { |
|
1556 int bi = (*_node_index)[base]; |
|
1557 Value pot = (*_node_data)[bi].pot; |
|
1558 |
|
1559 _matching->set(base, matching); |
|
1560 _blossom_node_list.push_back(base); |
|
1561 _node_potential->set(base, pot); |
|
1562 } else { |
|
1563 |
|
1564 Value pot = (*_blossom_data)[blossom].pot; |
|
1565 int bn = _blossom_node_list.size(); |
|
1566 |
|
1567 std::vector<int> subblossoms; |
|
1568 _blossom_set->split(blossom, std::back_inserter(subblossoms)); |
|
1569 int b = _blossom_set->find(base); |
|
1570 int ib = -1; |
|
1571 for (int i = 0; i < int(subblossoms.size()); ++i) { |
|
1572 if (subblossoms[i] == b) { ib = i; break; } |
|
1573 } |
|
1574 |
|
1575 for (int i = 1; i < int(subblossoms.size()); i += 2) { |
|
1576 int sb = subblossoms[(ib + i) % subblossoms.size()]; |
|
1577 int tb = subblossoms[(ib + i + 1) % subblossoms.size()]; |
|
1578 |
|
1579 Arc m = (*_blossom_data)[tb].next; |
|
1580 extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m)); |
|
1581 extractBlossom(tb, _graph.source(m), m); |
|
1582 } |
|
1583 extractBlossom(subblossoms[ib], base, matching); |
|
1584 |
|
1585 int en = _blossom_node_list.size(); |
|
1586 |
|
1587 _blossom_potential.push_back(BlossomVariable(bn, en, pot)); |
|
1588 } |
|
1589 } |
|
1590 |
|
1591 void extractMatching() { |
|
1592 std::vector<int> blossoms; |
|
1593 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) { |
|
1594 blossoms.push_back(c); |
|
1595 } |
|
1596 |
|
1597 for (int i = 0; i < int(blossoms.size()); ++i) { |
|
1598 if ((*_blossom_data)[blossoms[i]].status == MATCHED) { |
|
1599 |
|
1600 Value offset = (*_blossom_data)[blossoms[i]].offset; |
|
1601 (*_blossom_data)[blossoms[i]].pot += 2 * offset; |
|
1602 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); |
|
1603 n != INVALID; ++n) { |
|
1604 (*_node_data)[(*_node_index)[n]].pot -= offset; |
|
1605 } |
|
1606 |
|
1607 Arc matching = (*_blossom_data)[blossoms[i]].next; |
|
1608 Node base = _graph.source(matching); |
|
1609 extractBlossom(blossoms[i], base, matching); |
|
1610 } else { |
|
1611 Node base = (*_blossom_data)[blossoms[i]].base; |
|
1612 extractBlossom(blossoms[i], base, INVALID); |
|
1613 } |
|
1614 } |
|
1615 } |
|
1616 |
|
1617 public: |
|
1618 |
|
1619 /// \brief Constructor |
|
1620 /// |
|
1621 /// Constructor. |
|
1622 MaxWeightedMatching(const Graph& graph, const WeightMap& weight) |
|
1623 : _graph(graph), _weight(weight), _matching(0), |
|
1624 _node_potential(0), _blossom_potential(), _blossom_node_list(), |
|
1625 _node_num(0), _blossom_num(0), |
|
1626 |
|
1627 _blossom_index(0), _blossom_set(0), _blossom_data(0), |
|
1628 _node_index(0), _node_heap_index(0), _node_data(0), |
|
1629 _tree_set_index(0), _tree_set(0), |
|
1630 |
|
1631 _delta1_index(0), _delta1(0), |
|
1632 _delta2_index(0), _delta2(0), |
|
1633 _delta3_index(0), _delta3(0), |
|
1634 _delta4_index(0), _delta4(0), |
|
1635 |
|
1636 _delta_sum() {} |
|
1637 |
|
1638 ~MaxWeightedMatching() { |
|
1639 destroyStructures(); |
|
1640 } |
|
1641 |
|
1642 /// \name Execution control |
|
1643 /// The simplest way to execute the algorithm is to use the member |
|
1644 /// \c run() member function. |
|
1645 |
|
1646 ///@{ |
|
1647 |
|
1648 /// \brief Initialize the algorithm |
|
1649 /// |
|
1650 /// Initialize the algorithm |
|
1651 void init() { |
|
1652 createStructures(); |
|
1653 |
|
1654 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
1655 _node_heap_index->set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP); |
|
1656 } |
|
1657 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
1658 _delta1_index->set(n, _delta1->PRE_HEAP); |
|
1659 } |
|
1660 for (EdgeIt e(_graph); e != INVALID; ++e) { |
|
1661 _delta3_index->set(e, _delta3->PRE_HEAP); |
|
1662 } |
|
1663 for (int i = 0; i < _blossom_num; ++i) { |
|
1664 _delta2_index->set(i, _delta2->PRE_HEAP); |
|
1665 _delta4_index->set(i, _delta4->PRE_HEAP); |
|
1666 } |
|
1667 |
|
1668 int index = 0; |
|
1669 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
1670 Value max = 0; |
|
1671 for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
|
1672 if (_graph.target(e) == n) continue; |
|
1673 if ((dualScale * _weight[e]) / 2 > max) { |
|
1674 max = (dualScale * _weight[e]) / 2; |
|
1675 } |
|
1676 } |
|
1677 _node_index->set(n, index); |
|
1678 (*_node_data)[index].pot = max; |
|
1679 _delta1->push(n, max); |
|
1680 int blossom = |
|
1681 _blossom_set->insert(n, std::numeric_limits<Value>::max()); |
|
1682 |
|
1683 _tree_set->insert(blossom); |
|
1684 |
|
1685 (*_blossom_data)[blossom].status = EVEN; |
|
1686 (*_blossom_data)[blossom].pred = INVALID; |
|
1687 (*_blossom_data)[blossom].next = INVALID; |
|
1688 (*_blossom_data)[blossom].pot = 0; |
|
1689 (*_blossom_data)[blossom].offset = 0; |
|
1690 ++index; |
|
1691 } |
|
1692 for (EdgeIt e(_graph); e != INVALID; ++e) { |
|
1693 int si = (*_node_index)[_graph.u(e)]; |
|
1694 int ti = (*_node_index)[_graph.v(e)]; |
|
1695 if (_graph.u(e) != _graph.v(e)) { |
|
1696 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - |
|
1697 dualScale * _weight[e]) / 2); |
|
1698 } |
|
1699 } |
|
1700 } |
|
1701 |
|
1702 /// \brief Starts the algorithm |
|
1703 /// |
|
1704 /// Starts the algorithm |
|
1705 void start() { |
|
1706 enum OpType { |
|
1707 D1, D2, D3, D4 |
|
1708 }; |
|
1709 |
|
1710 int unmatched = _node_num; |
|
1711 while (unmatched > 0) { |
|
1712 Value d1 = !_delta1->empty() ? |
|
1713 _delta1->prio() : std::numeric_limits<Value>::max(); |
|
1714 |
|
1715 Value d2 = !_delta2->empty() ? |
|
1716 _delta2->prio() : std::numeric_limits<Value>::max(); |
|
1717 |
|
1718 Value d3 = !_delta3->empty() ? |
|
1719 _delta3->prio() : std::numeric_limits<Value>::max(); |
|
1720 |
|
1721 Value d4 = !_delta4->empty() ? |
|
1722 _delta4->prio() : std::numeric_limits<Value>::max(); |
|
1723 |
|
1724 _delta_sum = d1; OpType ot = D1; |
|
1725 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } |
|
1726 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } |
|
1727 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } |
|
1728 |
|
1729 |
|
1730 switch (ot) { |
|
1731 case D1: |
|
1732 { |
|
1733 Node n = _delta1->top(); |
|
1734 unmatchNode(n); |
|
1735 --unmatched; |
|
1736 } |
|
1737 break; |
|
1738 case D2: |
|
1739 { |
|
1740 int blossom = _delta2->top(); |
|
1741 Node n = _blossom_set->classTop(blossom); |
|
1742 Arc e = (*_node_data)[(*_node_index)[n]].heap.top(); |
|
1743 extendOnArc(e); |
|
1744 } |
|
1745 break; |
|
1746 case D3: |
|
1747 { |
|
1748 Edge e = _delta3->top(); |
|
1749 |
|
1750 int left_blossom = _blossom_set->find(_graph.u(e)); |
|
1751 int right_blossom = _blossom_set->find(_graph.v(e)); |
|
1752 |
|
1753 if (left_blossom == right_blossom) { |
|
1754 _delta3->pop(); |
|
1755 } else { |
|
1756 int left_tree; |
|
1757 if ((*_blossom_data)[left_blossom].status == EVEN) { |
|
1758 left_tree = _tree_set->find(left_blossom); |
|
1759 } else { |
|
1760 left_tree = -1; |
|
1761 ++unmatched; |
|
1762 } |
|
1763 int right_tree; |
|
1764 if ((*_blossom_data)[right_blossom].status == EVEN) { |
|
1765 right_tree = _tree_set->find(right_blossom); |
|
1766 } else { |
|
1767 right_tree = -1; |
|
1768 ++unmatched; |
|
1769 } |
|
1770 |
|
1771 if (left_tree == right_tree) { |
|
1772 shrinkOnArc(e, left_tree); |
|
1773 } else { |
|
1774 augmentOnArc(e); |
|
1775 unmatched -= 2; |
|
1776 } |
|
1777 } |
|
1778 } break; |
|
1779 case D4: |
|
1780 splitBlossom(_delta4->top()); |
|
1781 break; |
|
1782 } |
|
1783 } |
|
1784 extractMatching(); |
|
1785 } |
|
1786 |
|
1787 /// \brief Runs %MaxWeightedMatching algorithm. |
|
1788 /// |
|
1789 /// This method runs the %MaxWeightedMatching algorithm. |
|
1790 /// |
|
1791 /// \note mwm.run() is just a shortcut of the following code. |
|
1792 /// \code |
|
1793 /// mwm.init(); |
|
1794 /// mwm.start(); |
|
1795 /// \endcode |
|
1796 void run() { |
|
1797 init(); |
|
1798 start(); |
|
1799 } |
|
1800 |
|
1801 /// @} |
|
1802 |
|
1803 /// \name Primal solution |
|
1804 /// Functions for get the primal solution, ie. the matching. |
|
1805 |
|
1806 /// @{ |
|
1807 |
|
1808 /// \brief Returns the matching value. |
|
1809 /// |
|
1810 /// Returns the matching value. |
|
1811 Value matchingValue() const { |
|
1812 Value sum = 0; |
|
1813 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
1814 if ((*_matching)[n] != INVALID) { |
|
1815 sum += _weight[(*_matching)[n]]; |
|
1816 } |
|
1817 } |
|
1818 return sum /= 2; |
|
1819 } |
|
1820 |
|
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); |
|
1826 } |
|
1827 |
|
1828 /// \brief Returns the incident matching arc. |
|
1829 /// |
|
1830 /// Returns the incident matching arc from given node. If the |
|
1831 /// node is not matched then it gives back \c INVALID. |
|
1832 Arc matching(const Node& node) const { |
|
1833 return (*_matching)[node]; |
|
1834 } |
|
1835 |
|
1836 /// \brief Returns the mate of the node. |
|
1837 /// |
|
1838 /// Returns the adjancent node in a mathcing arc. If the node is |
|
1839 /// not matched then it gives back \c INVALID. |
|
1840 Node mate(const Node& node) const { |
|
1841 return (*_matching)[node] != INVALID ? |
|
1842 _graph.target((*_matching)[node]) : INVALID; |
|
1843 } |
|
1844 |
|
1845 /// @} |
|
1846 |
|
1847 /// \name Dual solution |
|
1848 /// Functions for get the dual solution. |
|
1849 |
|
1850 /// @{ |
|
1851 |
|
1852 /// \brief Returns the value of the dual solution. |
|
1853 /// |
|
1854 /// Returns the value of the dual solution. It should be equal to |
|
1855 /// the primal value scaled by \ref dualScale "dual scale". |
|
1856 Value dualValue() const { |
|
1857 Value sum = 0; |
|
1858 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
1859 sum += nodeValue(n); |
|
1860 } |
|
1861 for (int i = 0; i < blossomNum(); ++i) { |
|
1862 sum += blossomValue(i) * (blossomSize(i) / 2); |
|
1863 } |
|
1864 return sum; |
|
1865 } |
|
1866 |
|
1867 /// \brief Returns the value of the node. |
|
1868 /// |
|
1869 /// Returns the the value of the node. |
|
1870 Value nodeValue(const Node& n) const { |
|
1871 return (*_node_potential)[n]; |
|
1872 } |
|
1873 |
|
1874 /// \brief Returns the number of the blossoms in the basis. |
|
1875 /// |
|
1876 /// Returns the number of the blossoms in the basis. |
|
1877 /// \see BlossomIt |
|
1878 int blossomNum() const { |
|
1879 return _blossom_potential.size(); |
|
1880 } |
|
1881 |
|
1882 |
|
1883 /// \brief Returns the number of the nodes in the blossom. |
|
1884 /// |
|
1885 /// Returns the number of the nodes in the blossom. |
|
1886 int blossomSize(int k) const { |
|
1887 return _blossom_potential[k].end - _blossom_potential[k].begin; |
|
1888 } |
|
1889 |
|
1890 /// \brief Returns the value of the blossom. |
|
1891 /// |
|
1892 /// Returns the the value of the blossom. |
|
1893 /// \see BlossomIt |
|
1894 Value blossomValue(int k) const { |
|
1895 return _blossom_potential[k].value; |
|
1896 } |
|
1897 |
|
1898 /// \brief Lemon iterator for get the items of the blossom. |
|
1899 /// |
|
1900 /// Lemon iterator for get the nodes of the blossom. This class |
|
1901 /// provides a common style lemon iterator which gives back a |
|
1902 /// subset of the nodes. |
|
1903 class BlossomIt { |
|
1904 public: |
|
1905 |
|
1906 /// \brief Constructor. |
|
1907 /// |
|
1908 /// Constructor for get the nodes of the variable. |
|
1909 BlossomIt(const MaxWeightedMatching& algorithm, int variable) |
|
1910 : _algorithm(&algorithm) |
|
1911 { |
|
1912 _index = _algorithm->_blossom_potential[variable].begin; |
|
1913 _last = _algorithm->_blossom_potential[variable].end; |
|
1914 } |
|
1915 |
|
1916 /// \brief Invalid constructor. |
|
1917 /// |
|
1918 /// Invalid constructor. |
|
1919 BlossomIt(Invalid) : _index(-1) {} |
|
1920 |
|
1921 /// \brief Conversion to node. |
|
1922 /// |
|
1923 /// Conversion to node. |
|
1924 operator Node() const { |
|
1925 return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID; |
|
1926 } |
|
1927 |
|
1928 /// \brief Increment operator. |
|
1929 /// |
|
1930 /// Increment operator. |
|
1931 BlossomIt& operator++() { |
|
1932 ++_index; |
|
1933 if (_index == _last) { |
|
1934 _index = -1; |
|
1935 } |
|
1936 return *this; |
|
1937 } |
|
1938 |
|
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 } |
|
1945 |
|
1946 private: |
|
1947 const MaxWeightedMatching* _algorithm; |
|
1948 int _last; |
|
1949 int _index; |
|
1950 }; |
|
1951 |
|
1952 /// @} |
|
1953 |
|
1954 }; |
|
1955 |
|
1956 /// \ingroup matching |
|
1957 /// |
|
1958 /// \brief Weighted perfect matching in general graphs |
|
1959 /// |
|
1960 /// This class provides an efficient implementation of Edmond's |
|
1961 /// maximum weighted perfecr matching algorithm. The implementation |
|
1962 /// is based on extensive use of priority queues and provides |
|
1963 /// \f$O(nm\log(n))\f$ time complexity. |
|
1964 /// |
|
1965 /// The maximum weighted matching problem is to find undirected |
|
1966 /// arcs in the digraph with maximum overall weight and no two of |
|
1967 /// them shares their endpoints and covers all nodes. The problem |
|
1968 /// can be formulated with the next linear program: |
|
1969 /// \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] |
|
1971 /// \f[x_e \ge 0\quad \forall e\in E\f] |
|
1972 /// \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 in |
|
1974 /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in |
|
1975 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of |
|
1976 /// the nodes. |
|
1977 /// |
|
1978 /// The algorithm calculates an optimal matching and a proof of the |
|
1979 /// 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] |
|
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] |
|
1984 /// |
|
1985 /// The algorithm can be executed with \c run() or the \c init() and |
|
1986 /// then the \c start() member functions. After it the matching can |
|
1987 /// be asked with \c matching() or mate() functions. The dual |
|
1988 /// solution can be get with \c nodeValue(), \c blossomNum() and \c |
|
1989 /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt |
|
1990 /// "BlossomIt" nested class which is able to iterate on the nodes |
|
1991 /// of a blossom. If the value type is integral then the dual |
|
1992 /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4". |
|
1993 template <typename _Graph, |
|
1994 typename _WeightMap = typename _Graph::template EdgeMap<int> > |
|
1995 class MaxWeightedPerfectMatching { |
|
1996 public: |
|
1997 |
|
1998 typedef _Graph Graph; |
|
1999 typedef _WeightMap WeightMap; |
|
2000 typedef typename WeightMap::Value Value; |
|
2001 |
|
2002 /// \brief Scaling factor for dual solution |
|
2003 /// |
|
2004 /// Scaling factor for dual solution, it is equal to 4 or 1 |
|
2005 /// according to the value type. |
|
2006 static const int dualScale = |
|
2007 std::numeric_limits<Value>::is_integer ? 4 : 1; |
|
2008 |
|
2009 typedef typename Graph::template NodeMap<typename Graph::Arc> |
|
2010 MatchingMap; |
|
2011 |
|
2012 private: |
|
2013 |
|
2014 TEMPLATE_GRAPH_TYPEDEFS(Graph); |
|
2015 |
|
2016 typedef typename Graph::template NodeMap<Value> NodePotential; |
|
2017 typedef std::vector<Node> BlossomNodeList; |
|
2018 |
|
2019 struct BlossomVariable { |
|
2020 int begin, end; |
|
2021 Value value; |
|
2022 |
|
2023 BlossomVariable(int _begin, int _end, Value _value) |
|
2024 : begin(_begin), end(_end), value(_value) {} |
|
2025 |
|
2026 }; |
|
2027 |
|
2028 typedef std::vector<BlossomVariable> BlossomPotential; |
|
2029 |
|
2030 const Graph& _graph; |
|
2031 const WeightMap& _weight; |
|
2032 |
|
2033 MatchingMap* _matching; |
|
2034 |
|
2035 NodePotential* _node_potential; |
|
2036 |
|
2037 BlossomPotential _blossom_potential; |
|
2038 BlossomNodeList _blossom_node_list; |
|
2039 |
|
2040 int _node_num; |
|
2041 int _blossom_num; |
|
2042 |
|
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 typedef RangeMap<int> IntIntMap; |
|
2047 |
|
2048 enum Status { |
|
2049 EVEN = -1, MATCHED = 0, ODD = 1 |
|
2050 }; |
|
2051 |
|
2052 typedef HeapUnionFind<Value, NodeIntMap> BlossomSet; |
|
2053 struct BlossomData { |
|
2054 int tree; |
|
2055 Status status; |
|
2056 Arc pred, next; |
|
2057 Value pot, offset; |
|
2058 }; |
|
2059 |
|
2060 NodeIntMap *_blossom_index; |
|
2061 BlossomSet *_blossom_set; |
|
2062 RangeMap<BlossomData>* _blossom_data; |
|
2063 |
|
2064 NodeIntMap *_node_index; |
|
2065 ArcIntMap *_node_heap_index; |
|
2066 |
|
2067 struct NodeData { |
|
2068 |
|
2069 NodeData(ArcIntMap& node_heap_index) |
|
2070 : heap(node_heap_index) {} |
|
2071 |
|
2072 int blossom; |
|
2073 Value pot; |
|
2074 BinHeap<Value, ArcIntMap> heap; |
|
2075 std::map<int, Arc> heap_index; |
|
2076 |
|
2077 int tree; |
|
2078 }; |
|
2079 |
|
2080 RangeMap<NodeData>* _node_data; |
|
2081 |
|
2082 typedef ExtendFindEnum<IntIntMap> TreeSet; |
|
2083 |
|
2084 IntIntMap *_tree_set_index; |
|
2085 TreeSet *_tree_set; |
|
2086 |
|
2087 IntIntMap *_delta2_index; |
|
2088 BinHeap<Value, IntIntMap> *_delta2; |
|
2089 |
|
2090 EdgeIntMap *_delta3_index; |
|
2091 BinHeap<Value, EdgeIntMap> *_delta3; |
|
2092 |
|
2093 IntIntMap *_delta4_index; |
|
2094 BinHeap<Value, IntIntMap> *_delta4; |
|
2095 |
|
2096 Value _delta_sum; |
|
2097 |
|
2098 void createStructures() { |
|
2099 _node_num = countNodes(_graph); |
|
2100 _blossom_num = _node_num * 3 / 2; |
|
2101 |
|
2102 if (!_matching) { |
|
2103 _matching = new MatchingMap(_graph); |
|
2104 } |
|
2105 if (!_node_potential) { |
|
2106 _node_potential = new NodePotential(_graph); |
|
2107 } |
|
2108 if (!_blossom_set) { |
|
2109 _blossom_index = new NodeIntMap(_graph); |
|
2110 _blossom_set = new BlossomSet(*_blossom_index); |
|
2111 _blossom_data = new RangeMap<BlossomData>(_blossom_num); |
|
2112 } |
|
2113 |
|
2114 if (!_node_index) { |
|
2115 _node_index = new NodeIntMap(_graph); |
|
2116 _node_heap_index = new ArcIntMap(_graph); |
|
2117 _node_data = new RangeMap<NodeData>(_node_num, |
|
2118 NodeData(*_node_heap_index)); |
|
2119 } |
|
2120 |
|
2121 if (!_tree_set) { |
|
2122 _tree_set_index = new IntIntMap(_blossom_num); |
|
2123 _tree_set = new TreeSet(*_tree_set_index); |
|
2124 } |
|
2125 if (!_delta2) { |
|
2126 _delta2_index = new IntIntMap(_blossom_num); |
|
2127 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index); |
|
2128 } |
|
2129 if (!_delta3) { |
|
2130 _delta3_index = new EdgeIntMap(_graph); |
|
2131 _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index); |
|
2132 } |
|
2133 if (!_delta4) { |
|
2134 _delta4_index = new IntIntMap(_blossom_num); |
|
2135 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index); |
|
2136 } |
|
2137 } |
|
2138 |
|
2139 void destroyStructures() { |
|
2140 _node_num = countNodes(_graph); |
|
2141 _blossom_num = _node_num * 3 / 2; |
|
2142 |
|
2143 if (_matching) { |
|
2144 delete _matching; |
|
2145 } |
|
2146 if (_node_potential) { |
|
2147 delete _node_potential; |
|
2148 } |
|
2149 if (_blossom_set) { |
|
2150 delete _blossom_index; |
|
2151 delete _blossom_set; |
|
2152 delete _blossom_data; |
|
2153 } |
|
2154 |
|
2155 if (_node_index) { |
|
2156 delete _node_index; |
|
2157 delete _node_heap_index; |
|
2158 delete _node_data; |
|
2159 } |
|
2160 |
|
2161 if (_tree_set) { |
|
2162 delete _tree_set_index; |
|
2163 delete _tree_set; |
|
2164 } |
|
2165 if (_delta2) { |
|
2166 delete _delta2_index; |
|
2167 delete _delta2; |
|
2168 } |
|
2169 if (_delta3) { |
|
2170 delete _delta3_index; |
|
2171 delete _delta3; |
|
2172 } |
|
2173 if (_delta4) { |
|
2174 delete _delta4_index; |
|
2175 delete _delta4; |
|
2176 } |
|
2177 } |
|
2178 |
|
2179 void matchedToEven(int blossom, int tree) { |
|
2180 if (_delta2->state(blossom) == _delta2->IN_HEAP) { |
|
2181 _delta2->erase(blossom); |
|
2182 } |
|
2183 |
|
2184 if (!_blossom_set->trivial(blossom)) { |
|
2185 (*_blossom_data)[blossom].pot -= |
|
2186 2 * (_delta_sum - (*_blossom_data)[blossom].offset); |
|
2187 } |
|
2188 |
|
2189 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
|
2190 n != INVALID; ++n) { |
|
2191 |
|
2192 _blossom_set->increase(n, std::numeric_limits<Value>::max()); |
|
2193 int ni = (*_node_index)[n]; |
|
2194 |
|
2195 (*_node_data)[ni].heap.clear(); |
|
2196 (*_node_data)[ni].heap_index.clear(); |
|
2197 |
|
2198 (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset; |
|
2199 |
|
2200 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
2201 Node v = _graph.source(e); |
|
2202 int vb = _blossom_set->find(v); |
|
2203 int vi = (*_node_index)[v]; |
|
2204 |
|
2205 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - |
|
2206 dualScale * _weight[e]; |
|
2207 |
|
2208 if ((*_blossom_data)[vb].status == EVEN) { |
|
2209 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { |
|
2210 _delta3->push(e, rw / 2); |
|
2211 } |
|
2212 } else { |
|
2213 typename std::map<int, Arc>::iterator it = |
|
2214 (*_node_data)[vi].heap_index.find(tree); |
|
2215 |
|
2216 if (it != (*_node_data)[vi].heap_index.end()) { |
|
2217 if ((*_node_data)[vi].heap[it->second] > rw) { |
|
2218 (*_node_data)[vi].heap.replace(it->second, e); |
|
2219 (*_node_data)[vi].heap.decrease(e, rw); |
|
2220 it->second = e; |
|
2221 } |
|
2222 } else { |
|
2223 (*_node_data)[vi].heap.push(e, rw); |
|
2224 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); |
|
2225 } |
|
2226 |
|
2227 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { |
|
2228 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); |
|
2229 |
|
2230 if ((*_blossom_data)[vb].status == MATCHED) { |
|
2231 if (_delta2->state(vb) != _delta2->IN_HEAP) { |
|
2232 _delta2->push(vb, _blossom_set->classPrio(vb) - |
|
2233 (*_blossom_data)[vb].offset); |
|
2234 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - |
|
2235 (*_blossom_data)[vb].offset){ |
|
2236 _delta2->decrease(vb, _blossom_set->classPrio(vb) - |
|
2237 (*_blossom_data)[vb].offset); |
|
2238 } |
|
2239 } |
|
2240 } |
|
2241 } |
|
2242 } |
|
2243 } |
|
2244 (*_blossom_data)[blossom].offset = 0; |
|
2245 } |
|
2246 |
|
2247 void matchedToOdd(int blossom) { |
|
2248 if (_delta2->state(blossom) == _delta2->IN_HEAP) { |
|
2249 _delta2->erase(blossom); |
|
2250 } |
|
2251 (*_blossom_data)[blossom].offset += _delta_sum; |
|
2252 if (!_blossom_set->trivial(blossom)) { |
|
2253 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + |
|
2254 (*_blossom_data)[blossom].offset); |
|
2255 } |
|
2256 } |
|
2257 |
|
2258 void evenToMatched(int blossom, int tree) { |
|
2259 if (!_blossom_set->trivial(blossom)) { |
|
2260 (*_blossom_data)[blossom].pot += 2 * _delta_sum; |
|
2261 } |
|
2262 |
|
2263 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
|
2264 n != INVALID; ++n) { |
|
2265 int ni = (*_node_index)[n]; |
|
2266 (*_node_data)[ni].pot -= _delta_sum; |
|
2267 |
|
2268 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
2269 Node v = _graph.source(e); |
|
2270 int vb = _blossom_set->find(v); |
|
2271 int vi = (*_node_index)[v]; |
|
2272 |
|
2273 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - |
|
2274 dualScale * _weight[e]; |
|
2275 |
|
2276 if (vb == blossom) { |
|
2277 if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
2278 _delta3->erase(e); |
|
2279 } |
|
2280 } else if ((*_blossom_data)[vb].status == EVEN) { |
|
2281 |
|
2282 if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
2283 _delta3->erase(e); |
|
2284 } |
|
2285 |
|
2286 int vt = _tree_set->find(vb); |
|
2287 |
|
2288 if (vt != tree) { |
|
2289 |
|
2290 Arc r = _graph.oppositeArc(e); |
|
2291 |
|
2292 typename std::map<int, Arc>::iterator it = |
|
2293 (*_node_data)[ni].heap_index.find(vt); |
|
2294 |
|
2295 if (it != (*_node_data)[ni].heap_index.end()) { |
|
2296 if ((*_node_data)[ni].heap[it->second] > rw) { |
|
2297 (*_node_data)[ni].heap.replace(it->second, r); |
|
2298 (*_node_data)[ni].heap.decrease(r, rw); |
|
2299 it->second = r; |
|
2300 } |
|
2301 } else { |
|
2302 (*_node_data)[ni].heap.push(r, rw); |
|
2303 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); |
|
2304 } |
|
2305 |
|
2306 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { |
|
2307 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); |
|
2308 |
|
2309 if (_delta2->state(blossom) != _delta2->IN_HEAP) { |
|
2310 _delta2->push(blossom, _blossom_set->classPrio(blossom) - |
|
2311 (*_blossom_data)[blossom].offset); |
|
2312 } else if ((*_delta2)[blossom] > |
|
2313 _blossom_set->classPrio(blossom) - |
|
2314 (*_blossom_data)[blossom].offset){ |
|
2315 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - |
|
2316 (*_blossom_data)[blossom].offset); |
|
2317 } |
|
2318 } |
|
2319 } |
|
2320 } else { |
|
2321 |
|
2322 typename std::map<int, Arc>::iterator it = |
|
2323 (*_node_data)[vi].heap_index.find(tree); |
|
2324 |
|
2325 if (it != (*_node_data)[vi].heap_index.end()) { |
|
2326 (*_node_data)[vi].heap.erase(it->second); |
|
2327 (*_node_data)[vi].heap_index.erase(it); |
|
2328 if ((*_node_data)[vi].heap.empty()) { |
|
2329 _blossom_set->increase(v, std::numeric_limits<Value>::max()); |
|
2330 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { |
|
2331 _blossom_set->increase(v, (*_node_data)[vi].heap.prio()); |
|
2332 } |
|
2333 |
|
2334 if ((*_blossom_data)[vb].status == MATCHED) { |
|
2335 if (_blossom_set->classPrio(vb) == |
|
2336 std::numeric_limits<Value>::max()) { |
|
2337 _delta2->erase(vb); |
|
2338 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - |
|
2339 (*_blossom_data)[vb].offset) { |
|
2340 _delta2->increase(vb, _blossom_set->classPrio(vb) - |
|
2341 (*_blossom_data)[vb].offset); |
|
2342 } |
|
2343 } |
|
2344 } |
|
2345 } |
|
2346 } |
|
2347 } |
|
2348 } |
|
2349 |
|
2350 void oddToMatched(int blossom) { |
|
2351 (*_blossom_data)[blossom].offset -= _delta_sum; |
|
2352 |
|
2353 if (_blossom_set->classPrio(blossom) != |
|
2354 std::numeric_limits<Value>::max()) { |
|
2355 _delta2->push(blossom, _blossom_set->classPrio(blossom) - |
|
2356 (*_blossom_data)[blossom].offset); |
|
2357 } |
|
2358 |
|
2359 if (!_blossom_set->trivial(blossom)) { |
|
2360 _delta4->erase(blossom); |
|
2361 } |
|
2362 } |
|
2363 |
|
2364 void oddToEven(int blossom, int tree) { |
|
2365 if (!_blossom_set->trivial(blossom)) { |
|
2366 _delta4->erase(blossom); |
|
2367 (*_blossom_data)[blossom].pot -= |
|
2368 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); |
|
2369 } |
|
2370 |
|
2371 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
|
2372 n != INVALID; ++n) { |
|
2373 int ni = (*_node_index)[n]; |
|
2374 |
|
2375 _blossom_set->increase(n, std::numeric_limits<Value>::max()); |
|
2376 |
|
2377 (*_node_data)[ni].heap.clear(); |
|
2378 (*_node_data)[ni].heap_index.clear(); |
|
2379 (*_node_data)[ni].pot += |
|
2380 2 * _delta_sum - (*_blossom_data)[blossom].offset; |
|
2381 |
|
2382 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
2383 Node v = _graph.source(e); |
|
2384 int vb = _blossom_set->find(v); |
|
2385 int vi = (*_node_index)[v]; |
|
2386 |
|
2387 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - |
|
2388 dualScale * _weight[e]; |
|
2389 |
|
2390 if ((*_blossom_data)[vb].status == EVEN) { |
|
2391 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { |
|
2392 _delta3->push(e, rw / 2); |
|
2393 } |
|
2394 } else { |
|
2395 |
|
2396 typename std::map<int, Arc>::iterator it = |
|
2397 (*_node_data)[vi].heap_index.find(tree); |
|
2398 |
|
2399 if (it != (*_node_data)[vi].heap_index.end()) { |
|
2400 if ((*_node_data)[vi].heap[it->second] > rw) { |
|
2401 (*_node_data)[vi].heap.replace(it->second, e); |
|
2402 (*_node_data)[vi].heap.decrease(e, rw); |
|
2403 it->second = e; |
|
2404 } |
|
2405 } else { |
|
2406 (*_node_data)[vi].heap.push(e, rw); |
|
2407 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); |
|
2408 } |
|
2409 |
|
2410 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { |
|
2411 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); |
|
2412 |
|
2413 if ((*_blossom_data)[vb].status == MATCHED) { |
|
2414 if (_delta2->state(vb) != _delta2->IN_HEAP) { |
|
2415 _delta2->push(vb, _blossom_set->classPrio(vb) - |
|
2416 (*_blossom_data)[vb].offset); |
|
2417 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - |
|
2418 (*_blossom_data)[vb].offset) { |
|
2419 _delta2->decrease(vb, _blossom_set->classPrio(vb) - |
|
2420 (*_blossom_data)[vb].offset); |
|
2421 } |
|
2422 } |
|
2423 } |
|
2424 } |
|
2425 } |
|
2426 } |
|
2427 (*_blossom_data)[blossom].offset = 0; |
|
2428 } |
|
2429 |
|
2430 void alternatePath(int even, int tree) { |
|
2431 int odd; |
|
2432 |
|
2433 evenToMatched(even, tree); |
|
2434 (*_blossom_data)[even].status = MATCHED; |
|
2435 |
|
2436 while ((*_blossom_data)[even].pred != INVALID) { |
|
2437 odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred)); |
|
2438 (*_blossom_data)[odd].status = MATCHED; |
|
2439 oddToMatched(odd); |
|
2440 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred; |
|
2441 |
|
2442 even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred)); |
|
2443 (*_blossom_data)[even].status = MATCHED; |
|
2444 evenToMatched(even, tree); |
|
2445 (*_blossom_data)[even].next = |
|
2446 _graph.oppositeArc((*_blossom_data)[odd].pred); |
|
2447 } |
|
2448 |
|
2449 } |
|
2450 |
|
2451 void destroyTree(int tree) { |
|
2452 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) { |
|
2453 if ((*_blossom_data)[b].status == EVEN) { |
|
2454 (*_blossom_data)[b].status = MATCHED; |
|
2455 evenToMatched(b, tree); |
|
2456 } else if ((*_blossom_data)[b].status == ODD) { |
|
2457 (*_blossom_data)[b].status = MATCHED; |
|
2458 oddToMatched(b); |
|
2459 } |
|
2460 } |
|
2461 _tree_set->eraseClass(tree); |
|
2462 } |
|
2463 |
|
2464 void augmentOnArc(const Edge& arc) { |
|
2465 |
|
2466 int left = _blossom_set->find(_graph.u(arc)); |
|
2467 int right = _blossom_set->find(_graph.v(arc)); |
|
2468 |
|
2469 int left_tree = _tree_set->find(left); |
|
2470 alternatePath(left, left_tree); |
|
2471 destroyTree(left_tree); |
|
2472 |
|
2473 int right_tree = _tree_set->find(right); |
|
2474 alternatePath(right, right_tree); |
|
2475 destroyTree(right_tree); |
|
2476 |
|
2477 (*_blossom_data)[left].next = _graph.direct(arc, true); |
|
2478 (*_blossom_data)[right].next = _graph.direct(arc, false); |
|
2479 } |
|
2480 |
|
2481 void extendOnArc(const Arc& arc) { |
|
2482 int base = _blossom_set->find(_graph.target(arc)); |
|
2483 int tree = _tree_set->find(base); |
|
2484 |
|
2485 int odd = _blossom_set->find(_graph.source(arc)); |
|
2486 _tree_set->insert(odd, tree); |
|
2487 (*_blossom_data)[odd].status = ODD; |
|
2488 matchedToOdd(odd); |
|
2489 (*_blossom_data)[odd].pred = arc; |
|
2490 |
|
2491 int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next)); |
|
2492 (*_blossom_data)[even].pred = (*_blossom_data)[even].next; |
|
2493 _tree_set->insert(even, tree); |
|
2494 (*_blossom_data)[even].status = EVEN; |
|
2495 matchedToEven(even, tree); |
|
2496 } |
|
2497 |
|
2498 void shrinkOnArc(const Edge& edge, int tree) { |
|
2499 int nca = -1; |
|
2500 std::vector<int> left_path, right_path; |
|
2501 |
|
2502 { |
|
2503 std::set<int> left_set, right_set; |
|
2504 int left = _blossom_set->find(_graph.u(edge)); |
|
2505 left_path.push_back(left); |
|
2506 left_set.insert(left); |
|
2507 |
|
2508 int right = _blossom_set->find(_graph.v(edge)); |
|
2509 right_path.push_back(right); |
|
2510 right_set.insert(right); |
|
2511 |
|
2512 while (true) { |
|
2513 |
|
2514 if ((*_blossom_data)[left].pred == INVALID) break; |
|
2515 |
|
2516 left = |
|
2517 _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); |
|
2518 left_path.push_back(left); |
|
2519 left = |
|
2520 _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); |
|
2521 left_path.push_back(left); |
|
2522 |
|
2523 left_set.insert(left); |
|
2524 |
|
2525 if (right_set.find(left) != right_set.end()) { |
|
2526 nca = left; |
|
2527 break; |
|
2528 } |
|
2529 |
|
2530 if ((*_blossom_data)[right].pred == INVALID) break; |
|
2531 |
|
2532 right = |
|
2533 _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); |
|
2534 right_path.push_back(right); |
|
2535 right = |
|
2536 _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); |
|
2537 right_path.push_back(right); |
|
2538 |
|
2539 right_set.insert(right); |
|
2540 |
|
2541 if (left_set.find(right) != left_set.end()) { |
|
2542 nca = right; |
|
2543 break; |
|
2544 } |
|
2545 |
|
2546 } |
|
2547 |
|
2548 if (nca == -1) { |
|
2549 if ((*_blossom_data)[left].pred == INVALID) { |
|
2550 nca = right; |
|
2551 while (left_set.find(nca) == left_set.end()) { |
|
2552 nca = |
|
2553 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); |
|
2554 right_path.push_back(nca); |
|
2555 nca = |
|
2556 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); |
|
2557 right_path.push_back(nca); |
|
2558 } |
|
2559 } else { |
|
2560 nca = left; |
|
2561 while (right_set.find(nca) == right_set.end()) { |
|
2562 nca = |
|
2563 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); |
|
2564 left_path.push_back(nca); |
|
2565 nca = |
|
2566 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); |
|
2567 left_path.push_back(nca); |
|
2568 } |
|
2569 } |
|
2570 } |
|
2571 } |
|
2572 |
|
2573 std::vector<int> subblossoms; |
|
2574 Arc prev; |
|
2575 |
|
2576 prev = _graph.direct(edge, true); |
|
2577 for (int i = 0; left_path[i] != nca; i += 2) { |
|
2578 subblossoms.push_back(left_path[i]); |
|
2579 (*_blossom_data)[left_path[i]].next = prev; |
|
2580 _tree_set->erase(left_path[i]); |
|
2581 |
|
2582 subblossoms.push_back(left_path[i + 1]); |
|
2583 (*_blossom_data)[left_path[i + 1]].status = EVEN; |
|
2584 oddToEven(left_path[i + 1], tree); |
|
2585 _tree_set->erase(left_path[i + 1]); |
|
2586 prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred); |
|
2587 } |
|
2588 |
|
2589 int k = 0; |
|
2590 while (right_path[k] != nca) ++k; |
|
2591 |
|
2592 subblossoms.push_back(nca); |
|
2593 (*_blossom_data)[nca].next = prev; |
|
2594 |
|
2595 for (int i = k - 2; i >= 0; i -= 2) { |
|
2596 subblossoms.push_back(right_path[i + 1]); |
|
2597 (*_blossom_data)[right_path[i + 1]].status = EVEN; |
|
2598 oddToEven(right_path[i + 1], tree); |
|
2599 _tree_set->erase(right_path[i + 1]); |
|
2600 |
|
2601 (*_blossom_data)[right_path[i + 1]].next = |
|
2602 (*_blossom_data)[right_path[i + 1]].pred; |
|
2603 |
|
2604 subblossoms.push_back(right_path[i]); |
|
2605 _tree_set->erase(right_path[i]); |
|
2606 } |
|
2607 |
|
2608 int surface = |
|
2609 _blossom_set->join(subblossoms.begin(), subblossoms.end()); |
|
2610 |
|
2611 for (int i = 0; i < int(subblossoms.size()); ++i) { |
|
2612 if (!_blossom_set->trivial(subblossoms[i])) { |
|
2613 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum; |
|
2614 } |
|
2615 (*_blossom_data)[subblossoms[i]].status = MATCHED; |
|
2616 } |
|
2617 |
|
2618 (*_blossom_data)[surface].pot = -2 * _delta_sum; |
|
2619 (*_blossom_data)[surface].offset = 0; |
|
2620 (*_blossom_data)[surface].status = EVEN; |
|
2621 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred; |
|
2622 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred; |
|
2623 |
|
2624 _tree_set->insert(surface, tree); |
|
2625 _tree_set->erase(nca); |
|
2626 } |
|
2627 |
|
2628 void splitBlossom(int blossom) { |
|
2629 Arc next = (*_blossom_data)[blossom].next; |
|
2630 Arc pred = (*_blossom_data)[blossom].pred; |
|
2631 |
|
2632 int tree = _tree_set->find(blossom); |
|
2633 |
|
2634 (*_blossom_data)[blossom].status = MATCHED; |
|
2635 oddToMatched(blossom); |
|
2636 if (_delta2->state(blossom) == _delta2->IN_HEAP) { |
|
2637 _delta2->erase(blossom); |
|
2638 } |
|
2639 |
|
2640 std::vector<int> subblossoms; |
|
2641 _blossom_set->split(blossom, std::back_inserter(subblossoms)); |
|
2642 |
|
2643 Value offset = (*_blossom_data)[blossom].offset; |
|
2644 int b = _blossom_set->find(_graph.source(pred)); |
|
2645 int d = _blossom_set->find(_graph.source(next)); |
|
2646 |
|
2647 int ib = -1, id = -1; |
|
2648 for (int i = 0; i < int(subblossoms.size()); ++i) { |
|
2649 if (subblossoms[i] == b) ib = i; |
|
2650 if (subblossoms[i] == d) id = i; |
|
2651 |
|
2652 (*_blossom_data)[subblossoms[i]].offset = offset; |
|
2653 if (!_blossom_set->trivial(subblossoms[i])) { |
|
2654 (*_blossom_data)[subblossoms[i]].pot -= 2 * offset; |
|
2655 } |
|
2656 if (_blossom_set->classPrio(subblossoms[i]) != |
|
2657 std::numeric_limits<Value>::max()) { |
|
2658 _delta2->push(subblossoms[i], |
|
2659 _blossom_set->classPrio(subblossoms[i]) - |
|
2660 (*_blossom_data)[subblossoms[i]].offset); |
|
2661 } |
|
2662 } |
|
2663 |
|
2664 if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) { |
|
2665 for (int i = (id + 1) % subblossoms.size(); |
|
2666 i != ib; i = (i + 2) % subblossoms.size()) { |
|
2667 int sb = subblossoms[i]; |
|
2668 int tb = subblossoms[(i + 1) % subblossoms.size()]; |
|
2669 (*_blossom_data)[sb].next = |
|
2670 _graph.oppositeArc((*_blossom_data)[tb].next); |
|
2671 } |
|
2672 |
|
2673 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) { |
|
2674 int sb = subblossoms[i]; |
|
2675 int tb = subblossoms[(i + 1) % subblossoms.size()]; |
|
2676 int ub = subblossoms[(i + 2) % subblossoms.size()]; |
|
2677 |
|
2678 (*_blossom_data)[sb].status = ODD; |
|
2679 matchedToOdd(sb); |
|
2680 _tree_set->insert(sb, tree); |
|
2681 (*_blossom_data)[sb].pred = pred; |
|
2682 (*_blossom_data)[sb].next = |
|
2683 _graph.oppositeArc((*_blossom_data)[tb].next); |
|
2684 |
|
2685 pred = (*_blossom_data)[ub].next; |
|
2686 |
|
2687 (*_blossom_data)[tb].status = EVEN; |
|
2688 matchedToEven(tb, tree); |
|
2689 _tree_set->insert(tb, tree); |
|
2690 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next; |
|
2691 } |
|
2692 |
|
2693 (*_blossom_data)[subblossoms[id]].status = ODD; |
|
2694 matchedToOdd(subblossoms[id]); |
|
2695 _tree_set->insert(subblossoms[id], tree); |
|
2696 (*_blossom_data)[subblossoms[id]].next = next; |
|
2697 (*_blossom_data)[subblossoms[id]].pred = pred; |
|
2698 |
|
2699 } else { |
|
2700 |
|
2701 for (int i = (ib + 1) % subblossoms.size(); |
|
2702 i != id; i = (i + 2) % subblossoms.size()) { |
|
2703 int sb = subblossoms[i]; |
|
2704 int tb = subblossoms[(i + 1) % subblossoms.size()]; |
|
2705 (*_blossom_data)[sb].next = |
|
2706 _graph.oppositeArc((*_blossom_data)[tb].next); |
|
2707 } |
|
2708 |
|
2709 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) { |
|
2710 int sb = subblossoms[i]; |
|
2711 int tb = subblossoms[(i + 1) % subblossoms.size()]; |
|
2712 int ub = subblossoms[(i + 2) % subblossoms.size()]; |
|
2713 |
|
2714 (*_blossom_data)[sb].status = ODD; |
|
2715 matchedToOdd(sb); |
|
2716 _tree_set->insert(sb, tree); |
|
2717 (*_blossom_data)[sb].next = next; |
|
2718 (*_blossom_data)[sb].pred = |
|
2719 _graph.oppositeArc((*_blossom_data)[tb].next); |
|
2720 |
|
2721 (*_blossom_data)[tb].status = EVEN; |
|
2722 matchedToEven(tb, tree); |
|
2723 _tree_set->insert(tb, tree); |
|
2724 (*_blossom_data)[tb].pred = |
|
2725 (*_blossom_data)[tb].next = |
|
2726 _graph.oppositeArc((*_blossom_data)[ub].next); |
|
2727 next = (*_blossom_data)[ub].next; |
|
2728 } |
|
2729 |
|
2730 (*_blossom_data)[subblossoms[ib]].status = ODD; |
|
2731 matchedToOdd(subblossoms[ib]); |
|
2732 _tree_set->insert(subblossoms[ib], tree); |
|
2733 (*_blossom_data)[subblossoms[ib]].next = next; |
|
2734 (*_blossom_data)[subblossoms[ib]].pred = pred; |
|
2735 } |
|
2736 _tree_set->erase(blossom); |
|
2737 } |
|
2738 |
|
2739 void extractBlossom(int blossom, const Node& base, const Arc& matching) { |
|
2740 if (_blossom_set->trivial(blossom)) { |
|
2741 int bi = (*_node_index)[base]; |
|
2742 Value pot = (*_node_data)[bi].pot; |
|
2743 |
|
2744 _matching->set(base, matching); |
|
2745 _blossom_node_list.push_back(base); |
|
2746 _node_potential->set(base, pot); |
|
2747 } else { |
|
2748 |
|
2749 Value pot = (*_blossom_data)[blossom].pot; |
|
2750 int bn = _blossom_node_list.size(); |
|
2751 |
|
2752 std::vector<int> subblossoms; |
|
2753 _blossom_set->split(blossom, std::back_inserter(subblossoms)); |
|
2754 int b = _blossom_set->find(base); |
|
2755 int ib = -1; |
|
2756 for (int i = 0; i < int(subblossoms.size()); ++i) { |
|
2757 if (subblossoms[i] == b) { ib = i; break; } |
|
2758 } |
|
2759 |
|
2760 for (int i = 1; i < int(subblossoms.size()); i += 2) { |
|
2761 int sb = subblossoms[(ib + i) % subblossoms.size()]; |
|
2762 int tb = subblossoms[(ib + i + 1) % subblossoms.size()]; |
|
2763 |
|
2764 Arc m = (*_blossom_data)[tb].next; |
|
2765 extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m)); |
|
2766 extractBlossom(tb, _graph.source(m), m); |
|
2767 } |
|
2768 extractBlossom(subblossoms[ib], base, matching); |
|
2769 |
|
2770 int en = _blossom_node_list.size(); |
|
2771 |
|
2772 _blossom_potential.push_back(BlossomVariable(bn, en, pot)); |
|
2773 } |
|
2774 } |
|
2775 |
|
2776 void extractMatching() { |
|
2777 std::vector<int> blossoms; |
|
2778 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) { |
|
2779 blossoms.push_back(c); |
|
2780 } |
|
2781 |
|
2782 for (int i = 0; i < int(blossoms.size()); ++i) { |
|
2783 |
|
2784 Value offset = (*_blossom_data)[blossoms[i]].offset; |
|
2785 (*_blossom_data)[blossoms[i]].pot += 2 * offset; |
|
2786 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); |
|
2787 n != INVALID; ++n) { |
|
2788 (*_node_data)[(*_node_index)[n]].pot -= offset; |
|
2789 } |
|
2790 |
|
2791 Arc matching = (*_blossom_data)[blossoms[i]].next; |
|
2792 Node base = _graph.source(matching); |
|
2793 extractBlossom(blossoms[i], base, matching); |
|
2794 } |
|
2795 } |
|
2796 |
|
2797 public: |
|
2798 |
|
2799 /// \brief Constructor |
|
2800 /// |
|
2801 /// Constructor. |
|
2802 MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight) |
|
2803 : _graph(graph), _weight(weight), _matching(0), |
|
2804 _node_potential(0), _blossom_potential(), _blossom_node_list(), |
|
2805 _node_num(0), _blossom_num(0), |
|
2806 |
|
2807 _blossom_index(0), _blossom_set(0), _blossom_data(0), |
|
2808 _node_index(0), _node_heap_index(0), _node_data(0), |
|
2809 _tree_set_index(0), _tree_set(0), |
|
2810 |
|
2811 _delta2_index(0), _delta2(0), |
|
2812 _delta3_index(0), _delta3(0), |
|
2813 _delta4_index(0), _delta4(0), |
|
2814 |
|
2815 _delta_sum() {} |
|
2816 |
|
2817 ~MaxWeightedPerfectMatching() { |
|
2818 destroyStructures(); |
|
2819 } |
|
2820 |
|
2821 /// \name Execution control |
|
2822 /// The simplest way to execute the algorithm is to use the member |
|
2823 /// \c run() member function. |
|
2824 |
|
2825 ///@{ |
|
2826 |
|
2827 /// \brief Initialize the algorithm |
|
2828 /// |
|
2829 /// Initialize the algorithm |
|
2830 void init() { |
|
2831 createStructures(); |
|
2832 |
|
2833 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
2834 _node_heap_index->set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP); |
|
2835 } |
|
2836 for (EdgeIt e(_graph); e != INVALID; ++e) { |
|
2837 _delta3_index->set(e, _delta3->PRE_HEAP); |
|
2838 } |
|
2839 for (int i = 0; i < _blossom_num; ++i) { |
|
2840 _delta2_index->set(i, _delta2->PRE_HEAP); |
|
2841 _delta4_index->set(i, _delta4->PRE_HEAP); |
|
2842 } |
|
2843 |
|
2844 int index = 0; |
|
2845 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
2846 Value max = - std::numeric_limits<Value>::max(); |
|
2847 for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
|
2848 if (_graph.target(e) == n) continue; |
|
2849 if ((dualScale * _weight[e]) / 2 > max) { |
|
2850 max = (dualScale * _weight[e]) / 2; |
|
2851 } |
|
2852 } |
|
2853 _node_index->set(n, index); |
|
2854 (*_node_data)[index].pot = max; |
|
2855 int blossom = |
|
2856 _blossom_set->insert(n, std::numeric_limits<Value>::max()); |
|
2857 |
|
2858 _tree_set->insert(blossom); |
|
2859 |
|
2860 (*_blossom_data)[blossom].status = EVEN; |
|
2861 (*_blossom_data)[blossom].pred = INVALID; |
|
2862 (*_blossom_data)[blossom].next = INVALID; |
|
2863 (*_blossom_data)[blossom].pot = 0; |
|
2864 (*_blossom_data)[blossom].offset = 0; |
|
2865 ++index; |
|
2866 } |
|
2867 for (EdgeIt e(_graph); e != INVALID; ++e) { |
|
2868 int si = (*_node_index)[_graph.u(e)]; |
|
2869 int ti = (*_node_index)[_graph.v(e)]; |
|
2870 if (_graph.u(e) != _graph.v(e)) { |
|
2871 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - |
|
2872 dualScale * _weight[e]) / 2); |
|
2873 } |
|
2874 } |
|
2875 } |
|
2876 |
|
2877 /// \brief Starts the algorithm |
|
2878 /// |
|
2879 /// Starts the algorithm |
|
2880 bool start() { |
|
2881 enum OpType { |
|
2882 D2, D3, D4 |
|
2883 }; |
|
2884 |
|
2885 int unmatched = _node_num; |
|
2886 while (unmatched > 0) { |
|
2887 Value d2 = !_delta2->empty() ? |
|
2888 _delta2->prio() : std::numeric_limits<Value>::max(); |
|
2889 |
|
2890 Value d3 = !_delta3->empty() ? |
|
2891 _delta3->prio() : std::numeric_limits<Value>::max(); |
|
2892 |
|
2893 Value d4 = !_delta4->empty() ? |
|
2894 _delta4->prio() : std::numeric_limits<Value>::max(); |
|
2895 |
|
2896 _delta_sum = d2; OpType ot = D2; |
|
2897 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } |
|
2898 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } |
|
2899 |
|
2900 if (_delta_sum == std::numeric_limits<Value>::max()) { |
|
2901 return false; |
|
2902 } |
|
2903 |
|
2904 switch (ot) { |
|
2905 case D2: |
|
2906 { |
|
2907 int blossom = _delta2->top(); |
|
2908 Node n = _blossom_set->classTop(blossom); |
|
2909 Arc e = (*_node_data)[(*_node_index)[n]].heap.top(); |
|
2910 extendOnArc(e); |
|
2911 } |
|
2912 break; |
|
2913 case D3: |
|
2914 { |
|
2915 Edge e = _delta3->top(); |
|
2916 |
|
2917 int left_blossom = _blossom_set->find(_graph.u(e)); |
|
2918 int right_blossom = _blossom_set->find(_graph.v(e)); |
|
2919 |
|
2920 if (left_blossom == right_blossom) { |
|
2921 _delta3->pop(); |
|
2922 } else { |
|
2923 int left_tree = _tree_set->find(left_blossom); |
|
2924 int right_tree = _tree_set->find(right_blossom); |
|
2925 |
|
2926 if (left_tree == right_tree) { |
|
2927 shrinkOnArc(e, left_tree); |
|
2928 } else { |
|
2929 augmentOnArc(e); |
|
2930 unmatched -= 2; |
|
2931 } |
|
2932 } |
|
2933 } break; |
|
2934 case D4: |
|
2935 splitBlossom(_delta4->top()); |
|
2936 break; |
|
2937 } |
|
2938 } |
|
2939 extractMatching(); |
|
2940 return true; |
|
2941 } |
|
2942 |
|
2943 /// \brief Runs %MaxWeightedPerfectMatching algorithm. |
|
2944 /// |
|
2945 /// This method runs the %MaxWeightedPerfectMatching algorithm. |
|
2946 /// |
|
2947 /// \note mwm.run() is just a shortcut of the following code. |
|
2948 /// \code |
|
2949 /// mwm.init(); |
|
2950 /// mwm.start(); |
|
2951 /// \endcode |
|
2952 bool run() { |
|
2953 init(); |
|
2954 return start(); |
|
2955 } |
|
2956 |
|
2957 /// @} |
|
2958 |
|
2959 /// \name Primal solution |
|
2960 /// Functions for get the primal solution, ie. the matching. |
|
2961 |
|
2962 /// @{ |
|
2963 |
|
2964 /// \brief Returns the matching value. |
|
2965 /// |
|
2966 /// Returns the matching value. |
|
2967 Value matchingValue() const { |
|
2968 Value sum = 0; |
|
2969 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
2970 if ((*_matching)[n] != INVALID) { |
|
2971 sum += _weight[(*_matching)[n]]; |
|
2972 } |
|
2973 } |
|
2974 return sum /= 2; |
|
2975 } |
|
2976 |
|
2977 /// \brief Returns true when the arc is in the matching. |
|
2978 /// |
|
2979 /// Returns true when the arc is 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. |
|
2987 Arc matching(const Node& node) const { |
|
2988 return (*_matching)[node]; |
|
2989 } |
|
2990 |
|
2991 /// \brief Returns the mate of the node. |
|
2992 /// |
|
2993 /// Returns the adjancent node in a mathcing arc. |
|
2994 Node mate(const Node& node) const { |
|
2995 return _graph.target((*_matching)[node]); |
|
2996 } |
|
2997 |
|
2998 /// @} |
|
2999 |
|
3000 /// \name Dual solution |
|
3001 /// Functions for get the dual solution. |
|
3002 |
|
3003 /// @{ |
|
3004 |
|
3005 /// \brief Returns the value of the dual solution. |
|
3006 /// |
|
3007 /// Returns the value of the dual solution. It should be equal to |
|
3008 /// the primal value scaled by \ref dualScale "dual scale". |
|
3009 Value dualValue() const { |
|
3010 Value sum = 0; |
|
3011 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
3012 sum += nodeValue(n); |
|
3013 } |
|
3014 for (int i = 0; i < blossomNum(); ++i) { |
|
3015 sum += blossomValue(i) * (blossomSize(i) / 2); |
|
3016 } |
|
3017 return sum; |
|
3018 } |
|
3019 |
|
3020 /// \brief Returns the value of the node. |
|
3021 /// |
|
3022 /// Returns the the value of the node. |
|
3023 Value nodeValue(const Node& n) const { |
|
3024 return (*_node_potential)[n]; |
|
3025 } |
|
3026 |
|
3027 /// \brief Returns the number of the blossoms in the basis. |
|
3028 /// |
|
3029 /// Returns the number of the blossoms in the basis. |
|
3030 /// \see BlossomIt |
|
3031 int blossomNum() const { |
|
3032 return _blossom_potential.size(); |
|
3033 } |
|
3034 |
|
3035 |
|
3036 /// \brief Returns the number of the nodes in the blossom. |
|
3037 /// |
|
3038 /// Returns the number of the nodes in the blossom. |
|
3039 int blossomSize(int k) const { |
|
3040 return _blossom_potential[k].end - _blossom_potential[k].begin; |
|
3041 } |
|
3042 |
|
3043 /// \brief Returns the value of the blossom. |
|
3044 /// |
|
3045 /// Returns the the value of the blossom. |
|
3046 /// \see BlossomIt |
|
3047 Value blossomValue(int k) const { |
|
3048 return _blossom_potential[k].value; |
|
3049 } |
|
3050 |
|
3051 /// \brief Lemon iterator for get the items of the blossom. |
|
3052 /// |
|
3053 /// Lemon iterator for get the nodes of the blossom. This class |
|
3054 /// provides a common style lemon iterator which gives back a |
|
3055 /// subset of the nodes. |
|
3056 class BlossomIt { |
|
3057 public: |
|
3058 |
|
3059 /// \brief Constructor. |
|
3060 /// |
|
3061 /// Constructor for get the nodes of the variable. |
|
3062 BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable) |
|
3063 : _algorithm(&algorithm) |
|
3064 { |
|
3065 _index = _algorithm->_blossom_potential[variable].begin; |
|
3066 _last = _algorithm->_blossom_potential[variable].end; |
|
3067 } |
|
3068 |
|
3069 /// \brief Invalid constructor. |
|
3070 /// |
|
3071 /// Invalid constructor. |
|
3072 BlossomIt(Invalid) : _index(-1) {} |
|
3073 |
|
3074 /// \brief Conversion to node. |
|
3075 /// |
|
3076 /// Conversion to node. |
|
3077 operator Node() const { |
|
3078 return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID; |
|
3079 } |
|
3080 |
|
3081 /// \brief Increment operator. |
|
3082 /// |
|
3083 /// Increment operator. |
|
3084 BlossomIt& operator++() { |
|
3085 ++_index; |
|
3086 if (_index == _last) { |
|
3087 _index = -1; |
|
3088 } |
|
3089 return *this; |
|
3090 } |
|
3091 |
|
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 } |
|
3098 |
|
3099 private: |
|
3100 const MaxWeightedPerfectMatching* _algorithm; |
|
3101 int _last; |
|
3102 int _index; |
|
3103 }; |
|
3104 |
|
3105 /// @} |
|
3106 |
|
3107 }; |
|
3108 |
|
3109 |
|
3110 } //END OF NAMESPACE LEMON |
|
3111 |
|
3112 #endif //LEMON_MAX_MATCHING_H |