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-2009 |
|
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 general graphs. |
|
35 |
|
36 namespace lemon { |
|
37 |
|
38 /// \ingroup matching |
|
39 /// |
|
40 /// \brief Edmonds' alternating forest maximum matching algorithm. |
|
41 /// |
|
42 /// This class implements Edmonds' alternating forest matching |
|
43 /// algorithm. The algorithm can be started from an arbitrary initial |
|
44 /// matching (the default is the empty one) |
|
45 /// |
|
46 /// The dual solution of the problem is a map of the nodes to |
|
47 /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c |
|
48 /// MATCHED/C showing the Gallai-Edmonds decomposition of the |
|
49 /// graph. The nodes in \c EVEN/D induce a graph with |
|
50 /// factor-critical components, the nodes in \c ODD/A form the |
|
51 /// barrier, and the nodes in \c MATCHED/C induce a graph having a |
|
52 /// perfect matching. The number of the factor-critical components |
|
53 /// minus the number of barrier nodes is a lower bound on the |
|
54 /// unmatched nodes, and the matching is optimal if and only if this bound is |
|
55 /// tight. This decomposition can be attained by calling \c |
|
56 /// decomposition() after running the algorithm. |
|
57 /// |
|
58 /// \param GR The graph type the algorithm runs on. |
|
59 template <typename GR> |
|
60 class MaxMatching { |
|
61 public: |
|
62 |
|
63 typedef GR Graph; |
|
64 typedef typename Graph::template NodeMap<typename Graph::Arc> |
|
65 MatchingMap; |
|
66 |
|
67 ///\brief Indicates the Gallai-Edmonds decomposition of the graph. |
|
68 /// |
|
69 ///Indicates the Gallai-Edmonds decomposition of the graph. The |
|
70 ///nodes with Status \c EVEN/D induce a graph with factor-critical |
|
71 ///components, the nodes in \c ODD/A form the canonical barrier, |
|
72 ///and the nodes in \c MATCHED/C induce a graph having a perfect |
|
73 ///matching. |
|
74 enum Status { |
|
75 EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2 |
|
76 }; |
|
77 |
|
78 typedef typename Graph::template NodeMap<Status> StatusMap; |
|
79 |
|
80 private: |
|
81 |
|
82 TEMPLATE_GRAPH_TYPEDEFS(Graph); |
|
83 |
|
84 typedef UnionFindEnum<IntNodeMap> BlossomSet; |
|
85 typedef ExtendFindEnum<IntNodeMap> TreeSet; |
|
86 typedef RangeMap<Node> NodeIntMap; |
|
87 typedef MatchingMap EarMap; |
|
88 typedef std::vector<Node> NodeQueue; |
|
89 |
|
90 const Graph& _graph; |
|
91 MatchingMap* _matching; |
|
92 StatusMap* _status; |
|
93 |
|
94 EarMap* _ear; |
|
95 |
|
96 IntNodeMap* _blossom_set_index; |
|
97 BlossomSet* _blossom_set; |
|
98 NodeIntMap* _blossom_rep; |
|
99 |
|
100 IntNodeMap* _tree_set_index; |
|
101 TreeSet* _tree_set; |
|
102 |
|
103 NodeQueue _node_queue; |
|
104 int _process, _postpone, _last; |
|
105 |
|
106 int _node_num; |
|
107 |
|
108 private: |
|
109 |
|
110 void createStructures() { |
|
111 _node_num = countNodes(_graph); |
|
112 if (!_matching) { |
|
113 _matching = new MatchingMap(_graph); |
|
114 } |
|
115 if (!_status) { |
|
116 _status = new StatusMap(_graph); |
|
117 } |
|
118 if (!_ear) { |
|
119 _ear = new EarMap(_graph); |
|
120 } |
|
121 if (!_blossom_set) { |
|
122 _blossom_set_index = new IntNodeMap(_graph); |
|
123 _blossom_set = new BlossomSet(*_blossom_set_index); |
|
124 } |
|
125 if (!_blossom_rep) { |
|
126 _blossom_rep = new NodeIntMap(_node_num); |
|
127 } |
|
128 if (!_tree_set) { |
|
129 _tree_set_index = new IntNodeMap(_graph); |
|
130 _tree_set = new TreeSet(*_tree_set_index); |
|
131 } |
|
132 _node_queue.resize(_node_num); |
|
133 } |
|
134 |
|
135 void destroyStructures() { |
|
136 if (_matching) { |
|
137 delete _matching; |
|
138 } |
|
139 if (_status) { |
|
140 delete _status; |
|
141 } |
|
142 if (_ear) { |
|
143 delete _ear; |
|
144 } |
|
145 if (_blossom_set) { |
|
146 delete _blossom_set; |
|
147 delete _blossom_set_index; |
|
148 } |
|
149 if (_blossom_rep) { |
|
150 delete _blossom_rep; |
|
151 } |
|
152 if (_tree_set) { |
|
153 delete _tree_set_index; |
|
154 delete _tree_set; |
|
155 } |
|
156 } |
|
157 |
|
158 void processDense(const Node& n) { |
|
159 _process = _postpone = _last = 0; |
|
160 _node_queue[_last++] = n; |
|
161 |
|
162 while (_process != _last) { |
|
163 Node u = _node_queue[_process++]; |
|
164 for (OutArcIt a(_graph, u); a != INVALID; ++a) { |
|
165 Node v = _graph.target(a); |
|
166 if ((*_status)[v] == MATCHED) { |
|
167 extendOnArc(a); |
|
168 } else if ((*_status)[v] == UNMATCHED) { |
|
169 augmentOnArc(a); |
|
170 return; |
|
171 } |
|
172 } |
|
173 } |
|
174 |
|
175 while (_postpone != _last) { |
|
176 Node u = _node_queue[_postpone++]; |
|
177 |
|
178 for (OutArcIt a(_graph, u); a != INVALID ; ++a) { |
|
179 Node v = _graph.target(a); |
|
180 |
|
181 if ((*_status)[v] == EVEN) { |
|
182 if (_blossom_set->find(u) != _blossom_set->find(v)) { |
|
183 shrinkOnEdge(a); |
|
184 } |
|
185 } |
|
186 |
|
187 while (_process != _last) { |
|
188 Node w = _node_queue[_process++]; |
|
189 for (OutArcIt b(_graph, w); b != INVALID; ++b) { |
|
190 Node x = _graph.target(b); |
|
191 if ((*_status)[x] == MATCHED) { |
|
192 extendOnArc(b); |
|
193 } else if ((*_status)[x] == UNMATCHED) { |
|
194 augmentOnArc(b); |
|
195 return; |
|
196 } |
|
197 } |
|
198 } |
|
199 } |
|
200 } |
|
201 } |
|
202 |
|
203 void processSparse(const Node& n) { |
|
204 _process = _last = 0; |
|
205 _node_queue[_last++] = n; |
|
206 while (_process != _last) { |
|
207 Node u = _node_queue[_process++]; |
|
208 for (OutArcIt a(_graph, u); a != INVALID; ++a) { |
|
209 Node v = _graph.target(a); |
|
210 |
|
211 if ((*_status)[v] == EVEN) { |
|
212 if (_blossom_set->find(u) != _blossom_set->find(v)) { |
|
213 shrinkOnEdge(a); |
|
214 } |
|
215 } else if ((*_status)[v] == MATCHED) { |
|
216 extendOnArc(a); |
|
217 } else if ((*_status)[v] == UNMATCHED) { |
|
218 augmentOnArc(a); |
|
219 return; |
|
220 } |
|
221 } |
|
222 } |
|
223 } |
|
224 |
|
225 void shrinkOnEdge(const Edge& e) { |
|
226 Node nca = INVALID; |
|
227 |
|
228 { |
|
229 std::set<Node> left_set, right_set; |
|
230 |
|
231 Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))]; |
|
232 left_set.insert(left); |
|
233 |
|
234 Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))]; |
|
235 right_set.insert(right); |
|
236 |
|
237 while (true) { |
|
238 if ((*_matching)[left] == INVALID) break; |
|
239 left = _graph.target((*_matching)[left]); |
|
240 left = (*_blossom_rep)[_blossom_set-> |
|
241 find(_graph.target((*_ear)[left]))]; |
|
242 if (right_set.find(left) != right_set.end()) { |
|
243 nca = left; |
|
244 break; |
|
245 } |
|
246 left_set.insert(left); |
|
247 |
|
248 if ((*_matching)[right] == INVALID) break; |
|
249 right = _graph.target((*_matching)[right]); |
|
250 right = (*_blossom_rep)[_blossom_set-> |
|
251 find(_graph.target((*_ear)[right]))]; |
|
252 if (left_set.find(right) != left_set.end()) { |
|
253 nca = right; |
|
254 break; |
|
255 } |
|
256 right_set.insert(right); |
|
257 } |
|
258 |
|
259 if (nca == INVALID) { |
|
260 if ((*_matching)[left] == INVALID) { |
|
261 nca = right; |
|
262 while (left_set.find(nca) == left_set.end()) { |
|
263 nca = _graph.target((*_matching)[nca]); |
|
264 nca =(*_blossom_rep)[_blossom_set-> |
|
265 find(_graph.target((*_ear)[nca]))]; |
|
266 } |
|
267 } else { |
|
268 nca = left; |
|
269 while (right_set.find(nca) == right_set.end()) { |
|
270 nca = _graph.target((*_matching)[nca]); |
|
271 nca = (*_blossom_rep)[_blossom_set-> |
|
272 find(_graph.target((*_ear)[nca]))]; |
|
273 } |
|
274 } |
|
275 } |
|
276 } |
|
277 |
|
278 { |
|
279 |
|
280 Node node = _graph.u(e); |
|
281 Arc arc = _graph.direct(e, true); |
|
282 Node base = (*_blossom_rep)[_blossom_set->find(node)]; |
|
283 |
|
284 while (base != nca) { |
|
285 (*_ear)[node] = arc; |
|
286 |
|
287 Node n = node; |
|
288 while (n != base) { |
|
289 n = _graph.target((*_matching)[n]); |
|
290 Arc a = (*_ear)[n]; |
|
291 n = _graph.target(a); |
|
292 (*_ear)[n] = _graph.oppositeArc(a); |
|
293 } |
|
294 node = _graph.target((*_matching)[base]); |
|
295 _tree_set->erase(base); |
|
296 _tree_set->erase(node); |
|
297 _blossom_set->insert(node, _blossom_set->find(base)); |
|
298 (*_status)[node] = EVEN; |
|
299 _node_queue[_last++] = node; |
|
300 arc = _graph.oppositeArc((*_ear)[node]); |
|
301 node = _graph.target((*_ear)[node]); |
|
302 base = (*_blossom_rep)[_blossom_set->find(node)]; |
|
303 _blossom_set->join(_graph.target(arc), base); |
|
304 } |
|
305 } |
|
306 |
|
307 (*_blossom_rep)[_blossom_set->find(nca)] = nca; |
|
308 |
|
309 { |
|
310 |
|
311 Node node = _graph.v(e); |
|
312 Arc arc = _graph.direct(e, false); |
|
313 Node base = (*_blossom_rep)[_blossom_set->find(node)]; |
|
314 |
|
315 while (base != nca) { |
|
316 (*_ear)[node] = arc; |
|
317 |
|
318 Node n = node; |
|
319 while (n != base) { |
|
320 n = _graph.target((*_matching)[n]); |
|
321 Arc a = (*_ear)[n]; |
|
322 n = _graph.target(a); |
|
323 (*_ear)[n] = _graph.oppositeArc(a); |
|
324 } |
|
325 node = _graph.target((*_matching)[base]); |
|
326 _tree_set->erase(base); |
|
327 _tree_set->erase(node); |
|
328 _blossom_set->insert(node, _blossom_set->find(base)); |
|
329 (*_status)[node] = EVEN; |
|
330 _node_queue[_last++] = node; |
|
331 arc = _graph.oppositeArc((*_ear)[node]); |
|
332 node = _graph.target((*_ear)[node]); |
|
333 base = (*_blossom_rep)[_blossom_set->find(node)]; |
|
334 _blossom_set->join(_graph.target(arc), base); |
|
335 } |
|
336 } |
|
337 |
|
338 (*_blossom_rep)[_blossom_set->find(nca)] = nca; |
|
339 } |
|
340 |
|
341 |
|
342 |
|
343 void extendOnArc(const Arc& a) { |
|
344 Node base = _graph.source(a); |
|
345 Node odd = _graph.target(a); |
|
346 |
|
347 (*_ear)[odd] = _graph.oppositeArc(a); |
|
348 Node even = _graph.target((*_matching)[odd]); |
|
349 (*_blossom_rep)[_blossom_set->insert(even)] = even; |
|
350 (*_status)[odd] = ODD; |
|
351 (*_status)[even] = EVEN; |
|
352 int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]); |
|
353 _tree_set->insert(odd, tree); |
|
354 _tree_set->insert(even, tree); |
|
355 _node_queue[_last++] = even; |
|
356 |
|
357 } |
|
358 |
|
359 void augmentOnArc(const Arc& a) { |
|
360 Node even = _graph.source(a); |
|
361 Node odd = _graph.target(a); |
|
362 |
|
363 int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]); |
|
364 |
|
365 (*_matching)[odd] = _graph.oppositeArc(a); |
|
366 (*_status)[odd] = MATCHED; |
|
367 |
|
368 Arc arc = (*_matching)[even]; |
|
369 (*_matching)[even] = a; |
|
370 |
|
371 while (arc != INVALID) { |
|
372 odd = _graph.target(arc); |
|
373 arc = (*_ear)[odd]; |
|
374 even = _graph.target(arc); |
|
375 (*_matching)[odd] = arc; |
|
376 arc = (*_matching)[even]; |
|
377 (*_matching)[even] = _graph.oppositeArc((*_matching)[odd]); |
|
378 } |
|
379 |
|
380 for (typename TreeSet::ItemIt it(*_tree_set, tree); |
|
381 it != INVALID; ++it) { |
|
382 if ((*_status)[it] == ODD) { |
|
383 (*_status)[it] = MATCHED; |
|
384 } else { |
|
385 int blossom = _blossom_set->find(it); |
|
386 for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom); |
|
387 jt != INVALID; ++jt) { |
|
388 (*_status)[jt] = MATCHED; |
|
389 } |
|
390 _blossom_set->eraseClass(blossom); |
|
391 } |
|
392 } |
|
393 _tree_set->eraseClass(tree); |
|
394 |
|
395 } |
|
396 |
|
397 public: |
|
398 |
|
399 /// \brief Constructor |
|
400 /// |
|
401 /// Constructor. |
|
402 MaxMatching(const Graph& graph) |
|
403 : _graph(graph), _matching(0), _status(0), _ear(0), |
|
404 _blossom_set_index(0), _blossom_set(0), _blossom_rep(0), |
|
405 _tree_set_index(0), _tree_set(0) {} |
|
406 |
|
407 ~MaxMatching() { |
|
408 destroyStructures(); |
|
409 } |
|
410 |
|
411 /// \name Execution control |
|
412 /// The simplest way to execute the algorithm is to use the |
|
413 /// \c run() member function. |
|
414 /// \n |
|
415 |
|
416 /// If you need better control on the execution, you must call |
|
417 /// \ref init(), \ref greedyInit() or \ref matchingInit() |
|
418 /// functions first, then you can start the algorithm with the \ref |
|
419 /// startSparse() or startDense() functions. |
|
420 |
|
421 ///@{ |
|
422 |
|
423 /// \brief Sets the actual matching to the empty matching. |
|
424 /// |
|
425 /// Sets the actual matching to the empty matching. |
|
426 /// |
|
427 void init() { |
|
428 createStructures(); |
|
429 for(NodeIt n(_graph); n != INVALID; ++n) { |
|
430 (*_matching)[n] = INVALID; |
|
431 (*_status)[n] = UNMATCHED; |
|
432 } |
|
433 } |
|
434 |
|
435 ///\brief Finds an initial matching in a greedy way |
|
436 /// |
|
437 ///It finds an initial matching in a greedy way. |
|
438 void greedyInit() { |
|
439 createStructures(); |
|
440 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
441 (*_matching)[n] = INVALID; |
|
442 (*_status)[n] = UNMATCHED; |
|
443 } |
|
444 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
445 if ((*_matching)[n] == INVALID) { |
|
446 for (OutArcIt a(_graph, n); a != INVALID ; ++a) { |
|
447 Node v = _graph.target(a); |
|
448 if ((*_matching)[v] == INVALID && v != n) { |
|
449 (*_matching)[n] = a; |
|
450 (*_status)[n] = MATCHED; |
|
451 (*_matching)[v] = _graph.oppositeArc(a); |
|
452 (*_status)[v] = MATCHED; |
|
453 break; |
|
454 } |
|
455 } |
|
456 } |
|
457 } |
|
458 } |
|
459 |
|
460 |
|
461 /// \brief Initialize the matching from a map containing. |
|
462 /// |
|
463 /// Initialize the matching from a \c bool valued \c Edge map. This |
|
464 /// map must have the property that there are no two incident edges |
|
465 /// with true value, ie. it contains a matching. |
|
466 /// \return \c true if the map contains a matching. |
|
467 template <typename MatchingMap> |
|
468 bool matchingInit(const MatchingMap& matching) { |
|
469 createStructures(); |
|
470 |
|
471 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
472 (*_matching)[n] = INVALID; |
|
473 (*_status)[n] = UNMATCHED; |
|
474 } |
|
475 for(EdgeIt e(_graph); e!=INVALID; ++e) { |
|
476 if (matching[e]) { |
|
477 |
|
478 Node u = _graph.u(e); |
|
479 if ((*_matching)[u] != INVALID) return false; |
|
480 (*_matching)[u] = _graph.direct(e, true); |
|
481 (*_status)[u] = MATCHED; |
|
482 |
|
483 Node v = _graph.v(e); |
|
484 if ((*_matching)[v] != INVALID) return false; |
|
485 (*_matching)[v] = _graph.direct(e, false); |
|
486 (*_status)[v] = MATCHED; |
|
487 } |
|
488 } |
|
489 return true; |
|
490 } |
|
491 |
|
492 /// \brief Starts Edmonds' algorithm |
|
493 /// |
|
494 /// If runs the original Edmonds' algorithm. |
|
495 void startSparse() { |
|
496 for(NodeIt n(_graph); n != INVALID; ++n) { |
|
497 if ((*_status)[n] == UNMATCHED) { |
|
498 (*_blossom_rep)[_blossom_set->insert(n)] = n; |
|
499 _tree_set->insert(n); |
|
500 (*_status)[n] = EVEN; |
|
501 processSparse(n); |
|
502 } |
|
503 } |
|
504 } |
|
505 |
|
506 /// \brief Starts Edmonds' algorithm. |
|
507 /// |
|
508 /// It runs Edmonds' algorithm with a heuristic of postponing |
|
509 /// shrinks, therefore resulting in a faster algorithm for dense graphs. |
|
510 void startDense() { |
|
511 for(NodeIt n(_graph); n != INVALID; ++n) { |
|
512 if ((*_status)[n] == UNMATCHED) { |
|
513 (*_blossom_rep)[_blossom_set->insert(n)] = n; |
|
514 _tree_set->insert(n); |
|
515 (*_status)[n] = EVEN; |
|
516 processDense(n); |
|
517 } |
|
518 } |
|
519 } |
|
520 |
|
521 |
|
522 /// \brief Runs Edmonds' algorithm |
|
523 /// |
|
524 /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>) |
|
525 /// or Edmonds' algorithm with a heuristic of |
|
526 /// postponing shrinks for dense graphs. |
|
527 void run() { |
|
528 if (countEdges(_graph) < 2 * countNodes(_graph)) { |
|
529 greedyInit(); |
|
530 startSparse(); |
|
531 } else { |
|
532 init(); |
|
533 startDense(); |
|
534 } |
|
535 } |
|
536 |
|
537 /// @} |
|
538 |
|
539 /// \name Primal solution |
|
540 /// Functions to get the primal solution, ie. the matching. |
|
541 |
|
542 /// @{ |
|
543 |
|
544 ///\brief Returns the size of the current matching. |
|
545 /// |
|
546 ///Returns the size of the current matching. After \ref |
|
547 ///run() it returns the size of the maximum matching in the graph. |
|
548 int matchingSize() const { |
|
549 int size = 0; |
|
550 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
551 if ((*_matching)[n] != INVALID) { |
|
552 ++size; |
|
553 } |
|
554 } |
|
555 return size / 2; |
|
556 } |
|
557 |
|
558 /// \brief Returns true when the edge is in the matching. |
|
559 /// |
|
560 /// Returns true when the edge is in the matching. |
|
561 bool matching(const Edge& edge) const { |
|
562 return edge == (*_matching)[_graph.u(edge)]; |
|
563 } |
|
564 |
|
565 /// \brief Returns the matching edge incident to the given node. |
|
566 /// |
|
567 /// Returns the matching edge of a \c node in the actual matching or |
|
568 /// INVALID if the \c node is not covered by the actual matching. |
|
569 Arc matching(const Node& n) const { |
|
570 return (*_matching)[n]; |
|
571 } |
|
572 |
|
573 ///\brief Returns the mate of a node in the actual matching. |
|
574 /// |
|
575 ///Returns the mate of a \c node in the actual matching or |
|
576 ///INVALID if the \c node is not covered by the actual matching. |
|
577 Node mate(const Node& n) const { |
|
578 return (*_matching)[n] != INVALID ? |
|
579 _graph.target((*_matching)[n]) : INVALID; |
|
580 } |
|
581 |
|
582 /// @} |
|
583 |
|
584 /// \name Dual solution |
|
585 /// Functions to get the dual solution, ie. the decomposition. |
|
586 |
|
587 /// @{ |
|
588 |
|
589 /// \brief Returns the class of the node in the Edmonds-Gallai |
|
590 /// decomposition. |
|
591 /// |
|
592 /// Returns the class of the node in the Edmonds-Gallai |
|
593 /// decomposition. |
|
594 Status decomposition(const Node& n) const { |
|
595 return (*_status)[n]; |
|
596 } |
|
597 |
|
598 /// \brief Returns true when the node is in the barrier. |
|
599 /// |
|
600 /// Returns true when the node is in the barrier. |
|
601 bool barrier(const Node& n) const { |
|
602 return (*_status)[n] == ODD; |
|
603 } |
|
604 |
|
605 /// @} |
|
606 |
|
607 }; |
|
608 |
|
609 /// \ingroup matching |
|
610 /// |
|
611 /// \brief Weighted matching in general graphs |
|
612 /// |
|
613 /// This class provides an efficient implementation of Edmond's |
|
614 /// maximum weighted matching algorithm. The implementation is based |
|
615 /// on extensive use of priority queues and provides |
|
616 /// \f$O(nm\log n)\f$ time complexity. |
|
617 /// |
|
618 /// The maximum weighted matching problem is to find undirected |
|
619 /// edges in the graph with maximum overall weight and no two of |
|
620 /// them shares their ends. The problem can be formulated with the |
|
621 /// following linear program. |
|
622 /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f] |
|
623 /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} |
|
624 \quad \forall B\in\mathcal{O}\f] */ |
|
625 /// \f[x_e \ge 0\quad \forall e\in E\f] |
|
626 /// \f[\max \sum_{e\in E}x_ew_e\f] |
|
627 /// where \f$\delta(X)\f$ is the set of edges incident to a node in |
|
628 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in |
|
629 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality |
|
630 /// subsets of the nodes. |
|
631 /// |
|
632 /// The algorithm calculates an optimal matching and a proof of the |
|
633 /// optimality. The solution of the dual problem can be used to check |
|
634 /// the result of the algorithm. The dual linear problem is the |
|
635 /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)} |
|
636 z_B \ge w_{uv} \quad \forall uv\in E\f] */ |
|
637 /// \f[y_u \ge 0 \quad \forall u \in V\f] |
|
638 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] |
|
639 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} |
|
640 \frac{\vert B \vert - 1}{2}z_B\f] */ |
|
641 /// |
|
642 /// The algorithm can be executed with \c run() or the \c init() and |
|
643 /// then the \c start() member functions. After it the matching can |
|
644 /// be asked with \c matching() or mate() functions. The dual |
|
645 /// solution can be get with \c nodeValue(), \c blossomNum() and \c |
|
646 /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt |
|
647 /// "BlossomIt" nested class, which is able to iterate on the nodes |
|
648 /// of a blossom. If the value type is integral then the dual |
|
649 /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4". |
|
650 template <typename GR, |
|
651 typename WM = typename GR::template EdgeMap<int> > |
|
652 class MaxWeightedMatching { |
|
653 public: |
|
654 |
|
655 ///\e |
|
656 typedef GR Graph; |
|
657 ///\e |
|
658 typedef WM WeightMap; |
|
659 ///\e |
|
660 typedef typename WeightMap::Value Value; |
|
661 |
|
662 /// \brief Scaling factor for dual solution |
|
663 /// |
|
664 /// Scaling factor for dual solution, it is equal to 4 or 1 |
|
665 /// according to the value type. |
|
666 static const int dualScale = |
|
667 std::numeric_limits<Value>::is_integer ? 4 : 1; |
|
668 |
|
669 typedef typename Graph::template NodeMap<typename Graph::Arc> |
|
670 MatchingMap; |
|
671 |
|
672 private: |
|
673 |
|
674 TEMPLATE_GRAPH_TYPEDEFS(Graph); |
|
675 |
|
676 typedef typename Graph::template NodeMap<Value> NodePotential; |
|
677 typedef std::vector<Node> BlossomNodeList; |
|
678 |
|
679 struct BlossomVariable { |
|
680 int begin, end; |
|
681 Value value; |
|
682 |
|
683 BlossomVariable(int _begin, int _end, Value _value) |
|
684 : begin(_begin), end(_end), value(_value) {} |
|
685 |
|
686 }; |
|
687 |
|
688 typedef std::vector<BlossomVariable> BlossomPotential; |
|
689 |
|
690 const Graph& _graph; |
|
691 const WeightMap& _weight; |
|
692 |
|
693 MatchingMap* _matching; |
|
694 |
|
695 NodePotential* _node_potential; |
|
696 |
|
697 BlossomPotential _blossom_potential; |
|
698 BlossomNodeList _blossom_node_list; |
|
699 |
|
700 int _node_num; |
|
701 int _blossom_num; |
|
702 |
|
703 typedef RangeMap<int> IntIntMap; |
|
704 |
|
705 enum Status { |
|
706 EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2 |
|
707 }; |
|
708 |
|
709 typedef HeapUnionFind<Value, IntNodeMap> BlossomSet; |
|
710 struct BlossomData { |
|
711 int tree; |
|
712 Status status; |
|
713 Arc pred, next; |
|
714 Value pot, offset; |
|
715 Node base; |
|
716 }; |
|
717 |
|
718 IntNodeMap *_blossom_index; |
|
719 BlossomSet *_blossom_set; |
|
720 RangeMap<BlossomData>* _blossom_data; |
|
721 |
|
722 IntNodeMap *_node_index; |
|
723 IntArcMap *_node_heap_index; |
|
724 |
|
725 struct NodeData { |
|
726 |
|
727 NodeData(IntArcMap& node_heap_index) |
|
728 : heap(node_heap_index) {} |
|
729 |
|
730 int blossom; |
|
731 Value pot; |
|
732 BinHeap<Value, IntArcMap> heap; |
|
733 std::map<int, Arc> heap_index; |
|
734 |
|
735 int tree; |
|
736 }; |
|
737 |
|
738 RangeMap<NodeData>* _node_data; |
|
739 |
|
740 typedef ExtendFindEnum<IntIntMap> TreeSet; |
|
741 |
|
742 IntIntMap *_tree_set_index; |
|
743 TreeSet *_tree_set; |
|
744 |
|
745 IntNodeMap *_delta1_index; |
|
746 BinHeap<Value, IntNodeMap> *_delta1; |
|
747 |
|
748 IntIntMap *_delta2_index; |
|
749 BinHeap<Value, IntIntMap> *_delta2; |
|
750 |
|
751 IntEdgeMap *_delta3_index; |
|
752 BinHeap<Value, IntEdgeMap> *_delta3; |
|
753 |
|
754 IntIntMap *_delta4_index; |
|
755 BinHeap<Value, IntIntMap> *_delta4; |
|
756 |
|
757 Value _delta_sum; |
|
758 |
|
759 void createStructures() { |
|
760 _node_num = countNodes(_graph); |
|
761 _blossom_num = _node_num * 3 / 2; |
|
762 |
|
763 if (!_matching) { |
|
764 _matching = new MatchingMap(_graph); |
|
765 } |
|
766 if (!_node_potential) { |
|
767 _node_potential = new NodePotential(_graph); |
|
768 } |
|
769 if (!_blossom_set) { |
|
770 _blossom_index = new IntNodeMap(_graph); |
|
771 _blossom_set = new BlossomSet(*_blossom_index); |
|
772 _blossom_data = new RangeMap<BlossomData>(_blossom_num); |
|
773 } |
|
774 |
|
775 if (!_node_index) { |
|
776 _node_index = new IntNodeMap(_graph); |
|
777 _node_heap_index = new IntArcMap(_graph); |
|
778 _node_data = new RangeMap<NodeData>(_node_num, |
|
779 NodeData(*_node_heap_index)); |
|
780 } |
|
781 |
|
782 if (!_tree_set) { |
|
783 _tree_set_index = new IntIntMap(_blossom_num); |
|
784 _tree_set = new TreeSet(*_tree_set_index); |
|
785 } |
|
786 if (!_delta1) { |
|
787 _delta1_index = new IntNodeMap(_graph); |
|
788 _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index); |
|
789 } |
|
790 if (!_delta2) { |
|
791 _delta2_index = new IntIntMap(_blossom_num); |
|
792 _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index); |
|
793 } |
|
794 if (!_delta3) { |
|
795 _delta3_index = new IntEdgeMap(_graph); |
|
796 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index); |
|
797 } |
|
798 if (!_delta4) { |
|
799 _delta4_index = new IntIntMap(_blossom_num); |
|
800 _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index); |
|
801 } |
|
802 } |
|
803 |
|
804 void destroyStructures() { |
|
805 _node_num = countNodes(_graph); |
|
806 _blossom_num = _node_num * 3 / 2; |
|
807 |
|
808 if (_matching) { |
|
809 delete _matching; |
|
810 } |
|
811 if (_node_potential) { |
|
812 delete _node_potential; |
|
813 } |
|
814 if (_blossom_set) { |
|
815 delete _blossom_index; |
|
816 delete _blossom_set; |
|
817 delete _blossom_data; |
|
818 } |
|
819 |
|
820 if (_node_index) { |
|
821 delete _node_index; |
|
822 delete _node_heap_index; |
|
823 delete _node_data; |
|
824 } |
|
825 |
|
826 if (_tree_set) { |
|
827 delete _tree_set_index; |
|
828 delete _tree_set; |
|
829 } |
|
830 if (_delta1) { |
|
831 delete _delta1_index; |
|
832 delete _delta1; |
|
833 } |
|
834 if (_delta2) { |
|
835 delete _delta2_index; |
|
836 delete _delta2; |
|
837 } |
|
838 if (_delta3) { |
|
839 delete _delta3_index; |
|
840 delete _delta3; |
|
841 } |
|
842 if (_delta4) { |
|
843 delete _delta4_index; |
|
844 delete _delta4; |
|
845 } |
|
846 } |
|
847 |
|
848 void matchedToEven(int blossom, int tree) { |
|
849 if (_delta2->state(blossom) == _delta2->IN_HEAP) { |
|
850 _delta2->erase(blossom); |
|
851 } |
|
852 |
|
853 if (!_blossom_set->trivial(blossom)) { |
|
854 (*_blossom_data)[blossom].pot -= |
|
855 2 * (_delta_sum - (*_blossom_data)[blossom].offset); |
|
856 } |
|
857 |
|
858 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
|
859 n != INVALID; ++n) { |
|
860 |
|
861 _blossom_set->increase(n, std::numeric_limits<Value>::max()); |
|
862 int ni = (*_node_index)[n]; |
|
863 |
|
864 (*_node_data)[ni].heap.clear(); |
|
865 (*_node_data)[ni].heap_index.clear(); |
|
866 |
|
867 (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset; |
|
868 |
|
869 _delta1->push(n, (*_node_data)[ni].pot); |
|
870 |
|
871 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
872 Node v = _graph.source(e); |
|
873 int vb = _blossom_set->find(v); |
|
874 int vi = (*_node_index)[v]; |
|
875 |
|
876 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - |
|
877 dualScale * _weight[e]; |
|
878 |
|
879 if ((*_blossom_data)[vb].status == EVEN) { |
|
880 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { |
|
881 _delta3->push(e, rw / 2); |
|
882 } |
|
883 } else if ((*_blossom_data)[vb].status == UNMATCHED) { |
|
884 if (_delta3->state(e) != _delta3->IN_HEAP) { |
|
885 _delta3->push(e, rw); |
|
886 } |
|
887 } else { |
|
888 typename std::map<int, Arc>::iterator it = |
|
889 (*_node_data)[vi].heap_index.find(tree); |
|
890 |
|
891 if (it != (*_node_data)[vi].heap_index.end()) { |
|
892 if ((*_node_data)[vi].heap[it->second] > rw) { |
|
893 (*_node_data)[vi].heap.replace(it->second, e); |
|
894 (*_node_data)[vi].heap.decrease(e, rw); |
|
895 it->second = e; |
|
896 } |
|
897 } else { |
|
898 (*_node_data)[vi].heap.push(e, rw); |
|
899 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); |
|
900 } |
|
901 |
|
902 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { |
|
903 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); |
|
904 |
|
905 if ((*_blossom_data)[vb].status == MATCHED) { |
|
906 if (_delta2->state(vb) != _delta2->IN_HEAP) { |
|
907 _delta2->push(vb, _blossom_set->classPrio(vb) - |
|
908 (*_blossom_data)[vb].offset); |
|
909 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - |
|
910 (*_blossom_data)[vb].offset){ |
|
911 _delta2->decrease(vb, _blossom_set->classPrio(vb) - |
|
912 (*_blossom_data)[vb].offset); |
|
913 } |
|
914 } |
|
915 } |
|
916 } |
|
917 } |
|
918 } |
|
919 (*_blossom_data)[blossom].offset = 0; |
|
920 } |
|
921 |
|
922 void matchedToOdd(int blossom) { |
|
923 if (_delta2->state(blossom) == _delta2->IN_HEAP) { |
|
924 _delta2->erase(blossom); |
|
925 } |
|
926 (*_blossom_data)[blossom].offset += _delta_sum; |
|
927 if (!_blossom_set->trivial(blossom)) { |
|
928 _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + |
|
929 (*_blossom_data)[blossom].offset); |
|
930 } |
|
931 } |
|
932 |
|
933 void evenToMatched(int blossom, int tree) { |
|
934 if (!_blossom_set->trivial(blossom)) { |
|
935 (*_blossom_data)[blossom].pot += 2 * _delta_sum; |
|
936 } |
|
937 |
|
938 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
|
939 n != INVALID; ++n) { |
|
940 int ni = (*_node_index)[n]; |
|
941 (*_node_data)[ni].pot -= _delta_sum; |
|
942 |
|
943 _delta1->erase(n); |
|
944 |
|
945 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
946 Node v = _graph.source(e); |
|
947 int vb = _blossom_set->find(v); |
|
948 int vi = (*_node_index)[v]; |
|
949 |
|
950 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - |
|
951 dualScale * _weight[e]; |
|
952 |
|
953 if (vb == blossom) { |
|
954 if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
955 _delta3->erase(e); |
|
956 } |
|
957 } else if ((*_blossom_data)[vb].status == EVEN) { |
|
958 |
|
959 if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
960 _delta3->erase(e); |
|
961 } |
|
962 |
|
963 int vt = _tree_set->find(vb); |
|
964 |
|
965 if (vt != tree) { |
|
966 |
|
967 Arc r = _graph.oppositeArc(e); |
|
968 |
|
969 typename std::map<int, Arc>::iterator it = |
|
970 (*_node_data)[ni].heap_index.find(vt); |
|
971 |
|
972 if (it != (*_node_data)[ni].heap_index.end()) { |
|
973 if ((*_node_data)[ni].heap[it->second] > rw) { |
|
974 (*_node_data)[ni].heap.replace(it->second, r); |
|
975 (*_node_data)[ni].heap.decrease(r, rw); |
|
976 it->second = r; |
|
977 } |
|
978 } else { |
|
979 (*_node_data)[ni].heap.push(r, rw); |
|
980 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); |
|
981 } |
|
982 |
|
983 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { |
|
984 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); |
|
985 |
|
986 if (_delta2->state(blossom) != _delta2->IN_HEAP) { |
|
987 _delta2->push(blossom, _blossom_set->classPrio(blossom) - |
|
988 (*_blossom_data)[blossom].offset); |
|
989 } else if ((*_delta2)[blossom] > |
|
990 _blossom_set->classPrio(blossom) - |
|
991 (*_blossom_data)[blossom].offset){ |
|
992 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - |
|
993 (*_blossom_data)[blossom].offset); |
|
994 } |
|
995 } |
|
996 } |
|
997 |
|
998 } else if ((*_blossom_data)[vb].status == UNMATCHED) { |
|
999 if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
1000 _delta3->erase(e); |
|
1001 } |
|
1002 } else { |
|
1003 |
|
1004 typename std::map<int, Arc>::iterator it = |
|
1005 (*_node_data)[vi].heap_index.find(tree); |
|
1006 |
|
1007 if (it != (*_node_data)[vi].heap_index.end()) { |
|
1008 (*_node_data)[vi].heap.erase(it->second); |
|
1009 (*_node_data)[vi].heap_index.erase(it); |
|
1010 if ((*_node_data)[vi].heap.empty()) { |
|
1011 _blossom_set->increase(v, std::numeric_limits<Value>::max()); |
|
1012 } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) { |
|
1013 _blossom_set->increase(v, (*_node_data)[vi].heap.prio()); |
|
1014 } |
|
1015 |
|
1016 if ((*_blossom_data)[vb].status == MATCHED) { |
|
1017 if (_blossom_set->classPrio(vb) == |
|
1018 std::numeric_limits<Value>::max()) { |
|
1019 _delta2->erase(vb); |
|
1020 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) - |
|
1021 (*_blossom_data)[vb].offset) { |
|
1022 _delta2->increase(vb, _blossom_set->classPrio(vb) - |
|
1023 (*_blossom_data)[vb].offset); |
|
1024 } |
|
1025 } |
|
1026 } |
|
1027 } |
|
1028 } |
|
1029 } |
|
1030 } |
|
1031 |
|
1032 void oddToMatched(int blossom) { |
|
1033 (*_blossom_data)[blossom].offset -= _delta_sum; |
|
1034 |
|
1035 if (_blossom_set->classPrio(blossom) != |
|
1036 std::numeric_limits<Value>::max()) { |
|
1037 _delta2->push(blossom, _blossom_set->classPrio(blossom) - |
|
1038 (*_blossom_data)[blossom].offset); |
|
1039 } |
|
1040 |
|
1041 if (!_blossom_set->trivial(blossom)) { |
|
1042 _delta4->erase(blossom); |
|
1043 } |
|
1044 } |
|
1045 |
|
1046 void oddToEven(int blossom, int tree) { |
|
1047 if (!_blossom_set->trivial(blossom)) { |
|
1048 _delta4->erase(blossom); |
|
1049 (*_blossom_data)[blossom].pot -= |
|
1050 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset); |
|
1051 } |
|
1052 |
|
1053 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
|
1054 n != INVALID; ++n) { |
|
1055 int ni = (*_node_index)[n]; |
|
1056 |
|
1057 _blossom_set->increase(n, std::numeric_limits<Value>::max()); |
|
1058 |
|
1059 (*_node_data)[ni].heap.clear(); |
|
1060 (*_node_data)[ni].heap_index.clear(); |
|
1061 (*_node_data)[ni].pot += |
|
1062 2 * _delta_sum - (*_blossom_data)[blossom].offset; |
|
1063 |
|
1064 _delta1->push(n, (*_node_data)[ni].pot); |
|
1065 |
|
1066 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
1067 Node v = _graph.source(e); |
|
1068 int vb = _blossom_set->find(v); |
|
1069 int vi = (*_node_index)[v]; |
|
1070 |
|
1071 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - |
|
1072 dualScale * _weight[e]; |
|
1073 |
|
1074 if ((*_blossom_data)[vb].status == EVEN) { |
|
1075 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) { |
|
1076 _delta3->push(e, rw / 2); |
|
1077 } |
|
1078 } else if ((*_blossom_data)[vb].status == UNMATCHED) { |
|
1079 if (_delta3->state(e) != _delta3->IN_HEAP) { |
|
1080 _delta3->push(e, rw); |
|
1081 } |
|
1082 } else { |
|
1083 |
|
1084 typename std::map<int, Arc>::iterator it = |
|
1085 (*_node_data)[vi].heap_index.find(tree); |
|
1086 |
|
1087 if (it != (*_node_data)[vi].heap_index.end()) { |
|
1088 if ((*_node_data)[vi].heap[it->second] > rw) { |
|
1089 (*_node_data)[vi].heap.replace(it->second, e); |
|
1090 (*_node_data)[vi].heap.decrease(e, rw); |
|
1091 it->second = e; |
|
1092 } |
|
1093 } else { |
|
1094 (*_node_data)[vi].heap.push(e, rw); |
|
1095 (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e)); |
|
1096 } |
|
1097 |
|
1098 if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) { |
|
1099 _blossom_set->decrease(v, (*_node_data)[vi].heap.prio()); |
|
1100 |
|
1101 if ((*_blossom_data)[vb].status == MATCHED) { |
|
1102 if (_delta2->state(vb) != _delta2->IN_HEAP) { |
|
1103 _delta2->push(vb, _blossom_set->classPrio(vb) - |
|
1104 (*_blossom_data)[vb].offset); |
|
1105 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - |
|
1106 (*_blossom_data)[vb].offset) { |
|
1107 _delta2->decrease(vb, _blossom_set->classPrio(vb) - |
|
1108 (*_blossom_data)[vb].offset); |
|
1109 } |
|
1110 } |
|
1111 } |
|
1112 } |
|
1113 } |
|
1114 } |
|
1115 (*_blossom_data)[blossom].offset = 0; |
|
1116 } |
|
1117 |
|
1118 |
|
1119 void matchedToUnmatched(int blossom) { |
|
1120 if (_delta2->state(blossom) == _delta2->IN_HEAP) { |
|
1121 _delta2->erase(blossom); |
|
1122 } |
|
1123 |
|
1124 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
|
1125 n != INVALID; ++n) { |
|
1126 int ni = (*_node_index)[n]; |
|
1127 |
|
1128 _blossom_set->increase(n, std::numeric_limits<Value>::max()); |
|
1129 |
|
1130 (*_node_data)[ni].heap.clear(); |
|
1131 (*_node_data)[ni].heap_index.clear(); |
|
1132 |
|
1133 for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
|
1134 Node v = _graph.target(e); |
|
1135 int vb = _blossom_set->find(v); |
|
1136 int vi = (*_node_index)[v]; |
|
1137 |
|
1138 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - |
|
1139 dualScale * _weight[e]; |
|
1140 |
|
1141 if ((*_blossom_data)[vb].status == EVEN) { |
|
1142 if (_delta3->state(e) != _delta3->IN_HEAP) { |
|
1143 _delta3->push(e, rw); |
|
1144 } |
|
1145 } |
|
1146 } |
|
1147 } |
|
1148 } |
|
1149 |
|
1150 void unmatchedToMatched(int blossom) { |
|
1151 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); |
|
1152 n != INVALID; ++n) { |
|
1153 int ni = (*_node_index)[n]; |
|
1154 |
|
1155 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
1156 Node v = _graph.source(e); |
|
1157 int vb = _blossom_set->find(v); |
|
1158 int vi = (*_node_index)[v]; |
|
1159 |
|
1160 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - |
|
1161 dualScale * _weight[e]; |
|
1162 |
|
1163 if (vb == blossom) { |
|
1164 if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
1165 _delta3->erase(e); |
|
1166 } |
|
1167 } else if ((*_blossom_data)[vb].status == EVEN) { |
|
1168 |
|
1169 if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
1170 _delta3->erase(e); |
|
1171 } |
|
1172 |
|
1173 int vt = _tree_set->find(vb); |
|
1174 |
|
1175 Arc r = _graph.oppositeArc(e); |
|
1176 |
|
1177 typename std::map<int, Arc>::iterator it = |
|
1178 (*_node_data)[ni].heap_index.find(vt); |
|
1179 |
|
1180 if (it != (*_node_data)[ni].heap_index.end()) { |
|
1181 if ((*_node_data)[ni].heap[it->second] > rw) { |
|
1182 (*_node_data)[ni].heap.replace(it->second, r); |
|
1183 (*_node_data)[ni].heap.decrease(r, rw); |
|
1184 it->second = r; |
|
1185 } |
|
1186 } else { |
|
1187 (*_node_data)[ni].heap.push(r, rw); |
|
1188 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r)); |
|
1189 } |
|
1190 |
|
1191 if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) { |
|
1192 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio()); |
|
1193 |
|
1194 if (_delta2->state(blossom) != _delta2->IN_HEAP) { |
|
1195 _delta2->push(blossom, _blossom_set->classPrio(blossom) - |
|
1196 (*_blossom_data)[blossom].offset); |
|
1197 } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)- |
|
1198 (*_blossom_data)[blossom].offset){ |
|
1199 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) - |
|
1200 (*_blossom_data)[blossom].offset); |
|
1201 } |
|
1202 } |
|
1203 |
|
1204 } else if ((*_blossom_data)[vb].status == UNMATCHED) { |
|
1205 if (_delta3->state(e) == _delta3->IN_HEAP) { |
|
1206 _delta3->erase(e); |
|
1207 } |
|
1208 } |
|
1209 } |
|
1210 } |
|
1211 } |
|
1212 |
|
1213 void alternatePath(int even, int tree) { |
|
1214 int odd; |
|
1215 |
|
1216 evenToMatched(even, tree); |
|
1217 (*_blossom_data)[even].status = MATCHED; |
|
1218 |
|
1219 while ((*_blossom_data)[even].pred != INVALID) { |
|
1220 odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred)); |
|
1221 (*_blossom_data)[odd].status = MATCHED; |
|
1222 oddToMatched(odd); |
|
1223 (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred; |
|
1224 |
|
1225 even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred)); |
|
1226 (*_blossom_data)[even].status = MATCHED; |
|
1227 evenToMatched(even, tree); |
|
1228 (*_blossom_data)[even].next = |
|
1229 _graph.oppositeArc((*_blossom_data)[odd].pred); |
|
1230 } |
|
1231 |
|
1232 } |
|
1233 |
|
1234 void destroyTree(int tree) { |
|
1235 for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) { |
|
1236 if ((*_blossom_data)[b].status == EVEN) { |
|
1237 (*_blossom_data)[b].status = MATCHED; |
|
1238 evenToMatched(b, tree); |
|
1239 } else if ((*_blossom_data)[b].status == ODD) { |
|
1240 (*_blossom_data)[b].status = MATCHED; |
|
1241 oddToMatched(b); |
|
1242 } |
|
1243 } |
|
1244 _tree_set->eraseClass(tree); |
|
1245 } |
|
1246 |
|
1247 |
|
1248 void unmatchNode(const Node& node) { |
|
1249 int blossom = _blossom_set->find(node); |
|
1250 int tree = _tree_set->find(blossom); |
|
1251 |
|
1252 alternatePath(blossom, tree); |
|
1253 destroyTree(tree); |
|
1254 |
|
1255 (*_blossom_data)[blossom].status = UNMATCHED; |
|
1256 (*_blossom_data)[blossom].base = node; |
|
1257 matchedToUnmatched(blossom); |
|
1258 } |
|
1259 |
|
1260 |
|
1261 void augmentOnEdge(const Edge& edge) { |
|
1262 |
|
1263 int left = _blossom_set->find(_graph.u(edge)); |
|
1264 int right = _blossom_set->find(_graph.v(edge)); |
|
1265 |
|
1266 if ((*_blossom_data)[left].status == EVEN) { |
|
1267 int left_tree = _tree_set->find(left); |
|
1268 alternatePath(left, left_tree); |
|
1269 destroyTree(left_tree); |
|
1270 } else { |
|
1271 (*_blossom_data)[left].status = MATCHED; |
|
1272 unmatchedToMatched(left); |
|
1273 } |
|
1274 |
|
1275 if ((*_blossom_data)[right].status == EVEN) { |
|
1276 int right_tree = _tree_set->find(right); |
|
1277 alternatePath(right, right_tree); |
|
1278 destroyTree(right_tree); |
|
1279 } else { |
|
1280 (*_blossom_data)[right].status = MATCHED; |
|
1281 unmatchedToMatched(right); |
|
1282 } |
|
1283 |
|
1284 (*_blossom_data)[left].next = _graph.direct(edge, true); |
|
1285 (*_blossom_data)[right].next = _graph.direct(edge, false); |
|
1286 } |
|
1287 |
|
1288 void extendOnArc(const Arc& arc) { |
|
1289 int base = _blossom_set->find(_graph.target(arc)); |
|
1290 int tree = _tree_set->find(base); |
|
1291 |
|
1292 int odd = _blossom_set->find(_graph.source(arc)); |
|
1293 _tree_set->insert(odd, tree); |
|
1294 (*_blossom_data)[odd].status = ODD; |
|
1295 matchedToOdd(odd); |
|
1296 (*_blossom_data)[odd].pred = arc; |
|
1297 |
|
1298 int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next)); |
|
1299 (*_blossom_data)[even].pred = (*_blossom_data)[even].next; |
|
1300 _tree_set->insert(even, tree); |
|
1301 (*_blossom_data)[even].status = EVEN; |
|
1302 matchedToEven(even, tree); |
|
1303 } |
|
1304 |
|
1305 void shrinkOnEdge(const Edge& edge, int tree) { |
|
1306 int nca = -1; |
|
1307 std::vector<int> left_path, right_path; |
|
1308 |
|
1309 { |
|
1310 std::set<int> left_set, right_set; |
|
1311 int left = _blossom_set->find(_graph.u(edge)); |
|
1312 left_path.push_back(left); |
|
1313 left_set.insert(left); |
|
1314 |
|
1315 int right = _blossom_set->find(_graph.v(edge)); |
|
1316 right_path.push_back(right); |
|
1317 right_set.insert(right); |
|
1318 |
|
1319 while (true) { |
|
1320 |
|
1321 if ((*_blossom_data)[left].pred == INVALID) break; |
|
1322 |
|
1323 left = |
|
1324 _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); |
|
1325 left_path.push_back(left); |
|
1326 left = |
|
1327 _blossom_set->find(_graph.target((*_blossom_data)[left].pred)); |
|
1328 left_path.push_back(left); |
|
1329 |
|
1330 left_set.insert(left); |
|
1331 |
|
1332 if (right_set.find(left) != right_set.end()) { |
|
1333 nca = left; |
|
1334 break; |
|
1335 } |
|
1336 |
|
1337 if ((*_blossom_data)[right].pred == INVALID) break; |
|
1338 |
|
1339 right = |
|
1340 _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); |
|
1341 right_path.push_back(right); |
|
1342 right = |
|
1343 _blossom_set->find(_graph.target((*_blossom_data)[right].pred)); |
|
1344 right_path.push_back(right); |
|
1345 |
|
1346 right_set.insert(right); |
|
1347 |
|
1348 if (left_set.find(right) != left_set.end()) { |
|
1349 nca = right; |
|
1350 break; |
|
1351 } |
|
1352 |
|
1353 } |
|
1354 |
|
1355 if (nca == -1) { |
|
1356 if ((*_blossom_data)[left].pred == INVALID) { |
|
1357 nca = right; |
|
1358 while (left_set.find(nca) == left_set.end()) { |
|
1359 nca = |
|
1360 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); |
|
1361 right_path.push_back(nca); |
|
1362 nca = |
|
1363 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); |
|
1364 right_path.push_back(nca); |
|
1365 } |
|
1366 } else { |
|
1367 nca = left; |
|
1368 while (right_set.find(nca) == right_set.end()) { |
|
1369 nca = |
|
1370 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); |
|
1371 left_path.push_back(nca); |
|
1372 nca = |
|
1373 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred)); |
|
1374 left_path.push_back(nca); |
|
1375 } |
|
1376 } |
|
1377 } |
|
1378 } |
|
1379 |
|
1380 std::vector<int> subblossoms; |
|
1381 Arc prev; |
|
1382 |
|
1383 prev = _graph.direct(edge, true); |
|
1384 for (int i = 0; left_path[i] != nca; i += 2) { |
|
1385 subblossoms.push_back(left_path[i]); |
|
1386 (*_blossom_data)[left_path[i]].next = prev; |
|
1387 _tree_set->erase(left_path[i]); |
|
1388 |
|
1389 subblossoms.push_back(left_path[i + 1]); |
|
1390 (*_blossom_data)[left_path[i + 1]].status = EVEN; |
|
1391 oddToEven(left_path[i + 1], tree); |
|
1392 _tree_set->erase(left_path[i + 1]); |
|
1393 prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred); |
|
1394 } |
|
1395 |
|
1396 int k = 0; |
|
1397 while (right_path[k] != nca) ++k; |
|
1398 |
|
1399 subblossoms.push_back(nca); |
|
1400 (*_blossom_data)[nca].next = prev; |
|
1401 |
|
1402 for (int i = k - 2; i >= 0; i -= 2) { |
|
1403 subblossoms.push_back(right_path[i + 1]); |
|
1404 (*_blossom_data)[right_path[i + 1]].status = EVEN; |
|
1405 oddToEven(right_path[i + 1], tree); |
|
1406 _tree_set->erase(right_path[i + 1]); |
|
1407 |
|
1408 (*_blossom_data)[right_path[i + 1]].next = |
|
1409 (*_blossom_data)[right_path[i + 1]].pred; |
|
1410 |
|
1411 subblossoms.push_back(right_path[i]); |
|
1412 _tree_set->erase(right_path[i]); |
|
1413 } |
|
1414 |
|
1415 int surface = |
|
1416 _blossom_set->join(subblossoms.begin(), subblossoms.end()); |
|
1417 |
|
1418 for (int i = 0; i < int(subblossoms.size()); ++i) { |
|
1419 if (!_blossom_set->trivial(subblossoms[i])) { |
|
1420 (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum; |
|
1421 } |
|
1422 (*_blossom_data)[subblossoms[i]].status = MATCHED; |
|
1423 } |
|
1424 |
|
1425 (*_blossom_data)[surface].pot = -2 * _delta_sum; |
|
1426 (*_blossom_data)[surface].offset = 0; |
|
1427 (*_blossom_data)[surface].status = EVEN; |
|
1428 (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred; |
|
1429 (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred; |
|
1430 |
|
1431 _tree_set->insert(surface, tree); |
|
1432 _tree_set->erase(nca); |
|
1433 } |
|
1434 |
|
1435 void splitBlossom(int blossom) { |
|
1436 Arc next = (*_blossom_data)[blossom].next; |
|
1437 Arc pred = (*_blossom_data)[blossom].pred; |
|
1438 |
|
1439 int tree = _tree_set->find(blossom); |
|
1440 |
|
1441 (*_blossom_data)[blossom].status = MATCHED; |
|
1442 oddToMatched(blossom); |
|
1443 if (_delta2->state(blossom) == _delta2->IN_HEAP) { |
|
1444 _delta2->erase(blossom); |
|
1445 } |
|
1446 |
|
1447 std::vector<int> subblossoms; |
|
1448 _blossom_set->split(blossom, std::back_inserter(subblossoms)); |
|
1449 |
|
1450 Value offset = (*_blossom_data)[blossom].offset; |
|
1451 int b = _blossom_set->find(_graph.source(pred)); |
|
1452 int d = _blossom_set->find(_graph.source(next)); |
|
1453 |
|
1454 int ib = -1, id = -1; |
|
1455 for (int i = 0; i < int(subblossoms.size()); ++i) { |
|
1456 if (subblossoms[i] == b) ib = i; |
|
1457 if (subblossoms[i] == d) id = i; |
|
1458 |
|
1459 (*_blossom_data)[subblossoms[i]].offset = offset; |
|
1460 if (!_blossom_set->trivial(subblossoms[i])) { |
|
1461 (*_blossom_data)[subblossoms[i]].pot -= 2 * offset; |
|
1462 } |
|
1463 if (_blossom_set->classPrio(subblossoms[i]) != |
|
1464 std::numeric_limits<Value>::max()) { |
|
1465 _delta2->push(subblossoms[i], |
|
1466 _blossom_set->classPrio(subblossoms[i]) - |
|
1467 (*_blossom_data)[subblossoms[i]].offset); |
|
1468 } |
|
1469 } |
|
1470 |
|
1471 if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) { |
|
1472 for (int i = (id + 1) % subblossoms.size(); |
|
1473 i != ib; i = (i + 2) % subblossoms.size()) { |
|
1474 int sb = subblossoms[i]; |
|
1475 int tb = subblossoms[(i + 1) % subblossoms.size()]; |
|
1476 (*_blossom_data)[sb].next = |
|
1477 _graph.oppositeArc((*_blossom_data)[tb].next); |
|
1478 } |
|
1479 |
|
1480 for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) { |
|
1481 int sb = subblossoms[i]; |
|
1482 int tb = subblossoms[(i + 1) % subblossoms.size()]; |
|
1483 int ub = subblossoms[(i + 2) % subblossoms.size()]; |
|
1484 |
|
1485 (*_blossom_data)[sb].status = ODD; |
|
1486 matchedToOdd(sb); |
|
1487 _tree_set->insert(sb, tree); |
|
1488 (*_blossom_data)[sb].pred = pred; |
|
1489 (*_blossom_data)[sb].next = |
|
1490 _graph.oppositeArc((*_blossom_data)[tb].next); |
|
1491 |
|
1492 pred = (*_blossom_data)[ub].next; |
|
1493 |
|
1494 (*_blossom_data)[tb].status = EVEN; |
|
1495 matchedToEven(tb, tree); |
|
1496 _tree_set->insert(tb, tree); |
|
1497 (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next; |
|
1498 } |
|
1499 |
|
1500 (*_blossom_data)[subblossoms[id]].status = ODD; |
|
1501 matchedToOdd(subblossoms[id]); |
|
1502 _tree_set->insert(subblossoms[id], tree); |
|
1503 (*_blossom_data)[subblossoms[id]].next = next; |
|
1504 (*_blossom_data)[subblossoms[id]].pred = pred; |
|
1505 |
|
1506 } else { |
|
1507 |
|
1508 for (int i = (ib + 1) % subblossoms.size(); |
|
1509 i != id; i = (i + 2) % subblossoms.size()) { |
|
1510 int sb = subblossoms[i]; |
|
1511 int tb = subblossoms[(i + 1) % subblossoms.size()]; |
|
1512 (*_blossom_data)[sb].next = |
|
1513 _graph.oppositeArc((*_blossom_data)[tb].next); |
|
1514 } |
|
1515 |
|
1516 for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) { |
|
1517 int sb = subblossoms[i]; |
|
1518 int tb = subblossoms[(i + 1) % subblossoms.size()]; |
|
1519 int ub = subblossoms[(i + 2) % subblossoms.size()]; |
|
1520 |
|
1521 (*_blossom_data)[sb].status = ODD; |
|
1522 matchedToOdd(sb); |
|
1523 _tree_set->insert(sb, tree); |
|
1524 (*_blossom_data)[sb].next = next; |
|
1525 (*_blossom_data)[sb].pred = |
|
1526 _graph.oppositeArc((*_blossom_data)[tb].next); |
|
1527 |
|
1528 (*_blossom_data)[tb].status = EVEN; |
|
1529 matchedToEven(tb, tree); |
|
1530 _tree_set->insert(tb, tree); |
|
1531 (*_blossom_data)[tb].pred = |
|
1532 (*_blossom_data)[tb].next = |
|
1533 _graph.oppositeArc((*_blossom_data)[ub].next); |
|
1534 next = (*_blossom_data)[ub].next; |
|
1535 } |
|
1536 |
|
1537 (*_blossom_data)[subblossoms[ib]].status = ODD; |
|
1538 matchedToOdd(subblossoms[ib]); |
|
1539 _tree_set->insert(subblossoms[ib], tree); |
|
1540 (*_blossom_data)[subblossoms[ib]].next = next; |
|
1541 (*_blossom_data)[subblossoms[ib]].pred = pred; |
|
1542 } |
|
1543 _tree_set->erase(blossom); |
|
1544 } |
|
1545 |
|
1546 void extractBlossom(int blossom, const Node& base, const Arc& matching) { |
|
1547 if (_blossom_set->trivial(blossom)) { |
|
1548 int bi = (*_node_index)[base]; |
|
1549 Value pot = (*_node_data)[bi].pot; |
|
1550 |
|
1551 (*_matching)[base] = matching; |
|
1552 _blossom_node_list.push_back(base); |
|
1553 (*_node_potential)[base] = pot; |
|
1554 } else { |
|
1555 |
|
1556 Value pot = (*_blossom_data)[blossom].pot; |
|
1557 int bn = _blossom_node_list.size(); |
|
1558 |
|
1559 std::vector<int> subblossoms; |
|
1560 _blossom_set->split(blossom, std::back_inserter(subblossoms)); |
|
1561 int b = _blossom_set->find(base); |
|
1562 int ib = -1; |
|
1563 for (int i = 0; i < int(subblossoms.size()); ++i) { |
|
1564 if (subblossoms[i] == b) { ib = i; break; } |
|
1565 } |
|
1566 |
|
1567 for (int i = 1; i < int(subblossoms.size()); i += 2) { |
|
1568 int sb = subblossoms[(ib + i) % subblossoms.size()]; |
|
1569 int tb = subblossoms[(ib + i + 1) % subblossoms.size()]; |
|
1570 |
|
1571 Arc m = (*_blossom_data)[tb].next; |
|
1572 extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m)); |
|
1573 extractBlossom(tb, _graph.source(m), m); |
|
1574 } |
|
1575 extractBlossom(subblossoms[ib], base, matching); |
|
1576 |
|
1577 int en = _blossom_node_list.size(); |
|
1578 |
|
1579 _blossom_potential.push_back(BlossomVariable(bn, en, pot)); |
|
1580 } |
|
1581 } |
|
1582 |
|
1583 void extractMatching() { |
|
1584 std::vector<int> blossoms; |
|
1585 for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) { |
|
1586 blossoms.push_back(c); |
|
1587 } |
|
1588 |
|
1589 for (int i = 0; i < int(blossoms.size()); ++i) { |
|
1590 if ((*_blossom_data)[blossoms[i]].status == MATCHED) { |
|
1591 |
|
1592 Value offset = (*_blossom_data)[blossoms[i]].offset; |
|
1593 (*_blossom_data)[blossoms[i]].pot += 2 * offset; |
|
1594 for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); |
|
1595 n != INVALID; ++n) { |
|
1596 (*_node_data)[(*_node_index)[n]].pot -= offset; |
|
1597 } |
|
1598 |
|
1599 Arc matching = (*_blossom_data)[blossoms[i]].next; |
|
1600 Node base = _graph.source(matching); |
|
1601 extractBlossom(blossoms[i], base, matching); |
|
1602 } else { |
|
1603 Node base = (*_blossom_data)[blossoms[i]].base; |
|
1604 extractBlossom(blossoms[i], base, INVALID); |
|
1605 } |
|
1606 } |
|
1607 } |
|
1608 |
|
1609 public: |
|
1610 |
|
1611 /// \brief Constructor |
|
1612 /// |
|
1613 /// Constructor. |
|
1614 MaxWeightedMatching(const Graph& graph, const WeightMap& weight) |
|
1615 : _graph(graph), _weight(weight), _matching(0), |
|
1616 _node_potential(0), _blossom_potential(), _blossom_node_list(), |
|
1617 _node_num(0), _blossom_num(0), |
|
1618 |
|
1619 _blossom_index(0), _blossom_set(0), _blossom_data(0), |
|
1620 _node_index(0), _node_heap_index(0), _node_data(0), |
|
1621 _tree_set_index(0), _tree_set(0), |
|
1622 |
|
1623 _delta1_index(0), _delta1(0), |
|
1624 _delta2_index(0), _delta2(0), |
|
1625 _delta3_index(0), _delta3(0), |
|
1626 _delta4_index(0), _delta4(0), |
|
1627 |
|
1628 _delta_sum() {} |
|
1629 |
|
1630 ~MaxWeightedMatching() { |
|
1631 destroyStructures(); |
|
1632 } |
|
1633 |
|
1634 /// \name Execution control |
|
1635 /// The simplest way to execute the algorithm is to use the |
|
1636 /// \c run() member function. |
|
1637 |
|
1638 ///@{ |
|
1639 |
|
1640 /// \brief Initialize the algorithm |
|
1641 /// |
|
1642 /// Initialize the algorithm |
|
1643 void init() { |
|
1644 createStructures(); |
|
1645 |
|
1646 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
1647 (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP; |
|
1648 } |
|
1649 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
1650 (*_delta1_index)[n] = _delta1->PRE_HEAP; |
|
1651 } |
|
1652 for (EdgeIt e(_graph); e != INVALID; ++e) { |
|
1653 (*_delta3_index)[e] = _delta3->PRE_HEAP; |
|
1654 } |
|
1655 for (int i = 0; i < _blossom_num; ++i) { |
|
1656 (*_delta2_index)[i] = _delta2->PRE_HEAP; |
|
1657 (*_delta4_index)[i] = _delta4->PRE_HEAP; |
|
1658 } |
|
1659 |
|
1660 int index = 0; |
|
1661 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
1662 Value max = 0; |
|
1663 for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
|
1664 if (_graph.target(e) == n) continue; |
|
1665 if ((dualScale * _weight[e]) / 2 > max) { |
|
1666 max = (dualScale * _weight[e]) / 2; |
|
1667 } |
|
1668 } |
|
1669 (*_node_index)[n] = index; |
|
1670 (*_node_data)[index].pot = max; |
|
1671 _delta1->push(n, max); |
|
1672 int blossom = |
|
1673 _blossom_set->insert(n, std::numeric_limits<Value>::max()); |
|
1674 |
|
1675 _tree_set->insert(blossom); |
|
1676 |
|
1677 (*_blossom_data)[blossom].status = EVEN; |
|
1678 (*_blossom_data)[blossom].pred = INVALID; |
|
1679 (*_blossom_data)[blossom].next = INVALID; |
|
1680 (*_blossom_data)[blossom].pot = 0; |
|
1681 (*_blossom_data)[blossom].offset = 0; |
|
1682 ++index; |
|
1683 } |
|
1684 for (EdgeIt e(_graph); e != INVALID; ++e) { |
|
1685 int si = (*_node_index)[_graph.u(e)]; |
|
1686 int ti = (*_node_index)[_graph.v(e)]; |
|
1687 if (_graph.u(e) != _graph.v(e)) { |
|
1688 _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - |
|
1689 dualScale * _weight[e]) / 2); |
|
1690 } |
|
1691 } |
|
1692 } |
|
1693 |
|
1694 /// \brief Starts the algorithm |
|
1695 /// |
|
1696 /// Starts the algorithm |
|
1697 void start() { |
|
1698 enum OpType { |
|
1699 D1, D2, D3, D4 |
|
1700 }; |
|
1701 |
|
1702 int unmatched = _node_num; |
|
1703 while (unmatched > 0) { |
|
1704 Value d1 = !_delta1->empty() ? |
|
1705 _delta1->prio() : std::numeric_limits<Value>::max(); |
|
1706 |
|
1707 Value d2 = !_delta2->empty() ? |
|
1708 _delta2->prio() : std::numeric_limits<Value>::max(); |
|
1709 |
|
1710 Value d3 = !_delta3->empty() ? |
|
1711 _delta3->prio() : std::numeric_limits<Value>::max(); |
|
1712 |
|
1713 Value d4 = !_delta4->empty() ? |
|
1714 _delta4->prio() : std::numeric_limits<Value>::max(); |
|
1715 |
|
1716 _delta_sum = d1; OpType ot = D1; |
|
1717 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; } |
|
1718 if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; } |
|
1719 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; } |
|
1720 |
|
1721 |
|
1722 switch (ot) { |
|
1723 case D1: |
|
1724 { |
|
1725 Node n = _delta1->top(); |
|
1726 unmatchNode(n); |
|
1727 --unmatched; |
|
1728 } |
|
1729 break; |
|
1730 case D2: |
|
1731 { |
|
1732 int blossom = _delta2->top(); |
|
1733 Node n = _blossom_set->classTop(blossom); |
|
1734 Arc e = (*_node_data)[(*_node_index)[n]].heap.top(); |
|
1735 extendOnArc(e); |
|
1736 } |
|
1737 break; |
|
1738 case D3: |
|
1739 { |
|
1740 Edge e = _delta3->top(); |
|
1741 |
|
1742 int left_blossom = _blossom_set->find(_graph.u(e)); |
|
1743 int right_blossom = _blossom_set->find(_graph.v(e)); |
|
1744 |
|
1745 if (left_blossom == right_blossom) { |
|
1746 _delta3->pop(); |
|
1747 } else { |
|
1748 int left_tree; |
|
1749 if ((*_blossom_data)[left_blossom].status == EVEN) { |
|
1750 left_tree = _tree_set->find(left_blossom); |
|
1751 } else { |
|
1752 left_tree = -1; |
|
1753 ++unmatched; |
|
1754 } |
|
1755 int right_tree; |
|
1756 if ((*_blossom_data)[right_blossom].status == EVEN) { |
|
1757 right_tree = _tree_set->find(right_blossom); |
|
1758 } else { |
|
1759 right_tree = -1; |
|
1760 ++unmatched; |
|
1761 } |
|
1762 |
|
1763 if (left_tree == right_tree) { |
|
1764 shrinkOnEdge(e, left_tree); |
|
1765 } else { |
|
1766 augmentOnEdge(e); |
|
1767 unmatched -= 2; |
|
1768 } |
|
1769 } |
|
1770 } break; |
|
1771 case D4: |
|
1772 splitBlossom(_delta4->top()); |
|
1773 break; |
|
1774 } |
|
1775 } |
|
1776 extractMatching(); |
|
1777 } |
|
1778 |
|
1779 /// \brief Runs %MaxWeightedMatching algorithm. |
|
1780 /// |
|
1781 /// This method runs the %MaxWeightedMatching algorithm. |
|
1782 /// |
|
1783 /// \note mwm.run() is just a shortcut of the following code. |
|
1784 /// \code |
|
1785 /// mwm.init(); |
|
1786 /// mwm.start(); |
|
1787 /// \endcode |
|
1788 void run() { |
|
1789 init(); |
|
1790 start(); |
|
1791 } |
|
1792 |
|
1793 /// @} |
|
1794 |
|
1795 /// \name Primal solution |
|
1796 /// Functions to get the primal solution, ie. the matching. |
|
1797 |
|
1798 /// @{ |
|
1799 |
|
1800 /// \brief Returns the weight of the matching. |
|
1801 /// |
|
1802 /// Returns the weight of the matching. |
|
1803 Value matchingValue() const { |
|
1804 Value sum = 0; |
|
1805 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
1806 if ((*_matching)[n] != INVALID) { |
|
1807 sum += _weight[(*_matching)[n]]; |
|
1808 } |
|
1809 } |
|
1810 return sum /= 2; |
|
1811 } |
|
1812 |
|
1813 /// \brief Returns the cardinality of the matching. |
|
1814 /// |
|
1815 /// Returns the cardinality of the matching. |
|
1816 int matchingSize() const { |
|
1817 int num = 0; |
|
1818 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
1819 if ((*_matching)[n] != INVALID) { |
|
1820 ++num; |
|
1821 } |
|
1822 } |
|
1823 return num /= 2; |
|
1824 } |
|
1825 |
|
1826 /// \brief Returns true when the edge is in the matching. |
|
1827 /// |
|
1828 /// Returns true when the edge is in the matching. |
|
1829 bool matching(const Edge& edge) const { |
|
1830 return edge == (*_matching)[_graph.u(edge)]; |
|
1831 } |
|
1832 |
|
1833 /// \brief Returns the incident matching arc. |
|
1834 /// |
|
1835 /// Returns the incident matching arc from given node. If the |
|
1836 /// node is not matched then it gives back \c INVALID. |
|
1837 Arc matching(const Node& node) const { |
|
1838 return (*_matching)[node]; |
|
1839 } |
|
1840 |
|
1841 /// \brief Returns the mate of the node. |
|
1842 /// |
|
1843 /// Returns the adjancent node in a mathcing arc. If the node is |
|
1844 /// not matched then it gives back \c INVALID. |
|
1845 Node mate(const Node& node) const { |
|
1846 return (*_matching)[node] != INVALID ? |
|
1847 _graph.target((*_matching)[node]) : INVALID; |
|
1848 } |
|
1849 |
|
1850 /// @} |
|
1851 |
|
1852 /// \name Dual solution |
|
1853 /// Functions to get the dual solution. |
|
1854 |
|
1855 /// @{ |
|
1856 |
|
1857 /// \brief Returns the value of the dual solution. |
|
1858 /// |
|
1859 /// Returns the value of the dual solution. It should be equal to |
|
1860 /// the primal value scaled by \ref dualScale "dual scale". |
|
1861 Value dualValue() const { |
|
1862 Value sum = 0; |
|
1863 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
1864 sum += nodeValue(n); |
|
1865 } |
|
1866 for (int i = 0; i < blossomNum(); ++i) { |
|
1867 sum += blossomValue(i) * (blossomSize(i) / 2); |
|
1868 } |
|
1869 return sum; |
|
1870 } |
|
1871 |
|
1872 /// \brief Returns the value of the node. |
|
1873 /// |
|
1874 /// Returns the the value of the node. |
|
1875 Value nodeValue(const Node& n) const { |
|
1876 return (*_node_potential)[n]; |
|
1877 } |
|
1878 |
|
1879 /// \brief Returns the number of the blossoms in the basis. |
|
1880 /// |
|
1881 /// Returns the number of the blossoms in the basis. |
|
1882 /// \see BlossomIt |
|
1883 int blossomNum() const { |
|
1884 return _blossom_potential.size(); |
|
1885 } |
|
1886 |
|
1887 |
|
1888 /// \brief Returns the number of the nodes in the blossom. |
|
1889 /// |
|
1890 /// Returns the number of the nodes in the blossom. |
|
1891 int blossomSize(int k) const { |
|
1892 return _blossom_potential[k].end - _blossom_potential[k].begin; |
|
1893 } |
|
1894 |
|
1895 /// \brief Returns the value of the blossom. |
|
1896 /// |
|
1897 /// Returns the the value of the blossom. |
|
1898 /// \see BlossomIt |
|
1899 Value blossomValue(int k) const { |
|
1900 return _blossom_potential[k].value; |
|
1901 } |
|
1902 |
|
1903 /// \brief Iterator for obtaining the nodes of the blossom. |
|
1904 /// |
|
1905 /// Iterator for obtaining the nodes of the blossom. This class |
|
1906 /// provides a common lemon style iterator for listing a |
|
1907 /// subset of the nodes. |
|
1908 class BlossomIt { |
|
1909 public: |
|
1910 |
|
1911 /// \brief Constructor. |
|
1912 /// |
|
1913 /// Constructor to get the nodes of the variable. |
|
1914 BlossomIt(const MaxWeightedMatching& algorithm, int variable) |
|
1915 : _algorithm(&algorithm) |
|
1916 { |
|
1917 _index = _algorithm->_blossom_potential[variable].begin; |
|
1918 _last = _algorithm->_blossom_potential[variable].end; |
|
1919 } |
|
1920 |
|
1921 /// \brief Conversion to node. |
|
1922 /// |
|
1923 /// Conversion to node. |
|
1924 operator Node() const { |
|
1925 return _algorithm->_blossom_node_list[_index]; |
|
1926 } |
|
1927 |
|
1928 /// \brief Increment operator. |
|
1929 /// |
|
1930 /// Increment operator. |
|
1931 BlossomIt& operator++() { |
|
1932 ++_index; |
|
1933 return *this; |
|
1934 } |
|
1935 |
|
1936 /// \brief Validity checking |
|
1937 /// |
|
1938 /// Checks whether the iterator is invalid. |
|
1939 bool operator==(Invalid) const { return _index == _last; } |
|
1940 |
|
1941 /// \brief Validity checking |
|
1942 /// |
|
1943 /// Checks whether the iterator is valid. |
|
1944 bool operator!=(Invalid) const { return _index != _last; } |
|
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 perfect 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 /// edges in the graph with maximum overall weight and no two of |
|
1967 /// them shares their ends and covers all nodes. The problem can be |
|
1968 /// formulated with the following 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} |
|
1971 \quad \forall B\in\mathcal{O}\f] */ |
|
1972 /// \f[x_e \ge 0\quad \forall e\in E\f] |
|
1973 /// \f[\max \sum_{e\in E}x_ew_e\f] |
|
1974 /// where \f$\delta(X)\f$ is the set of edges incident to a node in |
|
1975 /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in |
|
1976 /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality |
|
1977 /// subsets of the nodes. |
|
1978 /// |
|
1979 /// The algorithm calculates an optimal matching and a proof of the |
|
1980 /// optimality. The solution of the dual problem can be used to check |
|
1981 /// the result of the algorithm. The dual linear problem is the |
|
1982 /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge |
|
1983 w_{uv} \quad \forall uv\in E\f] */ |
|
1984 /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f] |
|
1985 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}} |
|
1986 \frac{\vert B \vert - 1}{2}z_B\f] */ |
|
1987 /// |
|
1988 /// The algorithm can be executed with \c run() or the \c init() and |
|
1989 /// then the \c start() member functions. After it the matching can |
|
1990 /// be asked with \c matching() or mate() functions. The dual |
|
1991 /// solution can be get with \c nodeValue(), \c blossomNum() and \c |
|
1992 /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt |
|
1993 /// "BlossomIt" nested class which is able to iterate on the nodes |
|
1994 /// of a blossom. If the value type is integral then the dual |
|
1995 /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4". |
|
1996 template <typename GR, |
|
1997 typename WM = typename GR::template EdgeMap<int> > |
|
1998 class MaxWeightedPerfectMatching { |
|
1999 public: |
|
2000 |
|
2001 typedef GR Graph; |
|
2002 typedef WM WeightMap; |
|
2003 typedef typename WeightMap::Value Value; |
|
2004 |
|
2005 /// \brief Scaling factor for dual solution |
|
2006 /// |
|
2007 /// Scaling factor for dual solution, it is equal to 4 or 1 |
|
2008 /// according to the value type. |
|
2009 static const int dualScale = |
|
2010 std::numeric_limits<Value>::is_integer ? 4 : 1; |
|
2011 |
|
2012 typedef typename Graph::template NodeMap<typename Graph::Arc> |
|
2013 MatchingMap; |
|
2014 |
|
2015 private: |
|
2016 |
|
2017 TEMPLATE_GRAPH_TYPEDEFS(Graph); |
|
2018 |
|
2019 typedef typename Graph::template NodeMap<Value> NodePotential; |
|
2020 typedef std::vector<Node> BlossomNodeList; |
|
2021 |
|
2022 struct BlossomVariable { |
|
2023 int begin, end; |
|
2024 Value value; |
|
2025 |
|
2026 BlossomVariable(int _begin, int _end, Value _value) |
|
2027 : begin(_begin), end(_end), value(_value) {} |
|
2028 |
|
2029 }; |
|
2030 |
|
2031 typedef std::vector<BlossomVariable> BlossomPotential; |
|
2032 |
|
2033 const Graph& _graph; |
|
2034 const WeightMap& _weight; |
|
2035 |
|
2036 MatchingMap* _matching; |
|
2037 |
|
2038 NodePotential* _node_potential; |
|
2039 |
|
2040 BlossomPotential _blossom_potential; |
|
2041 BlossomNodeList _blossom_node_list; |
|
2042 |
|
2043 int _node_num; |
|
2044 int _blossom_num; |
|
2045 |
|
2046 typedef RangeMap<int> IntIntMap; |
|
2047 |
|
2048 enum Status { |
|
2049 EVEN = -1, MATCHED = 0, ODD = 1 |
|
2050 }; |
|
2051 |
|
2052 typedef HeapUnionFind<Value, IntNodeMap> BlossomSet; |
|
2053 struct BlossomData { |
|
2054 int tree; |
|
2055 Status status; |
|
2056 Arc pred, next; |
|
2057 Value pot, offset; |
|
2058 }; |
|
2059 |
|
2060 IntNodeMap *_blossom_index; |
|
2061 BlossomSet *_blossom_set; |
|
2062 RangeMap<BlossomData>* _blossom_data; |
|
2063 |
|
2064 IntNodeMap *_node_index; |
|
2065 IntArcMap *_node_heap_index; |
|
2066 |
|
2067 struct NodeData { |
|
2068 |
|
2069 NodeData(IntArcMap& node_heap_index) |
|
2070 : heap(node_heap_index) {} |
|
2071 |
|
2072 int blossom; |
|
2073 Value pot; |
|
2074 BinHeap<Value, IntArcMap> 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 IntEdgeMap *_delta3_index; |
|
2091 BinHeap<Value, IntEdgeMap> *_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 IntNodeMap(_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 IntNodeMap(_graph); |
|
2116 _node_heap_index = new IntArcMap(_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 IntEdgeMap(_graph); |
|
2131 _delta3 = new BinHeap<Value, IntEdgeMap>(*_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 augmentOnEdge(const Edge& edge) { |
|
2465 |
|
2466 int left = _blossom_set->find(_graph.u(edge)); |
|
2467 int right = _blossom_set->find(_graph.v(edge)); |
|
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(edge, true); |
|
2478 (*_blossom_data)[right].next = _graph.direct(edge, 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 shrinkOnEdge(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)[base] = matching; |
|
2745 _blossom_node_list.push_back(base); |
|
2746 (*_node_potential)[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 |
|
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)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP; |
|
2835 } |
|
2836 for (EdgeIt e(_graph); e != INVALID; ++e) { |
|
2837 (*_delta3_index)[e] = _delta3->PRE_HEAP; |
|
2838 } |
|
2839 for (int i = 0; i < _blossom_num; ++i) { |
|
2840 (*_delta2_index)[i] = _delta2->PRE_HEAP; |
|
2841 (*_delta4_index)[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)[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 shrinkOnEdge(e, left_tree); |
|
2928 } else { |
|
2929 augmentOnEdge(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 to 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 edge is in the matching. |
|
2978 /// |
|
2979 /// Returns true when the edge is in the matching. |
|
2980 bool matching(const Edge& edge) const { |
|
2981 return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge; |
|
2982 } |
|
2983 |
|
2984 /// \brief Returns the incident matching edge. |
|
2985 /// |
|
2986 /// Returns the incident matching arc from given edge. |
|
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 to 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 Iterator for obtaining the nodes of the blossom. |
|
3052 /// |
|
3053 /// Iterator for obtaining the nodes of the blossom. This class |
|
3054 /// provides a common lemon style iterator for listing a |
|
3055 /// subset of the nodes. |
|
3056 class BlossomIt { |
|
3057 public: |
|
3058 |
|
3059 /// \brief Constructor. |
|
3060 /// |
|
3061 /// Constructor to 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 Conversion to node. |
|
3070 /// |
|
3071 /// Conversion to node. |
|
3072 operator Node() const { |
|
3073 return _algorithm->_blossom_node_list[_index]; |
|
3074 } |
|
3075 |
|
3076 /// \brief Increment operator. |
|
3077 /// |
|
3078 /// Increment operator. |
|
3079 BlossomIt& operator++() { |
|
3080 ++_index; |
|
3081 return *this; |
|
3082 } |
|
3083 |
|
3084 /// \brief Validity checking |
|
3085 /// |
|
3086 /// Checks whether the iterator is invalid. |
|
3087 bool operator==(Invalid) const { return _index == _last; } |
|
3088 |
|
3089 /// \brief Validity checking |
|
3090 /// |
|
3091 /// Checks whether the iterator is valid. |
|
3092 bool operator!=(Invalid) const { return _index != _last; } |
|
3093 |
|
3094 private: |
|
3095 const MaxWeightedPerfectMatching* _algorithm; |
|
3096 int _last; |
|
3097 int _index; |
|
3098 }; |
|
3099 |
|
3100 /// @} |
|
3101 |
|
3102 }; |
|
3103 |
|
3104 |
|
3105 } //END OF NAMESPACE LEMON |
|
3106 |
|
3107 #endif //LEMON_MAX_MATCHING_H |
|