29 #include <lemon/bin_heap.h> |
29 #include <lemon/bin_heap.h> |
30 #include <lemon/maps.h> |
30 #include <lemon/maps.h> |
31 |
31 |
32 ///\ingroup matching |
32 ///\ingroup matching |
33 ///\file |
33 ///\file |
34 ///\brief Maximum matching algorithms in graph. |
34 ///\brief Maximum matching algorithms in general graphs. |
35 |
35 |
36 namespace lemon { |
36 namespace lemon { |
37 |
37 |
38 ///\ingroup matching |
38 /// \ingroup matching |
39 /// |
39 /// |
40 ///\brief Edmonds' alternating forest maximum matching algorithm. |
40 /// \brief Edmonds' alternating forest maximum matching algorithm. |
41 /// |
41 /// |
42 ///This class provides Edmonds' alternating forest matching |
42 /// This class provides Edmonds' alternating forest matching |
43 ///algorithm. The starting matching (if any) can be passed to the |
43 /// algorithm. The starting matching (if any) can be passed to the |
44 ///algorithm using some of init functions. |
44 /// algorithm using some of init functions. |
45 /// |
45 /// |
46 ///The dual side of a matching is a map of the nodes to |
46 /// The dual side of a matching is a map of the nodes to |
47 ///MaxMatching::DecompType, having values \c D, \c A and \c C |
47 /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c |
48 ///showing the Gallai-Edmonds decomposition of the digraph. The nodes |
48 /// MATCHED/C showing the Gallai-Edmonds decomposition of the |
49 ///in \c D induce a digraph with factor-critical components, the nodes |
49 /// graph. The nodes in \c EVEN/D induce a graph with |
50 ///in \c A form the barrier, and the nodes in \c C induce a digraph |
50 /// factor-critical components, the nodes in \c ODD/A form the |
51 ///having a perfect matching. This decomposition can be attained by |
51 /// barrier, and the nodes in \c MATCHED/C induce a graph having a |
52 ///calling \c decomposition() after running the algorithm. |
52 /// perfect matching. The number of the fractor critical components |
|
53 /// minus the number of barrier nodes is a lower bound on the |
|
54 /// unmatched nodes, and if the matching is optimal this bound is |
|
55 /// tight. This decomposition can be attained by calling \c |
|
56 /// decomposition() after running the algorithm. |
53 /// |
57 /// |
54 ///\param Digraph The graph type the algorithm runs on. |
58 /// \param _Graph The graph type the algorithm runs on. |
55 template <typename Graph> |
59 template <typename _Graph> |
56 class MaxMatching { |
60 class MaxMatching { |
57 |
61 public: |
58 protected: |
62 |
|
63 typedef _Graph 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, which |
|
70 ///shows an upper bound on the size of a maximum matching. The |
|
71 ///nodes with Status \c EVEN/D induce a graph with factor-critical |
|
72 ///components, the nodes in \c ODD/A form the canonical barrier, |
|
73 ///and the nodes in \c MATCHED/C induce a graph having a perfect |
|
74 ///matching. |
|
75 enum Status { |
|
76 EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2 |
|
77 }; |
|
78 |
|
79 typedef typename Graph::template NodeMap<Status> StatusMap; |
|
80 |
|
81 private: |
59 |
82 |
60 TEMPLATE_GRAPH_TYPEDEFS(Graph); |
83 TEMPLATE_GRAPH_TYPEDEFS(Graph); |
61 |
84 |
62 typedef typename Graph::template NodeMap<int> UFECrossRef; |
85 typedef UnionFindEnum<IntNodeMap> BlossomSet; |
63 typedef UnionFindEnum<UFECrossRef> UFE; |
86 typedef ExtendFindEnum<IntNodeMap> TreeSet; |
64 typedef std::vector<Node> NV; |
87 typedef RangeMap<Node> NodeIntMap; |
65 |
88 typedef MatchingMap EarMap; |
66 typedef typename Graph::template NodeMap<int> EFECrossRef; |
89 typedef std::vector<Node> NodeQueue; |
67 typedef ExtendFindEnum<EFECrossRef> EFE; |
90 |
|
91 const Graph& _graph; |
|
92 MatchingMap* _matching; |
|
93 StatusMap* _status; |
|
94 |
|
95 EarMap* _ear; |
|
96 |
|
97 IntNodeMap* _blossom_set_index; |
|
98 BlossomSet* _blossom_set; |
|
99 NodeIntMap* _blossom_rep; |
|
100 |
|
101 IntNodeMap* _tree_set_index; |
|
102 TreeSet* _tree_set; |
|
103 |
|
104 NodeQueue _node_queue; |
|
105 int _process, _postpone, _last; |
|
106 |
|
107 int _node_num; |
|
108 |
|
109 private: |
|
110 |
|
111 void createStructures() { |
|
112 _node_num = countNodes(_graph); |
|
113 if (!_matching) { |
|
114 _matching = new MatchingMap(_graph); |
|
115 } |
|
116 if (!_status) { |
|
117 _status = new StatusMap(_graph); |
|
118 } |
|
119 if (!_ear) { |
|
120 _ear = new EarMap(_graph); |
|
121 } |
|
122 if (!_blossom_set) { |
|
123 _blossom_set_index = new IntNodeMap(_graph); |
|
124 _blossom_set = new BlossomSet(*_blossom_set_index); |
|
125 } |
|
126 if (!_blossom_rep) { |
|
127 _blossom_rep = new NodeIntMap(_node_num); |
|
128 } |
|
129 if (!_tree_set) { |
|
130 _tree_set_index = new IntNodeMap(_graph); |
|
131 _tree_set = new TreeSet(*_tree_set_index); |
|
132 } |
|
133 _node_queue.resize(_node_num); |
|
134 } |
|
135 |
|
136 void destroyStructures() { |
|
137 if (_matching) { |
|
138 delete _matching; |
|
139 } |
|
140 if (_status) { |
|
141 delete _status; |
|
142 } |
|
143 if (_ear) { |
|
144 delete _ear; |
|
145 } |
|
146 if (_blossom_set) { |
|
147 delete _blossom_set; |
|
148 delete _blossom_set_index; |
|
149 } |
|
150 if (_blossom_rep) { |
|
151 delete _blossom_rep; |
|
152 } |
|
153 if (_tree_set) { |
|
154 delete _tree_set_index; |
|
155 delete _tree_set; |
|
156 } |
|
157 } |
|
158 |
|
159 void processDense(const Node& n) { |
|
160 _process = _postpone = _last = 0; |
|
161 _node_queue[_last++] = n; |
|
162 |
|
163 while (_process != _last) { |
|
164 Node u = _node_queue[_process++]; |
|
165 for (OutArcIt a(_graph, u); a != INVALID; ++a) { |
|
166 Node v = _graph.target(a); |
|
167 if ((*_status)[v] == MATCHED) { |
|
168 extendOnArc(a); |
|
169 } else if ((*_status)[v] == UNMATCHED) { |
|
170 augmentOnArc(a); |
|
171 return; |
|
172 } |
|
173 } |
|
174 } |
|
175 |
|
176 while (_postpone != _last) { |
|
177 Node u = _node_queue[_postpone++]; |
|
178 |
|
179 for (OutArcIt a(_graph, u); a != INVALID ; ++a) { |
|
180 Node v = _graph.target(a); |
|
181 |
|
182 if ((*_status)[v] == EVEN) { |
|
183 if (_blossom_set->find(u) != _blossom_set->find(v)) { |
|
184 shrinkOnEdge(a); |
|
185 } |
|
186 } |
|
187 |
|
188 while (_process != _last) { |
|
189 Node w = _node_queue[_process++]; |
|
190 for (OutArcIt b(_graph, w); b != INVALID; ++b) { |
|
191 Node x = _graph.target(b); |
|
192 if ((*_status)[x] == MATCHED) { |
|
193 extendOnArc(b); |
|
194 } else if ((*_status)[x] == UNMATCHED) { |
|
195 augmentOnArc(b); |
|
196 return; |
|
197 } |
|
198 } |
|
199 } |
|
200 } |
|
201 } |
|
202 } |
|
203 |
|
204 void processSparse(const Node& n) { |
|
205 _process = _last = 0; |
|
206 _node_queue[_last++] = n; |
|
207 while (_process != _last) { |
|
208 Node u = _node_queue[_process++]; |
|
209 for (OutArcIt a(_graph, u); a != INVALID; ++a) { |
|
210 Node v = _graph.target(a); |
|
211 |
|
212 if ((*_status)[v] == EVEN) { |
|
213 if (_blossom_set->find(u) != _blossom_set->find(v)) { |
|
214 shrinkOnEdge(a); |
|
215 } |
|
216 } else if ((*_status)[v] == MATCHED) { |
|
217 extendOnArc(a); |
|
218 } else if ((*_status)[v] == UNMATCHED) { |
|
219 augmentOnArc(a); |
|
220 return; |
|
221 } |
|
222 } |
|
223 } |
|
224 } |
|
225 |
|
226 void shrinkOnEdge(const Edge& e) { |
|
227 Node nca = INVALID; |
|
228 |
|
229 { |
|
230 std::set<Node> left_set, right_set; |
|
231 |
|
232 Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))]; |
|
233 left_set.insert(left); |
|
234 |
|
235 Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))]; |
|
236 right_set.insert(right); |
|
237 |
|
238 while (true) { |
|
239 if ((*_matching)[left] == INVALID) break; |
|
240 left = _graph.target((*_matching)[left]); |
|
241 left = (*_blossom_rep)[_blossom_set-> |
|
242 find(_graph.target((*_ear)[left]))]; |
|
243 if (right_set.find(left) != right_set.end()) { |
|
244 nca = left; |
|
245 break; |
|
246 } |
|
247 left_set.insert(left); |
|
248 |
|
249 if ((*_matching)[right] == INVALID) break; |
|
250 right = _graph.target((*_matching)[right]); |
|
251 right = (*_blossom_rep)[_blossom_set-> |
|
252 find(_graph.target((*_ear)[right]))]; |
|
253 if (left_set.find(right) != left_set.end()) { |
|
254 nca = right; |
|
255 break; |
|
256 } |
|
257 right_set.insert(right); |
|
258 } |
|
259 |
|
260 if (nca == INVALID) { |
|
261 if ((*_matching)[left] == INVALID) { |
|
262 nca = right; |
|
263 while (left_set.find(nca) == left_set.end()) { |
|
264 nca = _graph.target((*_matching)[nca]); |
|
265 nca =(*_blossom_rep)[_blossom_set-> |
|
266 find(_graph.target((*_ear)[nca]))]; |
|
267 } |
|
268 } else { |
|
269 nca = left; |
|
270 while (right_set.find(nca) == right_set.end()) { |
|
271 nca = _graph.target((*_matching)[nca]); |
|
272 nca = (*_blossom_rep)[_blossom_set-> |
|
273 find(_graph.target((*_ear)[nca]))]; |
|
274 } |
|
275 } |
|
276 } |
|
277 } |
|
278 |
|
279 { |
|
280 |
|
281 Node node = _graph.u(e); |
|
282 Arc arc = _graph.direct(e, true); |
|
283 Node base = (*_blossom_rep)[_blossom_set->find(node)]; |
|
284 |
|
285 while (base != nca) { |
|
286 _ear->set(node, arc); |
|
287 |
|
288 Node n = node; |
|
289 while (n != base) { |
|
290 n = _graph.target((*_matching)[n]); |
|
291 Arc a = (*_ear)[n]; |
|
292 n = _graph.target(a); |
|
293 _ear->set(n, _graph.oppositeArc(a)); |
|
294 } |
|
295 node = _graph.target((*_matching)[base]); |
|
296 _tree_set->erase(base); |
|
297 _tree_set->erase(node); |
|
298 _blossom_set->insert(node, _blossom_set->find(base)); |
|
299 _status->set(node, EVEN); |
|
300 _node_queue[_last++] = node; |
|
301 arc = _graph.oppositeArc((*_ear)[node]); |
|
302 node = _graph.target((*_ear)[node]); |
|
303 base = (*_blossom_rep)[_blossom_set->find(node)]; |
|
304 _blossom_set->join(_graph.target(arc), base); |
|
305 } |
|
306 } |
|
307 |
|
308 _blossom_rep->set(_blossom_set->find(nca), nca); |
|
309 |
|
310 { |
|
311 |
|
312 Node node = _graph.v(e); |
|
313 Arc arc = _graph.direct(e, false); |
|
314 Node base = (*_blossom_rep)[_blossom_set->find(node)]; |
|
315 |
|
316 while (base != nca) { |
|
317 _ear->set(node, arc); |
|
318 |
|
319 Node n = node; |
|
320 while (n != base) { |
|
321 n = _graph.target((*_matching)[n]); |
|
322 Arc a = (*_ear)[n]; |
|
323 n = _graph.target(a); |
|
324 _ear->set(n, _graph.oppositeArc(a)); |
|
325 } |
|
326 node = _graph.target((*_matching)[base]); |
|
327 _tree_set->erase(base); |
|
328 _tree_set->erase(node); |
|
329 _blossom_set->insert(node, _blossom_set->find(base)); |
|
330 _status->set(node, EVEN); |
|
331 _node_queue[_last++] = node; |
|
332 arc = _graph.oppositeArc((*_ear)[node]); |
|
333 node = _graph.target((*_ear)[node]); |
|
334 base = (*_blossom_rep)[_blossom_set->find(node)]; |
|
335 _blossom_set->join(_graph.target(arc), base); |
|
336 } |
|
337 } |
|
338 |
|
339 _blossom_rep->set(_blossom_set->find(nca), nca); |
|
340 } |
|
341 |
|
342 |
|
343 |
|
344 void extendOnArc(const Arc& a) { |
|
345 Node base = _graph.source(a); |
|
346 Node odd = _graph.target(a); |
|
347 |
|
348 _ear->set(odd, _graph.oppositeArc(a)); |
|
349 Node even = _graph.target((*_matching)[odd]); |
|
350 _blossom_rep->set(_blossom_set->insert(even), even); |
|
351 _status->set(odd, ODD); |
|
352 _status->set(even, EVEN); |
|
353 int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]); |
|
354 _tree_set->insert(odd, tree); |
|
355 _tree_set->insert(even, tree); |
|
356 _node_queue[_last++] = even; |
|
357 |
|
358 } |
|
359 |
|
360 void augmentOnArc(const Arc& a) { |
|
361 Node even = _graph.source(a); |
|
362 Node odd = _graph.target(a); |
|
363 |
|
364 int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]); |
|
365 |
|
366 _matching->set(odd, _graph.oppositeArc(a)); |
|
367 _status->set(odd, MATCHED); |
|
368 |
|
369 Arc arc = (*_matching)[even]; |
|
370 _matching->set(even, a); |
|
371 |
|
372 while (arc != INVALID) { |
|
373 odd = _graph.target(arc); |
|
374 arc = (*_ear)[odd]; |
|
375 even = _graph.target(arc); |
|
376 _matching->set(odd, arc); |
|
377 arc = (*_matching)[even]; |
|
378 _matching->set(even, _graph.oppositeArc((*_matching)[odd])); |
|
379 } |
|
380 |
|
381 for (typename TreeSet::ItemIt it(*_tree_set, tree); |
|
382 it != INVALID; ++it) { |
|
383 if ((*_status)[it] == ODD) { |
|
384 _status->set(it, MATCHED); |
|
385 } else { |
|
386 int blossom = _blossom_set->find(it); |
|
387 for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom); |
|
388 jt != INVALID; ++jt) { |
|
389 _status->set(jt, MATCHED); |
|
390 } |
|
391 _blossom_set->eraseClass(blossom); |
|
392 } |
|
393 } |
|
394 _tree_set->eraseClass(tree); |
|
395 |
|
396 } |
68 |
397 |
69 public: |
398 public: |
70 |
399 |
71 ///\brief Indicates the Gallai-Edmonds decomposition of the digraph. |
400 /// \brief Constructor |
72 /// |
401 /// |
73 ///Indicates the Gallai-Edmonds decomposition of the digraph, which |
402 /// Constructor. |
74 ///shows an upper bound on the size of a maximum matching. The |
403 MaxMatching(const Graph& graph) |
75 ///nodes with DecompType \c D induce a digraph with factor-critical |
404 : _graph(graph), _matching(0), _status(0), _ear(0), |
76 ///components, the nodes in \c A form the canonical barrier, and the |
405 _blossom_set_index(0), _blossom_set(0), _blossom_rep(0), |
77 ///nodes in \c C induce a digraph having a perfect matching. |
406 _tree_set_index(0), _tree_set(0) {} |
78 enum DecompType { |
407 |
79 D=0, |
408 ~MaxMatching() { |
80 A=1, |
409 destroyStructures(); |
81 C=2 |
410 } |
82 }; |
411 |
83 |
412 /// \name Execution control |
84 protected: |
413 /// The simplest way to execute the algorithm is to use the member |
85 |
414 /// \c run() member function. |
86 static const int HEUR_density=2; |
415 /// \n |
87 const Graph& g; |
416 |
88 typename Graph::template NodeMap<Node> _mate; |
417 /// If you need more control on the execution, first you must call |
89 typename Graph::template NodeMap<DecompType> position; |
418 /// \ref init(), \ref greedyInit() or \ref matchingInit() |
90 |
419 /// functions, then you can start the algorithm with the \ref |
91 public: |
420 /// startParse() or startDense() functions. |
92 |
421 |
93 MaxMatching(const Graph& _g) |
422 ///@{ |
94 : g(_g), _mate(_g), position(_g) {} |
423 |
95 |
424 /// \brief Sets the actual matching to the empty matching. |
96 ///\brief Sets the actual matching to the empty matching. |
425 /// |
97 /// |
426 /// Sets the actual matching to the empty matching. |
98 ///Sets the actual matching to the empty matching. |
|
99 /// |
427 /// |
100 void init() { |
428 void init() { |
101 for(NodeIt v(g); v!=INVALID; ++v) { |
429 createStructures(); |
102 _mate.set(v,INVALID); |
430 for(NodeIt n(_graph); n != INVALID; ++n) { |
103 position.set(v,C); |
431 _matching->set(n, INVALID); |
|
432 _status->set(n, UNMATCHED); |
104 } |
433 } |
105 } |
434 } |
106 |
435 |
107 ///\brief Finds a greedy matching for initial matching. |
436 ///\brief Finds a greedy matching for initial matching. |
108 /// |
437 /// |
109 ///For initial matchig it finds a maximal greedy matching. |
438 ///For initial matchig it finds a maximal greedy matching. |
110 void greedyInit() { |
439 void greedyInit() { |
111 for(NodeIt v(g); v!=INVALID; ++v) { |
440 createStructures(); |
112 _mate.set(v,INVALID); |
441 for (NodeIt n(_graph); n != INVALID; ++n) { |
113 position.set(v,C); |
442 _matching->set(n, INVALID); |
114 } |
443 _status->set(n, UNMATCHED); |
115 for(NodeIt v(g); v!=INVALID; ++v) |
444 } |
116 if ( _mate[v]==INVALID ) { |
445 for (NodeIt n(_graph); n != INVALID; ++n) { |
117 for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) { |
446 if ((*_matching)[n] == INVALID) { |
118 Node y=g.runningNode(e); |
447 for (OutArcIt a(_graph, n); a != INVALID ; ++a) { |
119 if ( _mate[y]==INVALID && y!=v ) { |
448 Node v = _graph.target(a); |
120 _mate.set(v,y); |
449 if ((*_matching)[v] == INVALID && v != n) { |
121 _mate.set(y,v); |
450 _matching->set(n, a); |
|
451 _status->set(n, MATCHED); |
|
452 _matching->set(v, _graph.oppositeArc(a)); |
|
453 _status->set(v, MATCHED); |
122 break; |
454 break; |
123 } |
455 } |
124 } |
456 } |
125 } |
457 } |
126 } |
458 } |
127 |
459 } |
128 ///\brief Initialize the matching from each nodes' mate. |
460 |
129 /// |
461 |
130 ///Initialize the matching from a \c Node valued \c Node map. This |
462 /// \brief Initialize the matching from the map containing a matching. |
131 ///map must be \e symmetric, i.e. if \c map[u]==v then \c |
463 /// |
132 ///map[v]==u must hold, and \c uv will be an arc of the initial |
464 /// Initialize the matching from a \c bool valued \c Edge map. This |
133 ///matching. |
465 /// map must have the property that there are no two incident edges |
134 template <typename MateMap> |
466 /// with true value, ie. it contains a matching. |
135 void mateMapInit(MateMap& map) { |
467 /// \return %True if the map contains a matching. |
136 for(NodeIt v(g); v!=INVALID; ++v) { |
|
137 _mate.set(v,map[v]); |
|
138 position.set(v,C); |
|
139 } |
|
140 } |
|
141 |
|
142 ///\brief Initialize the matching from a node map with the |
|
143 ///incident matching arcs. |
|
144 /// |
|
145 ///Initialize the matching from an \c Edge valued \c Node map. \c |
|
146 ///map[v] must be an \c Edge incident to \c v. This map must have |
|
147 ///the property that if \c g.oppositeNode(u,map[u])==v then \c \c |
|
148 ///g.oppositeNode(v,map[v])==u holds, and now some arc joining \c |
|
149 ///u to \c v will be an arc of the matching. |
|
150 template<typename MatchingMap> |
|
151 void matchingMapInit(MatchingMap& map) { |
|
152 for(NodeIt v(g); v!=INVALID; ++v) { |
|
153 position.set(v,C); |
|
154 Edge e=map[v]; |
|
155 if ( e!=INVALID ) |
|
156 _mate.set(v,g.oppositeNode(v,e)); |
|
157 else |
|
158 _mate.set(v,INVALID); |
|
159 } |
|
160 } |
|
161 |
|
162 ///\brief Initialize the matching from the map containing the |
|
163 ///undirected matching arcs. |
|
164 /// |
|
165 ///Initialize the matching from a \c bool valued \c Edge map. This |
|
166 ///map must have the property that there are no two incident arcs |
|
167 ///\c e, \c f with \c map[e]==map[f]==true. The arcs \c e with \c |
|
168 ///map[e]==true form the matching. |
|
169 template <typename MatchingMap> |
468 template <typename MatchingMap> |
170 void matchingInit(MatchingMap& map) { |
469 bool matchingInit(const MatchingMap& matching) { |
171 for(NodeIt v(g); v!=INVALID; ++v) { |
470 createStructures(); |
172 _mate.set(v,INVALID); |
471 |
173 position.set(v,C); |
472 for (NodeIt n(_graph); n != INVALID; ++n) { |
174 } |
473 _matching->set(n, INVALID); |
175 for(EdgeIt e(g); e!=INVALID; ++e) { |
474 _status->set(n, UNMATCHED); |
176 if ( map[e] ) { |
475 } |
177 Node u=g.u(e); |
476 for(EdgeIt e(_graph); e!=INVALID; ++e) { |
178 Node v=g.v(e); |
477 if (matching[e]) { |
179 _mate.set(u,v); |
478 |
180 _mate.set(v,u); |
479 Node u = _graph.u(e); |
181 } |
480 if ((*_matching)[u] != INVALID) return false; |
182 } |
481 _matching->set(u, _graph.direct(e, true)); |
183 } |
482 _status->set(u, MATCHED); |
184 |
483 |
185 |
484 Node v = _graph.v(e); |
186 ///\brief Runs Edmonds' algorithm. |
485 if ((*_matching)[v] != INVALID) return false; |
187 /// |
486 _matching->set(v, _graph.direct(e, false)); |
188 ///Runs Edmonds' algorithm for sparse digraphs (number of arcs < |
487 _status->set(v, MATCHED); |
189 ///2*number of nodes), and a heuristical Edmonds' algorithm with a |
488 } |
190 ///heuristic of postponing shrinks for dense digraphs. |
489 } |
|
490 return true; |
|
491 } |
|
492 |
|
493 /// \brief Starts Edmonds' algorithm |
|
494 /// |
|
495 /// If runs the original Edmonds' algorithm. |
|
496 void startSparse() { |
|
497 for(NodeIt n(_graph); n != INVALID; ++n) { |
|
498 if ((*_status)[n] == UNMATCHED) { |
|
499 (*_blossom_rep)[_blossom_set->insert(n)] = n; |
|
500 _tree_set->insert(n); |
|
501 _status->set(n, EVEN); |
|
502 processSparse(n); |
|
503 } |
|
504 } |
|
505 } |
|
506 |
|
507 /// \brief Starts Edmonds' algorithm. |
|
508 /// |
|
509 /// It runs Edmonds' algorithm with a heuristic of postponing |
|
510 /// shrinks, giving a faster algorithm for dense graphs. |
|
511 void startDense() { |
|
512 for(NodeIt n(_graph); n != INVALID; ++n) { |
|
513 if ((*_status)[n] == UNMATCHED) { |
|
514 (*_blossom_rep)[_blossom_set->insert(n)] = n; |
|
515 _tree_set->insert(n); |
|
516 _status->set(n, EVEN); |
|
517 processDense(n); |
|
518 } |
|
519 } |
|
520 } |
|
521 |
|
522 |
|
523 /// \brief Runs Edmonds' algorithm |
|
524 /// |
|
525 /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>) |
|
526 /// or Edmonds' algorithm with a heuristic of |
|
527 /// postponing shrinks for dense graphs. |
191 void run() { |
528 void run() { |
192 if (countEdges(g) < HEUR_density * countNodes(g)) { |
529 if (countEdges(_graph) < 2 * countNodes(_graph)) { |
193 greedyInit(); |
530 greedyInit(); |
194 startSparse(); |
531 startSparse(); |
195 } else { |
532 } else { |
196 init(); |
533 init(); |
197 startDense(); |
534 startDense(); |
198 } |
535 } |
199 } |
536 } |
200 |
537 |
201 |
538 /// @} |
202 ///\brief Starts Edmonds' algorithm. |
539 |
203 /// |
540 /// \name Primal solution |
204 ///If runs the original Edmonds' algorithm. |
541 /// Functions for get the primal solution, ie. the matching. |
205 void startSparse() { |
542 |
206 |
543 /// @{ |
207 typename Graph::template NodeMap<Node> ear(g,INVALID); |
|
208 //undefined for the base nodes of the blossoms (i.e. for the |
|
209 //representative elements of UFE blossom) and for the nodes in C |
|
210 |
|
211 UFECrossRef blossom_base(g); |
|
212 UFE blossom(blossom_base); |
|
213 NV rep(countNodes(g)); |
|
214 |
|
215 EFECrossRef tree_base(g); |
|
216 EFE tree(tree_base); |
|
217 |
|
218 //If these UFE's would be members of the class then also |
|
219 //blossom_base and tree_base should be a member. |
|
220 |
|
221 //We build only one tree and the other vertices uncovered by the |
|
222 //matching belong to C. (They can be considered as singleton |
|
223 //trees.) If this tree can be augmented or no more |
|
224 //grow/augmentation/shrink is possible then we return to this |
|
225 //"for" cycle. |
|
226 for(NodeIt v(g); v!=INVALID; ++v) { |
|
227 if (position[v]==C && _mate[v]==INVALID) { |
|
228 rep[blossom.insert(v)] = v; |
|
229 tree.insert(v); |
|
230 position.set(v,D); |
|
231 normShrink(v, ear, blossom, rep, tree); |
|
232 } |
|
233 } |
|
234 } |
|
235 |
|
236 ///\brief Starts Edmonds' algorithm. |
|
237 /// |
|
238 ///It runs Edmonds' algorithm with a heuristic of postponing |
|
239 ///shrinks, giving a faster algorithm for dense digraphs. |
|
240 void startDense() { |
|
241 |
|
242 typename Graph::template NodeMap<Node> ear(g,INVALID); |
|
243 //undefined for the base nodes of the blossoms (i.e. for the |
|
244 //representative elements of UFE blossom) and for the nodes in C |
|
245 |
|
246 UFECrossRef blossom_base(g); |
|
247 UFE blossom(blossom_base); |
|
248 NV rep(countNodes(g)); |
|
249 |
|
250 EFECrossRef tree_base(g); |
|
251 EFE tree(tree_base); |
|
252 |
|
253 //If these UFE's would be members of the class then also |
|
254 //blossom_base and tree_base should be a member. |
|
255 |
|
256 //We build only one tree and the other vertices uncovered by the |
|
257 //matching belong to C. (They can be considered as singleton |
|
258 //trees.) If this tree can be augmented or no more |
|
259 //grow/augmentation/shrink is possible then we return to this |
|
260 //"for" cycle. |
|
261 for(NodeIt v(g); v!=INVALID; ++v) { |
|
262 if ( position[v]==C && _mate[v]==INVALID ) { |
|
263 rep[blossom.insert(v)] = v; |
|
264 tree.insert(v); |
|
265 position.set(v,D); |
|
266 lateShrink(v, ear, blossom, rep, tree); |
|
267 } |
|
268 } |
|
269 } |
|
270 |
|
271 |
|
272 |
544 |
273 ///\brief Returns the size of the actual matching stored. |
545 ///\brief Returns the size of the actual matching stored. |
274 /// |
546 /// |
275 ///Returns the size of the actual matching stored. After \ref |
547 ///Returns the size of the actual matching stored. After \ref |
276 ///run() it returns the size of a maximum matching in the digraph. |
548 ///run() it returns the size of the maximum matching in the graph. |
277 int size() const { |
549 int matchingSize() const { |
278 int s=0; |
550 int size = 0; |
279 for(NodeIt v(g); v!=INVALID; ++v) { |
551 for (NodeIt n(_graph); n != INVALID; ++n) { |
280 if ( _mate[v]!=INVALID ) { |
552 if ((*_matching)[n] != INVALID) { |
281 ++s; |
553 ++size; |
282 } |
554 } |
283 } |
555 } |
284 return s/2; |
556 return size / 2; |
285 } |
557 } |
286 |
558 |
|
559 /// \brief Returns true when the edge is in the matching. |
|
560 /// |
|
561 /// Returns true when the edge is in the matching. |
|
562 bool matching(const Edge& edge) const { |
|
563 return edge == (*_matching)[_graph.u(edge)]; |
|
564 } |
|
565 |
|
566 /// \brief Returns the matching edge incident to the given node. |
|
567 /// |
|
568 /// Returns the matching edge of a \c node in the actual matching or |
|
569 /// INVALID if the \c node is not covered by the actual matching. |
|
570 Arc matching(const Node& n) const { |
|
571 return (*_matching)[n]; |
|
572 } |
287 |
573 |
288 ///\brief Returns the mate of a node in the actual matching. |
574 ///\brief Returns the mate of a node in the actual matching. |
289 /// |
575 /// |
290 ///Returns the mate of a \c node in the actual matching. |
576 ///Returns the mate of a \c node in the actual matching or |
291 ///Returns INVALID if the \c node is not covered by the actual matching. |
577 ///INVALID if the \c node is not covered by the actual matching. |
292 Node mate(const Node& node) const { |
578 Node mate(const Node& n) const { |
293 return _mate[node]; |
579 return (*_matching)[n] != INVALID ? |
294 } |
580 _graph.target((*_matching)[n]) : INVALID; |
295 |
581 } |
296 ///\brief Returns the matching arc incident to the given node. |
582 |
297 /// |
583 /// @} |
298 ///Returns the matching arc of a \c node in the actual matching. |
584 |
299 ///Returns INVALID if the \c node is not covered by the actual matching. |
585 /// \name Dual solution |
300 Edge matchingArc(const Node& node) const { |
586 /// Functions for get the dual solution, ie. the decomposition. |
301 if (_mate[node] == INVALID) return INVALID; |
587 |
302 Node n = node < _mate[node] ? node : _mate[node]; |
588 /// @{ |
303 for (IncEdgeIt e(g, n); e != INVALID; ++e) { |
|
304 if (g.oppositeNode(n, e) == _mate[n]) { |
|
305 return e; |
|
306 } |
|
307 } |
|
308 return INVALID; |
|
309 } |
|
310 |
589 |
311 /// \brief Returns the class of the node in the Edmonds-Gallai |
590 /// \brief Returns the class of the node in the Edmonds-Gallai |
312 /// decomposition. |
591 /// decomposition. |
313 /// |
592 /// |
314 /// Returns the class of the node in the Edmonds-Gallai |
593 /// Returns the class of the node in the Edmonds-Gallai |
315 /// decomposition. |
594 /// decomposition. |
316 DecompType decomposition(const Node& n) { |
595 Status decomposition(const Node& n) const { |
317 return position[n] == A; |
596 return (*_status)[n]; |
318 } |
597 } |
319 |
598 |
320 /// \brief Returns true when the node is in the barrier. |
599 /// \brief Returns true when the node is in the barrier. |
321 /// |
600 /// |
322 /// Returns true when the node is in the barrier. |
601 /// Returns true when the node is in the barrier. |
323 bool barrier(const Node& n) { |
602 bool barrier(const Node& n) const { |
324 return position[n] == A; |
603 return (*_status)[n] == ODD; |
325 } |
604 } |
326 |
605 |
327 ///\brief Gives back the matching in a \c Node of mates. |
606 /// @} |
328 /// |
|
329 ///Writes the stored matching to a \c Node valued \c Node map. The |
|
330 ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c |
|
331 ///map[v]==u will hold, and now \c uv is an arc of the matching. |
|
332 template <typename MateMap> |
|
333 void mateMap(MateMap& map) const { |
|
334 for(NodeIt v(g); v!=INVALID; ++v) { |
|
335 map.set(v,_mate[v]); |
|
336 } |
|
337 } |
|
338 |
|
339 ///\brief Gives back the matching in an \c Edge valued \c Node |
|
340 ///map. |
|
341 /// |
|
342 ///Writes the stored matching to an \c Edge valued \c Node |
|
343 ///map. \c map[v] will be an \c Edge incident to \c v. This |
|
344 ///map will have the property that if \c g.oppositeNode(u,map[u]) |
|
345 ///== v then \c map[u]==map[v] holds, and now this arc is an arc |
|
346 ///of the matching. |
|
347 template <typename MatchingMap> |
|
348 void matchingMap(MatchingMap& map) const { |
|
349 typename Graph::template NodeMap<bool> todo(g,true); |
|
350 for(NodeIt v(g); v!=INVALID; ++v) { |
|
351 if (_mate[v]!=INVALID && v < _mate[v]) { |
|
352 Node u=_mate[v]; |
|
353 for(IncEdgeIt e(g,v); e!=INVALID; ++e) { |
|
354 if ( g.runningNode(e) == u ) { |
|
355 map.set(u,e); |
|
356 map.set(v,e); |
|
357 todo.set(u,false); |
|
358 todo.set(v,false); |
|
359 break; |
|
360 } |
|
361 } |
|
362 } |
|
363 } |
|
364 } |
|
365 |
|
366 |
|
367 ///\brief Gives back the matching in a \c bool valued \c Edge |
|
368 ///map. |
|
369 /// |
|
370 ///Writes the matching stored to a \c bool valued \c Arc |
|
371 ///map. This map will have the property that there are no two |
|
372 ///incident arcs \c e, \c f with \c map[e]==map[f]==true. The |
|
373 ///arcs \c e with \c map[e]==true form the matching. |
|
374 template<typename MatchingMap> |
|
375 void matching(MatchingMap& map) const { |
|
376 for(EdgeIt e(g); e!=INVALID; ++e) map.set(e,false); |
|
377 |
|
378 typename Graph::template NodeMap<bool> todo(g,true); |
|
379 for(NodeIt v(g); v!=INVALID; ++v) { |
|
380 if ( todo[v] && _mate[v]!=INVALID ) { |
|
381 Node u=_mate[v]; |
|
382 for(IncEdgeIt e(g,v); e!=INVALID; ++e) { |
|
383 if ( g.runningNode(e) == u ) { |
|
384 map.set(e,true); |
|
385 todo.set(u,false); |
|
386 todo.set(v,false); |
|
387 break; |
|
388 } |
|
389 } |
|
390 } |
|
391 } |
|
392 } |
|
393 |
|
394 |
|
395 ///\brief Returns the canonical decomposition of the digraph after running |
|
396 ///the algorithm. |
|
397 /// |
|
398 ///After calling any run methods of the class, it writes the |
|
399 ///Gallai-Edmonds canonical decomposition of the digraph. \c map |
|
400 ///must be a node map of \ref DecompType 's. |
|
401 template <typename DecompositionMap> |
|
402 void decomposition(DecompositionMap& map) const { |
|
403 for(NodeIt v(g); v!=INVALID; ++v) map.set(v,position[v]); |
|
404 } |
|
405 |
|
406 ///\brief Returns a barrier on the nodes. |
|
407 /// |
|
408 ///After calling any run methods of the class, it writes a |
|
409 ///canonical barrier on the nodes. The odd component number of the |
|
410 ///remaining digraph minus the barrier size is a lower bound for the |
|
411 ///uncovered nodes in the digraph. The \c map must be a node map of |
|
412 ///bools. |
|
413 template <typename BarrierMap> |
|
414 void barrier(BarrierMap& barrier) { |
|
415 for(NodeIt v(g); v!=INVALID; ++v) barrier.set(v,position[v] == A); |
|
416 } |
|
417 |
|
418 private: |
|
419 |
|
420 |
|
421 void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear, |
|
422 UFE& blossom, NV& rep, EFE& tree) { |
|
423 //We have one tree which we grow, and also shrink but only if it |
|
424 //cannot be postponed. If we augment then we return to the "for" |
|
425 //cycle of runEdmonds(). |
|
426 |
|
427 std::queue<Node> Q; //queue of the totally unscanned nodes |
|
428 Q.push(v); |
|
429 std::queue<Node> R; |
|
430 //queue of the nodes which must be scanned for a possible shrink |
|
431 |
|
432 while ( !Q.empty() ) { |
|
433 Node x=Q.front(); |
|
434 Q.pop(); |
|
435 for( IncEdgeIt e(g,x); e!= INVALID; ++e ) { |
|
436 Node y=g.runningNode(e); |
|
437 //growOrAugment grows if y is covered by the matching and |
|
438 //augments if not. In this latter case it returns 1. |
|
439 if (position[y]==C && |
|
440 growOrAugment(y, x, ear, blossom, rep, tree, Q)) return; |
|
441 } |
|
442 R.push(x); |
|
443 } |
|
444 |
|
445 while ( !R.empty() ) { |
|
446 Node x=R.front(); |
|
447 R.pop(); |
|
448 |
|
449 for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) { |
|
450 Node y=g.runningNode(e); |
|
451 |
|
452 if ( position[y] == D && blossom.find(x) != blossom.find(y) ) |
|
453 //Recall that we have only one tree. |
|
454 shrink( x, y, ear, blossom, rep, tree, Q); |
|
455 |
|
456 while ( !Q.empty() ) { |
|
457 Node z=Q.front(); |
|
458 Q.pop(); |
|
459 for( IncEdgeIt f(g,z); f!= INVALID; ++f ) { |
|
460 Node w=g.runningNode(f); |
|
461 //growOrAugment grows if y is covered by the matching and |
|
462 //augments if not. In this latter case it returns 1. |
|
463 if (position[w]==C && |
|
464 growOrAugment(w, z, ear, blossom, rep, tree, Q)) return; |
|
465 } |
|
466 R.push(z); |
|
467 } |
|
468 } //for e |
|
469 } // while ( !R.empty() ) |
|
470 } |
|
471 |
|
472 void normShrink(Node v, typename Graph::template NodeMap<Node>& ear, |
|
473 UFE& blossom, NV& rep, EFE& tree) { |
|
474 //We have one tree, which we grow and shrink. If we augment then we |
|
475 //return to the "for" cycle of runEdmonds(). |
|
476 |
|
477 std::queue<Node> Q; //queue of the unscanned nodes |
|
478 Q.push(v); |
|
479 while ( !Q.empty() ) { |
|
480 |
|
481 Node x=Q.front(); |
|
482 Q.pop(); |
|
483 |
|
484 for( IncEdgeIt e(g,x); e!=INVALID; ++e ) { |
|
485 Node y=g.runningNode(e); |
|
486 |
|
487 switch ( position[y] ) { |
|
488 case D: //x and y must be in the same tree |
|
489 if ( blossom.find(x) != blossom.find(y)) |
|
490 //x and y are in the same tree |
|
491 shrink(x, y, ear, blossom, rep, tree, Q); |
|
492 break; |
|
493 case C: |
|
494 //growOrAugment grows if y is covered by the matching and |
|
495 //augments if not. In this latter case it returns 1. |
|
496 if (growOrAugment(y, x, ear, blossom, rep, tree, Q)) return; |
|
497 break; |
|
498 default: break; |
|
499 } |
|
500 } |
|
501 } |
|
502 } |
|
503 |
|
504 void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear, |
|
505 UFE& blossom, NV& rep, EFE& tree,std::queue<Node>& Q) { |
|
506 //x and y are the two adjacent vertices in two blossoms. |
|
507 |
|
508 typename Graph::template NodeMap<bool> path(g,false); |
|
509 |
|
510 Node b=rep[blossom.find(x)]; |
|
511 path.set(b,true); |
|
512 b=_mate[b]; |
|
513 while ( b!=INVALID ) { |
|
514 b=rep[blossom.find(ear[b])]; |
|
515 path.set(b,true); |
|
516 b=_mate[b]; |
|
517 } //we go until the root through bases of blossoms and odd vertices |
|
518 |
|
519 Node top=y; |
|
520 Node middle=rep[blossom.find(top)]; |
|
521 Node bottom=x; |
|
522 while ( !path[middle] ) |
|
523 shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q); |
|
524 //Until we arrive to a node on the path, we update blossom, tree |
|
525 //and the positions of the odd nodes. |
|
526 |
|
527 Node base=middle; |
|
528 top=x; |
|
529 middle=rep[blossom.find(top)]; |
|
530 bottom=y; |
|
531 Node blossom_base=rep[blossom.find(base)]; |
|
532 while ( middle!=blossom_base ) |
|
533 shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q); |
|
534 //Until we arrive to a node on the path, we update blossom, tree |
|
535 //and the positions of the odd nodes. |
|
536 |
|
537 rep[blossom.find(base)] = base; |
|
538 } |
|
539 |
|
540 void shrinkStep(Node& top, Node& middle, Node& bottom, |
|
541 typename Graph::template NodeMap<Node>& ear, |
|
542 UFE& blossom, NV& rep, EFE& tree, std::queue<Node>& Q) { |
|
543 //We traverse a blossom and update everything. |
|
544 |
|
545 ear.set(top,bottom); |
|
546 Node t=top; |
|
547 while ( t!=middle ) { |
|
548 Node u=_mate[t]; |
|
549 t=ear[u]; |
|
550 ear.set(t,u); |
|
551 } |
|
552 bottom=_mate[middle]; |
|
553 position.set(bottom,D); |
|
554 Q.push(bottom); |
|
555 top=ear[bottom]; |
|
556 Node oldmiddle=middle; |
|
557 middle=rep[blossom.find(top)]; |
|
558 tree.erase(bottom); |
|
559 tree.erase(oldmiddle); |
|
560 blossom.insert(bottom); |
|
561 blossom.join(bottom, oldmiddle); |
|
562 blossom.join(top, oldmiddle); |
|
563 } |
|
564 |
|
565 |
|
566 |
|
567 bool growOrAugment(Node& y, Node& x, typename Graph::template |
|
568 NodeMap<Node>& ear, UFE& blossom, NV& rep, EFE& tree, |
|
569 std::queue<Node>& Q) { |
|
570 //x is in a blossom in the tree, y is outside. If y is covered by |
|
571 //the matching we grow, otherwise we augment. In this case we |
|
572 //return 1. |
|
573 |
|
574 if ( _mate[y]!=INVALID ) { //grow |
|
575 ear.set(y,x); |
|
576 Node w=_mate[y]; |
|
577 rep[blossom.insert(w)] = w; |
|
578 position.set(y,A); |
|
579 position.set(w,D); |
|
580 int t = tree.find(rep[blossom.find(x)]); |
|
581 tree.insert(y,t); |
|
582 tree.insert(w,t); |
|
583 Q.push(w); |
|
584 } else { //augment |
|
585 augment(x, ear, blossom, rep, tree); |
|
586 _mate.set(x,y); |
|
587 _mate.set(y,x); |
|
588 return true; |
|
589 } |
|
590 return false; |
|
591 } |
|
592 |
|
593 void augment(Node x, typename Graph::template NodeMap<Node>& ear, |
|
594 UFE& blossom, NV& rep, EFE& tree) { |
|
595 Node v=_mate[x]; |
|
596 while ( v!=INVALID ) { |
|
597 |
|
598 Node u=ear[v]; |
|
599 _mate.set(v,u); |
|
600 Node tmp=v; |
|
601 v=_mate[u]; |
|
602 _mate.set(u,tmp); |
|
603 } |
|
604 int y = tree.find(rep[blossom.find(x)]); |
|
605 for (typename EFE::ItemIt tit(tree, y); tit != INVALID; ++tit) { |
|
606 if ( position[tit] == D ) { |
|
607 int b = blossom.find(tit); |
|
608 for (typename UFE::ItemIt bit(blossom, b); bit != INVALID; ++bit) { |
|
609 position.set(bit, C); |
|
610 } |
|
611 blossom.eraseClass(b); |
|
612 } else position.set(tit, C); |
|
613 } |
|
614 tree.eraseClass(y); |
|
615 |
|
616 } |
|
617 |
607 |
618 }; |
608 }; |
619 |
609 |
620 /// \ingroup matching |
610 /// \ingroup matching |
621 /// |
611 /// |