1 //=======================================================================
2 // Copyright (c) 2018 Yi Ji
3 //
4 // Distributed under the Boost Software License, Version 1.0.
5 // (See accompanying file LICENSE_1_0.txt or copy at
6 // http://www.boost.org/LICENSE_1_0.txt)
7 //
8 //=======================================================================
9
10 #ifndef BOOST_GRAPH_MAXIMUM_WEIGHTED_MATCHING_HPP
11 #define BOOST_GRAPH_MAXIMUM_WEIGHTED_MATCHING_HPP
12
13 #include <algorithm> // for std::iter_swap
14 #include <boost/shared_ptr.hpp>
15 #include <boost/make_shared.hpp>
16 #include <boost/graph/max_cardinality_matching.hpp>
17
18 namespace boost
19 {
20 template < typename Graph, typename MateMap, typename VertexIndexMap >
21 typename property_traits<
22 typename property_map< Graph, edge_weight_t >::type >::value_type
matching_weight_sum(const Graph & g,MateMap mate,VertexIndexMap vm)23 matching_weight_sum(const Graph& g, MateMap mate, VertexIndexMap vm)
24 {
25 typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator_t;
26 typedef
27 typename graph_traits< Graph >::vertex_descriptor vertex_descriptor_t;
28 typedef typename property_traits< typename property_map< Graph,
29 edge_weight_t >::type >::value_type edge_property_t;
30
31 edge_property_t weight_sum = 0;
32 vertex_iterator_t vi, vi_end;
33
34 for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
35 {
36 vertex_descriptor_t v = *vi;
37 if (get(mate, v) != graph_traits< Graph >::null_vertex()
38 && get(vm, v) < get(vm, get(mate, v)))
39 weight_sum += get(edge_weight, g, edge(v, mate[v], g).first);
40 }
41 return weight_sum;
42 }
43
44 template < typename Graph, typename MateMap >
45 inline typename property_traits<
46 typename property_map< Graph, edge_weight_t >::type >::value_type
matching_weight_sum(const Graph & g,MateMap mate)47 matching_weight_sum(const Graph& g, MateMap mate)
48 {
49 return matching_weight_sum(g, mate, get(vertex_index, g));
50 }
51
52 template < typename Graph, typename MateMap, typename VertexIndexMap >
53 class weighted_augmenting_path_finder
54 {
55 public:
56 template < typename T > struct map_vertex_to_
57 {
58 typedef boost::iterator_property_map<
59 typename std::vector< T >::iterator, VertexIndexMap >
60 type;
61 };
62 typedef typename graph::detail::VERTEX_STATE vertex_state_t;
63 typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator_t;
64 typedef
65 typename graph_traits< Graph >::vertex_descriptor vertex_descriptor_t;
66 typedef typename std::vector< vertex_descriptor_t >::const_iterator
67 vertex_vec_iter_t;
68 typedef
69 typename graph_traits< Graph >::out_edge_iterator out_edge_iterator_t;
70 typedef typename graph_traits< Graph >::edge_descriptor edge_descriptor_t;
71 typedef typename graph_traits< Graph >::edge_iterator edge_iterator_t;
72 typedef typename property_traits< typename property_map< Graph,
73 edge_weight_t >::type >::value_type edge_property_t;
74 typedef std::deque< vertex_descriptor_t > vertex_list_t;
75 typedef std::vector< edge_descriptor_t > edge_list_t;
76 typedef typename map_vertex_to_< vertex_descriptor_t >::type
77 vertex_to_vertex_map_t;
78 typedef
79 typename map_vertex_to_< edge_property_t >::type vertex_to_weight_map_t;
80 typedef typename map_vertex_to_< bool >::type vertex_to_bool_map_t;
81 typedef typename map_vertex_to_< std::pair< vertex_descriptor_t,
82 vertex_descriptor_t > >::type vertex_to_pair_map_t;
83 typedef
84 typename map_vertex_to_< std::pair< edge_descriptor_t, bool > >::type
85 vertex_to_edge_map_t;
86 typedef typename map_vertex_to_< vertex_to_edge_map_t >::type
87 vertex_pair_to_edge_map_t;
88
89 class blossom
90 {
91 public:
92 typedef boost::shared_ptr< blossom > blossom_ptr_t;
93 std::vector< blossom_ptr_t > sub_blossoms;
94 edge_property_t dual_var;
95 blossom_ptr_t father;
96
blossom()97 blossom() : dual_var(0), father(blossom_ptr_t()) {}
98
99 // get the base vertex of a blossom by recursively getting
100 // its base sub-blossom, which is always the first one in
101 // sub_blossoms because of how we create and maintain blossoms
get_base() const102 virtual vertex_descriptor_t get_base() const
103 {
104 const blossom* b = this;
105 while (!b->sub_blossoms.empty())
106 b = b->sub_blossoms[0].get();
107 return b->get_base();
108 }
109
110 // set a sub-blossom as a blossom's base by exchanging it
111 // with its first sub-blossom
set_base(const blossom_ptr_t & sub)112 void set_base(const blossom_ptr_t& sub)
113 {
114 for (blossom_iterator_t bi = sub_blossoms.begin();
115 bi != sub_blossoms.end(); ++bi)
116 {
117 if (sub.get() == bi->get())
118 {
119 std::iter_swap(sub_blossoms.begin(), bi);
120 break;
121 }
122 }
123 }
124
125 // get all vertices inside recursively
vertices() const126 virtual std::vector< vertex_descriptor_t > vertices() const
127 {
128 std::vector< vertex_descriptor_t > all_vertices;
129 for (typename std::vector< blossom_ptr_t >::const_iterator bi
130 = sub_blossoms.begin();
131 bi != sub_blossoms.end(); ++bi)
132 {
133 std::vector< vertex_descriptor_t > some_vertices
134 = (*bi)->vertices();
135 all_vertices.insert(all_vertices.end(), some_vertices.begin(),
136 some_vertices.end());
137 }
138 return all_vertices;
139 }
140 };
141
142 // a trivial_blossom only has one vertex and no sub-blossom;
143 // for each vertex v, in_blossom[v] is the trivial_blossom that contains it
144 // directly
145 class trivial_blossom : public blossom
146 {
147 public:
trivial_blossom(vertex_descriptor_t v)148 trivial_blossom(vertex_descriptor_t v) : trivial_vertex(v) {}
get_base() const149 virtual vertex_descriptor_t get_base() const { return trivial_vertex; }
150
vertices() const151 virtual std::vector< vertex_descriptor_t > vertices() const
152 {
153 std::vector< vertex_descriptor_t > all_vertices;
154 all_vertices.push_back(trivial_vertex);
155 return all_vertices;
156 }
157
158 private:
159 vertex_descriptor_t trivial_vertex;
160 };
161
162 typedef boost::shared_ptr< blossom > blossom_ptr_t;
163 typedef typename std::vector< blossom_ptr_t >::iterator blossom_iterator_t;
164 typedef
165 typename map_vertex_to_< blossom_ptr_t >::type vertex_to_blossom_map_t;
166
weighted_augmenting_path_finder(const Graph & arg_g,MateMap arg_mate,VertexIndexMap arg_vm)167 weighted_augmenting_path_finder(
168 const Graph& arg_g, MateMap arg_mate, VertexIndexMap arg_vm)
169 : g(arg_g)
170 , vm(arg_vm)
171 , null_edge(std::pair< edge_descriptor_t, bool >(
172 num_edges(g) == 0 ? edge_descriptor_t() : *edges(g).first, false))
173 , mate_vector(num_vertices(g))
174 , label_S_vector(num_vertices(g), graph_traits< Graph >::null_vertex())
175 , label_T_vector(num_vertices(g), graph_traits< Graph >::null_vertex())
176 , outlet_vector(num_vertices(g), graph_traits< Graph >::null_vertex())
177 , tau_idx_vector(num_vertices(g), graph_traits< Graph >::null_vertex())
178 , dual_var_vector(std::vector< edge_property_t >(
179 num_vertices(g), std::numeric_limits< edge_property_t >::min()))
180 , pi_vector(std::vector< edge_property_t >(
181 num_vertices(g), std::numeric_limits< edge_property_t >::max()))
182 , gamma_vector(std::vector< edge_property_t >(
183 num_vertices(g), std::numeric_limits< edge_property_t >::max()))
184 , tau_vector(std::vector< edge_property_t >(
185 num_vertices(g), std::numeric_limits< edge_property_t >::max()))
186 , in_blossom_vector(num_vertices(g))
187 , old_label_vector(num_vertices(g))
188 , critical_edge_vectors(num_vertices(g),
189 std::vector< std::pair< edge_descriptor_t, bool > >(
190 num_vertices(g), null_edge))
191 ,
192
193 mate(mate_vector.begin(), vm)
194 , label_S(label_S_vector.begin(), vm)
195 , label_T(label_T_vector.begin(), vm)
196 , outlet(outlet_vector.begin(), vm)
197 , tau_idx(tau_idx_vector.begin(), vm)
198 , dual_var(dual_var_vector.begin(), vm)
199 , pi(pi_vector.begin(), vm)
200 , gamma(gamma_vector.begin(), vm)
201 , tau(tau_vector.begin(), vm)
202 , in_blossom(in_blossom_vector.begin(), vm)
203 , old_label(old_label_vector.begin(), vm)
204 {
205 vertex_iterator_t vi, vi_end;
206 edge_iterator_t ei, ei_end;
207
208 edge_property_t max_weight
209 = std::numeric_limits< edge_property_t >::min();
210 for (boost::tie(ei, ei_end) = edges(g); ei != ei_end; ++ei)
211 max_weight = std::max(max_weight, get(edge_weight, g, *ei));
212
213 typename std::vector<
214 std::vector< std::pair< edge_descriptor_t, bool > > >::iterator vei;
215
216 for (boost::tie(vi, vi_end) = vertices(g),
217 vei = critical_edge_vectors.begin();
218 vi != vi_end; ++vi, ++vei)
219 {
220 vertex_descriptor_t u = *vi;
221 mate[u] = get(arg_mate, u);
222 dual_var[u] = 2 * max_weight;
223 in_blossom[u] = boost::make_shared< trivial_blossom >(u);
224 outlet[u] = u;
225 critical_edge_vector.push_back(
226 vertex_to_edge_map_t(vei->begin(), vm));
227 }
228
229 critical_edge
230 = vertex_pair_to_edge_map_t(critical_edge_vector.begin(), vm);
231
232 init();
233 }
234
235 // return the top blossom where v is contained inside
in_top_blossom(vertex_descriptor_t v) const236 blossom_ptr_t in_top_blossom(vertex_descriptor_t v) const
237 {
238 blossom_ptr_t b = in_blossom[v];
239 while (b->father != blossom_ptr_t())
240 b = b->father;
241 return b;
242 }
243
244 // check if vertex v is in blossom b
is_in_blossom(blossom_ptr_t b,vertex_descriptor_t v) const245 bool is_in_blossom(blossom_ptr_t b, vertex_descriptor_t v) const
246 {
247 if (v == graph_traits< Graph >::null_vertex())
248 return false;
249 blossom_ptr_t vb = in_blossom[v]->father;
250 while (vb != blossom_ptr_t())
251 {
252 if (vb.get() == b.get())
253 return true;
254 vb = vb->father;
255 }
256 return false;
257 }
258
259 // return the base vertex of the top blossom that contains v
base_vertex(vertex_descriptor_t v) const260 inline vertex_descriptor_t base_vertex(vertex_descriptor_t v) const
261 {
262 return in_top_blossom(v)->get_base();
263 }
264
265 // add an existed top blossom of base vertex v into new top
266 // blossom b as its sub-blossom
add_sub_blossom(blossom_ptr_t b,vertex_descriptor_t v)267 void add_sub_blossom(blossom_ptr_t b, vertex_descriptor_t v)
268 {
269 blossom_ptr_t sub = in_top_blossom(v);
270 sub->father = b;
271 b->sub_blossoms.push_back(sub);
272 if (sub->sub_blossoms.empty())
273 return;
274 for (blossom_iterator_t bi = top_blossoms.begin();
275 bi != top_blossoms.end(); ++bi)
276 {
277 if (bi->get() == sub.get())
278 {
279 top_blossoms.erase(bi);
280 break;
281 }
282 }
283 }
284
285 // when a top blossom is created or its base vertex getting an S-label,
286 // add all edges incident to this blossom into even_edges
bloom(blossom_ptr_t b)287 void bloom(blossom_ptr_t b)
288 {
289 std::vector< vertex_descriptor_t > vertices_of_b = b->vertices();
290 vertex_vec_iter_t vi;
291 for (vi = vertices_of_b.begin(); vi != vertices_of_b.end(); ++vi)
292 {
293 out_edge_iterator_t oei, oei_end;
294 for (boost::tie(oei, oei_end) = out_edges(*vi, g); oei != oei_end;
295 ++oei)
296 {
297 if (target(*oei, g) != *vi && mate[*vi] != target(*oei, g))
298 even_edges.push_back(*oei);
299 }
300 }
301 }
302
303 // assigning a T-label to a non S-vertex, along with outlet and updating pi
304 // value if updated pi[v] equals zero, augment the matching from its mate
305 // vertex
put_T_label(vertex_descriptor_t v,vertex_descriptor_t T_label,vertex_descriptor_t outlet_v,edge_property_t pi_v)306 void put_T_label(vertex_descriptor_t v, vertex_descriptor_t T_label,
307 vertex_descriptor_t outlet_v, edge_property_t pi_v)
308 {
309 if (label_S[v] != graph_traits< Graph >::null_vertex())
310 return;
311
312 label_T[v] = T_label;
313 outlet[v] = outlet_v;
314 pi[v] = pi_v;
315
316 vertex_descriptor_t v_mate = mate[v];
317 if (pi[v] == 0)
318 {
319 label_T[v_mate] = graph_traits< Graph >::null_vertex();
320 label_S[v_mate] = v;
321 bloom(in_top_blossom(v_mate));
322 }
323 }
324
325 // get the missing T-label for a to-be-expanded base vertex
326 // the missing T-label is the last vertex of the path from outlet[v] to v
missing_label(vertex_descriptor_t b_base)327 std::pair< vertex_descriptor_t, vertex_descriptor_t > missing_label(
328 vertex_descriptor_t b_base)
329 {
330 vertex_descriptor_t missing_outlet = outlet[b_base];
331
332 if (outlet[b_base] == b_base)
333 return std::make_pair(
334 graph_traits< Graph >::null_vertex(), missing_outlet);
335
336 vertex_iterator_t vi, vi_end;
337 for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
338 old_label[*vi] = std::make_pair(label_T[*vi], outlet[*vi]);
339
340 std::pair< vertex_descriptor_t, vertex_state_t > child(
341 outlet[b_base], graph::detail::V_EVEN);
342 blossom_ptr_t b = in_blossom[child.first];
343 for (; b->father->father != blossom_ptr_t(); b = b->father)
344 ;
345 child.first = b->get_base();
346
347 if (child.first == b_base)
348 return std::make_pair(
349 graph_traits< Graph >::null_vertex(), missing_outlet);
350
351 while (true)
352 {
353 std::pair< vertex_descriptor_t, vertex_state_t > child_parent
354 = parent(child, true);
355
356 for (b = in_blossom[child_parent.first];
357 b->father->father != blossom_ptr_t(); b = b->father)
358 ;
359 missing_outlet = child_parent.first;
360 child_parent.first = b->get_base();
361
362 if (child_parent.first == b_base)
363 break;
364 else
365 child = child_parent;
366 }
367 return std::make_pair(child.first, missing_outlet);
368 }
369
370 // expand a top blossom, put all its non-trivial sub-blossoms into
371 // top_blossoms
expand_blossom(blossom_iterator_t bi,std::vector<blossom_ptr_t> & new_ones)372 blossom_iterator_t expand_blossom(
373 blossom_iterator_t bi, std::vector< blossom_ptr_t >& new_ones)
374 {
375 blossom_ptr_t b = *bi;
376 for (blossom_iterator_t i = b->sub_blossoms.begin();
377 i != b->sub_blossoms.end(); ++i)
378 {
379 blossom_ptr_t sub_blossom = *i;
380 vertex_descriptor_t sub_base = sub_blossom->get_base();
381 label_S[sub_base] = label_T[sub_base]
382 = graph_traits< Graph >::null_vertex();
383 outlet[sub_base] = sub_base;
384 sub_blossom->father = blossom_ptr_t();
385 // new top blossoms cannot be pushed back into top_blossoms
386 // immediately, because push_back() may cause reallocation and then
387 // invalid iterators
388 if (!sub_blossom->sub_blossoms.empty())
389 new_ones.push_back(sub_blossom);
390 }
391 return top_blossoms.erase(bi);
392 }
393
394 // when expanding a T-blossom with base v, it requires more operations:
395 // supply the missing T-labels for new base vertices by picking the minimum
396 // tau from vertices of each corresponding new top-blossoms; when label_T[v]
397 // is null or we have a smaller tau from missing_label(v), replace T-label
398 // and outlet of v (but don't bloom v)
expand_T_blossom(blossom_iterator_t bi,std::vector<blossom_ptr_t> & new_ones)399 blossom_iterator_t expand_T_blossom(
400 blossom_iterator_t bi, std::vector< blossom_ptr_t >& new_ones)
401 {
402 blossom_ptr_t b = *bi;
403
404 vertex_descriptor_t b_base = b->get_base();
405 std::pair< vertex_descriptor_t, vertex_descriptor_t > T_and_outlet
406 = missing_label(b_base);
407
408 blossom_iterator_t next_bi = expand_blossom(bi, new_ones);
409
410 for (blossom_iterator_t i = b->sub_blossoms.begin();
411 i != b->sub_blossoms.end(); ++i)
412 {
413 blossom_ptr_t sub_blossom = *i;
414 vertex_descriptor_t sub_base = sub_blossom->get_base();
415 vertex_descriptor_t min_tau_v
416 = graph_traits< Graph >::null_vertex();
417 edge_property_t min_tau
418 = std::numeric_limits< edge_property_t >::max();
419
420 std::vector< vertex_descriptor_t > sub_vertices
421 = sub_blossom->vertices();
422 for (vertex_vec_iter_t v = sub_vertices.begin();
423 v != sub_vertices.end(); ++v)
424 {
425 if (tau[*v] < min_tau)
426 {
427 min_tau = tau[*v];
428 min_tau_v = *v;
429 }
430 }
431
432 if (min_tau < std::numeric_limits< edge_property_t >::max())
433 put_T_label(
434 sub_base, tau_idx[min_tau_v], min_tau_v, tau[min_tau_v]);
435 }
436
437 if (label_T[b_base] == graph_traits< Graph >::null_vertex()
438 || tau[old_label[b_base].second] < pi[b_base])
439 boost::tie(label_T[b_base], outlet[b_base]) = T_and_outlet;
440
441 return next_bi;
442 }
443
444 // when vertices v and w are matched to each other by augmenting,
445 // we must set v/w as base vertex of any blossom who contains v/w and
446 // is a sub-blossom of their lowest (smallest) common blossom
adjust_blossom(vertex_descriptor_t v,vertex_descriptor_t w)447 void adjust_blossom(vertex_descriptor_t v, vertex_descriptor_t w)
448 {
449 blossom_ptr_t vb = in_blossom[v], wb = in_blossom[w],
450 lowest_common_blossom;
451 std::vector< blossom_ptr_t > v_ancestors, w_ancestors;
452
453 while (vb->father != blossom_ptr_t())
454 {
455 v_ancestors.push_back(vb->father);
456 vb = vb->father;
457 }
458 while (wb->father != blossom_ptr_t())
459 {
460 w_ancestors.push_back(wb->father);
461 wb = wb->father;
462 }
463
464 typename std::vector< blossom_ptr_t >::reverse_iterator i, j;
465 i = v_ancestors.rbegin();
466 j = w_ancestors.rbegin();
467 while (i != v_ancestors.rend() && j != w_ancestors.rend()
468 && i->get() == j->get())
469 {
470 lowest_common_blossom = *i;
471 ++i;
472 ++j;
473 }
474
475 vb = in_blossom[v];
476 wb = in_blossom[w];
477 while (vb->father != lowest_common_blossom)
478 {
479 vb->father->set_base(vb);
480 vb = vb->father;
481 }
482 while (wb->father != lowest_common_blossom)
483 {
484 wb->father->set_base(wb);
485 wb = wb->father;
486 }
487 }
488
489 // every edge weight is multiplied by 4 to ensure integer weights
490 // throughout the algorithm if all input weights are integers
slack(const edge_descriptor_t & e) const491 inline edge_property_t slack(const edge_descriptor_t& e) const
492 {
493 vertex_descriptor_t v, w;
494 v = source(e, g);
495 w = target(e, g);
496 return dual_var[v] + dual_var[w] - 4 * get(edge_weight, g, e);
497 }
498
499 // backtrace one step on vertex v along the augmenting path
500 // by its labels and its vertex state;
501 // boolean parameter "use_old" means whether we are updating labels,
502 // if we are, then we use old labels to backtrace and also we
503 // don't jump to its base vertex when we reach an odd vertex
parent(std::pair<vertex_descriptor_t,vertex_state_t> v,bool use_old=false) const504 std::pair< vertex_descriptor_t, vertex_state_t > parent(
505 std::pair< vertex_descriptor_t, vertex_state_t > v,
506 bool use_old = false) const
507 {
508 if (v.second == graph::detail::V_EVEN)
509 {
510 // a paranoid check: label_S shoule be the same as mate in
511 // backtracing
512 if (label_S[v.first] == graph_traits< Graph >::null_vertex())
513 label_S[v.first] = mate[v.first];
514 return std::make_pair(label_S[v.first], graph::detail::V_ODD);
515 }
516 else if (v.second == graph::detail::V_ODD)
517 {
518 vertex_descriptor_t w = use_old ? old_label[v.first].first
519 : base_vertex(label_T[v.first]);
520 return std::make_pair(w, graph::detail::V_EVEN);
521 }
522 return std::make_pair(v.first, graph::detail::V_UNREACHED);
523 }
524
525 // backtrace from vertices v and w to their free (unmatched) ancesters,
526 // return the nearest common ancestor (null_vertex if none) of v and w
nearest_common_ancestor(vertex_descriptor_t v,vertex_descriptor_t w,vertex_descriptor_t & v_free_ancestor,vertex_descriptor_t & w_free_ancestor) const527 vertex_descriptor_t nearest_common_ancestor(vertex_descriptor_t v,
528 vertex_descriptor_t w, vertex_descriptor_t& v_free_ancestor,
529 vertex_descriptor_t& w_free_ancestor) const
530 {
531 std::pair< vertex_descriptor_t, vertex_state_t > v_up(
532 v, graph::detail::V_EVEN);
533 std::pair< vertex_descriptor_t, vertex_state_t > w_up(
534 w, graph::detail::V_EVEN);
535 vertex_descriptor_t nca;
536 nca = w_free_ancestor = v_free_ancestor
537 = graph_traits< Graph >::null_vertex();
538
539 std::vector< bool > ancestor_of_w_vector(num_vertices(g), false);
540 std::vector< bool > ancestor_of_v_vector(num_vertices(g), false);
541 vertex_to_bool_map_t ancestor_of_w(ancestor_of_w_vector.begin(), vm);
542 vertex_to_bool_map_t ancestor_of_v(ancestor_of_v_vector.begin(), vm);
543
544 while (nca == graph_traits< Graph >::null_vertex()
545 && (v_free_ancestor == graph_traits< Graph >::null_vertex()
546 || w_free_ancestor == graph_traits< Graph >::null_vertex()))
547 {
548 ancestor_of_w[w_up.first] = true;
549 ancestor_of_v[v_up.first] = true;
550
551 if (w_free_ancestor == graph_traits< Graph >::null_vertex())
552 w_up = parent(w_up);
553 if (v_free_ancestor == graph_traits< Graph >::null_vertex())
554 v_up = parent(v_up);
555
556 if (mate[v_up.first] == graph_traits< Graph >::null_vertex())
557 v_free_ancestor = v_up.first;
558 if (mate[w_up.first] == graph_traits< Graph >::null_vertex())
559 w_free_ancestor = w_up.first;
560
561 if (ancestor_of_w[v_up.first] == true || v_up.first == w_up.first)
562 nca = v_up.first;
563 else if (ancestor_of_v[w_up.first] == true)
564 nca = w_up.first;
565 else if (v_free_ancestor == w_free_ancestor
566 && v_free_ancestor != graph_traits< Graph >::null_vertex())
567 nca = v_up.first;
568 }
569
570 return nca;
571 }
572
573 // when a new top blossom b is created by connecting (v, w), we add
574 // sub-blossoms into b along backtracing from v_prime and w_prime to
575 // stop_vertex (the base vertex); also, we set labels and outlet for each
576 // base vertex we pass by
make_blossom(blossom_ptr_t b,vertex_descriptor_t w_prime,vertex_descriptor_t v_prime,vertex_descriptor_t stop_vertex)577 void make_blossom(blossom_ptr_t b, vertex_descriptor_t w_prime,
578 vertex_descriptor_t v_prime, vertex_descriptor_t stop_vertex)
579 {
580 std::pair< vertex_descriptor_t, vertex_state_t > u(
581 v_prime, graph::detail::V_ODD);
582 std::pair< vertex_descriptor_t, vertex_state_t > u_up(
583 w_prime, graph::detail::V_EVEN);
584
585 for (; u_up.first != stop_vertex; u = u_up, u_up = parent(u))
586 {
587 if (u_up.second == graph::detail::V_EVEN)
588 {
589 if (!in_top_blossom(u_up.first)->sub_blossoms.empty())
590 outlet[u_up.first] = label_T[u.first];
591 label_T[u_up.first] = outlet[u.first];
592 }
593 else if (u_up.second == graph::detail::V_ODD)
594 label_S[u_up.first] = u.first;
595
596 add_sub_blossom(b, u_up.first);
597 }
598 }
599
600 // the design of recursively expanding augmenting path in
601 // (reversed_)retrieve_augmenting_path functions is inspired by same
602 // functions in max_cardinality_matching.hpp; except that in weighted
603 // matching, we use "outlet" vertices instead of "bridge" vertex pairs: if
604 // blossom b is the smallest non-trivial blossom that contains its base
605 // vertex v, then v and outlet[v] are where augmenting path enters and
606 // leaves b
retrieve_augmenting_path(vertex_descriptor_t v,vertex_descriptor_t w,vertex_state_t v_state)607 void retrieve_augmenting_path(
608 vertex_descriptor_t v, vertex_descriptor_t w, vertex_state_t v_state)
609 {
610 if (v == w)
611 aug_path.push_back(v);
612 else if (v_state == graph::detail::V_EVEN)
613 {
614 aug_path.push_back(v);
615 retrieve_augmenting_path(label_S[v], w, graph::detail::V_ODD);
616 }
617 else if (v_state == graph::detail::V_ODD)
618 {
619 if (outlet[v] == v)
620 aug_path.push_back(v);
621 else
622 reversed_retrieve_augmenting_path(
623 outlet[v], v, graph::detail::V_EVEN);
624 retrieve_augmenting_path(label_T[v], w, graph::detail::V_EVEN);
625 }
626 }
627
reversed_retrieve_augmenting_path(vertex_descriptor_t v,vertex_descriptor_t w,vertex_state_t v_state)628 void reversed_retrieve_augmenting_path(
629 vertex_descriptor_t v, vertex_descriptor_t w, vertex_state_t v_state)
630 {
631 if (v == w)
632 aug_path.push_back(v);
633 else if (v_state == graph::detail::V_EVEN)
634 {
635 reversed_retrieve_augmenting_path(
636 label_S[v], w, graph::detail::V_ODD);
637 aug_path.push_back(v);
638 }
639 else if (v_state == graph::detail::V_ODD)
640 {
641 reversed_retrieve_augmenting_path(
642 label_T[v], w, graph::detail::V_EVEN);
643 if (outlet[v] != v)
644 retrieve_augmenting_path(outlet[v], v, graph::detail::V_EVEN);
645 else
646 aug_path.push_back(v);
647 }
648 }
649
650 // correct labels for vertices in the augmenting path
relabel(vertex_descriptor_t v)651 void relabel(vertex_descriptor_t v)
652 {
653 blossom_ptr_t b = in_blossom[v]->father;
654
655 if (!is_in_blossom(b, mate[v]))
656 { // if v is a new base vertex
657 std::pair< vertex_descriptor_t, vertex_state_t > u(
658 v, graph::detail::V_EVEN);
659 while (label_S[u.first] != u.first
660 && is_in_blossom(b, label_S[u.first]))
661 u = parent(u, true);
662
663 vertex_descriptor_t old_base = u.first;
664 if (label_S[old_base] != old_base)
665 { // if old base is not exposed
666 label_T[v] = label_S[old_base];
667 outlet[v] = old_base;
668 }
669 else
670 { // if old base is exposed then new label_T[v] is not in b,
671 // we must (i) make b2 the smallest blossom containing v but not
672 // as base vertex (ii) backtrace from b2's new base vertex to b
673 label_T[v] = graph_traits< Graph >::null_vertex();
674 for (b = b->father; b != blossom_ptr_t() && b->get_base() == v;
675 b = b->father)
676 ;
677 if (b != blossom_ptr_t())
678 {
679 u = std::make_pair(b->get_base(), graph::detail::V_ODD);
680 while (!is_in_blossom(
681 in_blossom[v]->father, old_label[u.first].first))
682 u = parent(u, true);
683 label_T[v] = u.first;
684 outlet[v] = old_label[u.first].first;
685 }
686 }
687 }
688 else if (label_S[v] == v || !is_in_blossom(b, label_S[v]))
689 { // if v is an old base vertex
690 // let u be the new base vertex; backtrace from u's old T-label
691 std::pair< vertex_descriptor_t, vertex_state_t > u(
692 b->get_base(), graph::detail::V_ODD);
693 while (
694 old_label[u.first].first != graph_traits< Graph >::null_vertex()
695 && old_label[u.first].first != v)
696 u = parent(u, true);
697 label_T[v] = old_label[u.first].second;
698 outlet[v] = v;
699 }
700 else // if v is neither a new nor an old base vertex
701 label_T[v] = label_S[v];
702 }
703
augmenting(vertex_descriptor_t v,vertex_descriptor_t v_free_ancestor,vertex_descriptor_t w,vertex_descriptor_t w_free_ancestor)704 void augmenting(vertex_descriptor_t v, vertex_descriptor_t v_free_ancestor,
705 vertex_descriptor_t w, vertex_descriptor_t w_free_ancestor)
706 {
707 vertex_iterator_t vi, vi_end;
708
709 // retrieve the augmenting path and put it in aug_path
710 reversed_retrieve_augmenting_path(
711 v, v_free_ancestor, graph::detail::V_EVEN);
712 retrieve_augmenting_path(w, w_free_ancestor, graph::detail::V_EVEN);
713
714 // augment the matching along aug_path
715 vertex_descriptor_t a, b;
716 vertex_list_t reversed_aug_path;
717 while (!aug_path.empty())
718 {
719 a = aug_path.front();
720 aug_path.pop_front();
721 reversed_aug_path.push_back(a);
722 b = aug_path.front();
723 aug_path.pop_front();
724 reversed_aug_path.push_back(b);
725
726 mate[a] = b;
727 mate[b] = a;
728
729 // reset base vertex for every blossom in augment path
730 adjust_blossom(a, b);
731 }
732
733 for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
734 old_label[*vi] = std::make_pair(label_T[*vi], outlet[*vi]);
735
736 // correct labels for in-blossom vertices along aug_path
737 while (!reversed_aug_path.empty())
738 {
739 a = reversed_aug_path.front();
740 reversed_aug_path.pop_front();
741
742 if (in_blossom[a]->father != blossom_ptr_t())
743 relabel(a);
744 }
745
746 for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
747 {
748 vertex_descriptor_t u = *vi;
749 if (mate[u] != graph_traits< Graph >::null_vertex())
750 label_S[u] = mate[u];
751 }
752
753 // expand blossoms with zero dual variables
754 std::vector< blossom_ptr_t > new_top_blossoms;
755 for (blossom_iterator_t bi = top_blossoms.begin();
756 bi != top_blossoms.end();)
757 {
758 if ((*bi)->dual_var <= 0)
759 bi = expand_blossom(bi, new_top_blossoms);
760 else
761 ++bi;
762 }
763 top_blossoms.insert(top_blossoms.end(), new_top_blossoms.begin(),
764 new_top_blossoms.end());
765 init();
766 }
767
768 // create a new blossom and set labels for vertices inside
blossoming(vertex_descriptor_t v,vertex_descriptor_t v_prime,vertex_descriptor_t w,vertex_descriptor_t w_prime,vertex_descriptor_t nca)769 void blossoming(vertex_descriptor_t v, vertex_descriptor_t v_prime,
770 vertex_descriptor_t w, vertex_descriptor_t w_prime,
771 vertex_descriptor_t nca)
772 {
773 vertex_iterator_t vi, vi_end;
774
775 std::vector< bool > is_old_base_vector(num_vertices(g));
776 vertex_to_bool_map_t is_old_base(is_old_base_vector.begin(), vm);
777 for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
778 {
779 if (*vi == base_vertex(*vi))
780 is_old_base[*vi] = true;
781 }
782
783 blossom_ptr_t b = boost::make_shared< blossom >();
784 add_sub_blossom(b, nca);
785
786 label_T[w_prime] = v;
787 label_T[v_prime] = w;
788 outlet[w_prime] = w;
789 outlet[v_prime] = v;
790
791 make_blossom(b, w_prime, v_prime, nca);
792 make_blossom(b, v_prime, w_prime, nca);
793
794 label_T[nca] = graph_traits< Graph >::null_vertex();
795 outlet[nca] = nca;
796
797 top_blossoms.push_back(b);
798 bloom(b);
799
800 // set gamma[b_base] = min_slack{critical_edge(b_base, other_base)}
801 // where each critical edge is updated before, by
802 // argmin{slack(old_bases_in_b, other_base)};
803 vertex_vec_iter_t i, j;
804 std::vector< vertex_descriptor_t > b_vertices = b->vertices(),
805 old_base_in_b, other_base;
806 vertex_descriptor_t b_base = b->get_base();
807 for (i = b_vertices.begin(); i != b_vertices.end(); ++i)
808 {
809 if (is_old_base[*i])
810 old_base_in_b.push_back(*i);
811 }
812 for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
813 {
814 if (*vi != b_base && *vi == base_vertex(*vi))
815 other_base.push_back(*vi);
816 }
817 for (i = other_base.begin(); i != other_base.end(); ++i)
818 {
819 edge_property_t min_slack
820 = std::numeric_limits< edge_property_t >::max();
821 std::pair< edge_descriptor_t, bool > b_vi = null_edge;
822 for (j = old_base_in_b.begin(); j != old_base_in_b.end(); ++j)
823 {
824 if (critical_edge[*j][*i] != null_edge
825 && min_slack > slack(critical_edge[*j][*i].first))
826 {
827 min_slack = slack(critical_edge[*j][*i].first);
828 b_vi = critical_edge[*j][*i];
829 }
830 }
831 critical_edge[b_base][*i] = critical_edge[*i][b_base] = b_vi;
832 }
833 gamma[b_base] = std::numeric_limits< edge_property_t >::max();
834 for (i = other_base.begin(); i != other_base.end(); ++i)
835 {
836 if (critical_edge[b_base][*i] != null_edge)
837 gamma[b_base] = std::min(
838 gamma[b_base], slack(critical_edge[b_base][*i].first));
839 }
840 }
841
init()842 void init()
843 {
844 even_edges.clear();
845
846 vertex_iterator_t vi, vi_end;
847 typename std::vector<
848 std::vector< std::pair< edge_descriptor_t, bool > > >::iterator vei;
849
850 for (boost::tie(vi, vi_end) = vertices(g),
851 vei = critical_edge_vectors.begin();
852 vi != vi_end; ++vi, ++vei)
853 {
854 vertex_descriptor_t u = *vi;
855 out_edge_iterator_t ei, ei_end;
856
857 gamma[u] = tau[u] = pi[u]
858 = std::numeric_limits< edge_property_t >::max();
859 std::fill(vei->begin(), vei->end(), null_edge);
860
861 if (base_vertex(u) != u)
862 continue;
863
864 label_S[u] = label_T[u] = graph_traits< Graph >::null_vertex();
865 outlet[u] = u;
866
867 if (mate[u] == graph_traits< Graph >::null_vertex())
868 {
869 label_S[u] = u;
870 bloom(in_top_blossom(u));
871 }
872 }
873 }
874
augment_matching()875 bool augment_matching()
876 {
877 vertex_descriptor_t v, w, w_free_ancestor, v_free_ancestor;
878 v = w = w_free_ancestor = v_free_ancestor
879 = graph_traits< Graph >::null_vertex();
880 bool found_alternating_path = false;
881
882 // note that we only use edges of zero slack value for augmenting
883 while (!even_edges.empty() && !found_alternating_path)
884 {
885 // search for augmenting paths depth-first
886 edge_descriptor_t current_edge = even_edges.back();
887 even_edges.pop_back();
888
889 v = source(current_edge, g);
890 w = target(current_edge, g);
891
892 vertex_descriptor_t v_prime = base_vertex(v);
893 vertex_descriptor_t w_prime = base_vertex(w);
894
895 // w_prime == v_prime implies that we get an edge that has been
896 // shrunk into a blossom
897 if (v_prime == w_prime)
898 continue;
899
900 // a paranoid check
901 if (label_S[v_prime] == graph_traits< Graph >::null_vertex())
902 {
903 std::swap(v_prime, w_prime);
904 std::swap(v, w);
905 }
906
907 // w_prime may be unlabeled or have a T-label; replace the existed
908 // T-label if the edge slack is smaller than current pi[w_prime] and
909 // update it. Note that a T-label is "deserved" only when pi equals
910 // zero. also update tau and tau_idx so that tau_idx becomes T-label
911 // when a T-blossom is expanded
912 if (label_S[w_prime] == graph_traits< Graph >::null_vertex())
913 {
914 if (slack(current_edge) < pi[w_prime])
915 put_T_label(w_prime, v, w, slack(current_edge));
916 if (slack(current_edge) < tau[w])
917 {
918 if (in_blossom[w]->father == blossom_ptr_t()
919 || label_T[w_prime] == v
920 || label_T[w_prime]
921 == graph_traits< Graph >::null_vertex()
922 || nearest_common_ancestor(v_prime, label_T[w_prime],
923 v_free_ancestor, w_free_ancestor)
924 == graph_traits< Graph >::null_vertex())
925 {
926 tau[w] = slack(current_edge);
927 tau_idx[w] = v;
928 }
929 }
930 }
931
932 else
933 {
934 if (slack(current_edge) > 0)
935 {
936 // update gamma and critical_edges when we have a smaller
937 // edge slack
938 gamma[v_prime]
939 = std::min(gamma[v_prime], slack(current_edge));
940 gamma[w_prime]
941 = std::min(gamma[w_prime], slack(current_edge));
942 if (critical_edge[v_prime][w_prime] == null_edge
943 || slack(critical_edge[v_prime][w_prime].first)
944 > slack(current_edge))
945 {
946 critical_edge[v_prime][w_prime]
947 = std::pair< edge_descriptor_t, bool >(
948 current_edge, true);
949 critical_edge[w_prime][v_prime]
950 = std::pair< edge_descriptor_t, bool >(
951 current_edge, true);
952 }
953 continue;
954 }
955 else if (slack(current_edge) == 0)
956 {
957 // if nca is null_vertex then we have an augmenting path;
958 // otherwise we have a new top blossom with nca as its base
959 // vertex
960 vertex_descriptor_t nca = nearest_common_ancestor(
961 v_prime, w_prime, v_free_ancestor, w_free_ancestor);
962
963 if (nca == graph_traits< Graph >::null_vertex())
964 found_alternating_path
965 = true; // to break out of the loop
966 else
967 blossoming(v, v_prime, w, w_prime, nca);
968 }
969 }
970 }
971
972 if (!found_alternating_path)
973 return false;
974
975 augmenting(v, v_free_ancestor, w, w_free_ancestor);
976 return true;
977 }
978
979 // slack the vertex and blossom dual variables when there is no augmenting
980 // path found according to the primal-dual method
adjust_dual()981 bool adjust_dual()
982 {
983 edge_property_t delta1, delta2, delta3, delta4, delta;
984 delta1 = delta2 = delta3 = delta4
985 = std::numeric_limits< edge_property_t >::max();
986
987 vertex_iterator_t vi, vi_end;
988
989 for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
990 {
991 delta1 = std::min(delta1, dual_var[*vi]);
992 delta4 = pi[*vi] > 0 ? std::min(delta4, pi[*vi]) : delta4;
993 if (*vi == base_vertex(*vi))
994 delta3 = std::min(delta3, gamma[*vi] / 2);
995 }
996
997 for (blossom_iterator_t bi = top_blossoms.begin();
998 bi != top_blossoms.end(); ++bi)
999 {
1000 vertex_descriptor_t b_base = (*bi)->get_base();
1001 if (label_T[b_base] != graph_traits< Graph >::null_vertex()
1002 && pi[b_base] == 0)
1003 delta2 = std::min(delta2, (*bi)->dual_var / 2);
1004 }
1005
1006 delta = std::min(std::min(delta1, delta2), std::min(delta3, delta4));
1007
1008 // start updating dual variables, note that the order is important
1009
1010 for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
1011 {
1012 vertex_descriptor_t v = *vi, v_prime = base_vertex(v);
1013
1014 if (label_S[v_prime] != graph_traits< Graph >::null_vertex())
1015 dual_var[v] -= delta;
1016 else if (label_T[v_prime] != graph_traits< Graph >::null_vertex()
1017 && pi[v_prime] == 0)
1018 dual_var[v] += delta;
1019
1020 if (v == v_prime)
1021 gamma[v] -= 2 * delta;
1022 }
1023
1024 for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
1025 {
1026 vertex_descriptor_t v_prime = base_vertex(*vi);
1027 if (pi[v_prime] > 0)
1028 tau[*vi] -= delta;
1029 }
1030
1031 for (blossom_iterator_t bi = top_blossoms.begin();
1032 bi != top_blossoms.end(); ++bi)
1033 {
1034 vertex_descriptor_t b_base = (*bi)->get_base();
1035 if (label_T[b_base] != graph_traits< Graph >::null_vertex()
1036 && pi[b_base] == 0)
1037 (*bi)->dual_var -= 2 * delta;
1038 if (label_S[b_base] != graph_traits< Graph >::null_vertex())
1039 (*bi)->dual_var += 2 * delta;
1040 }
1041
1042 for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
1043 {
1044 vertex_descriptor_t v = *vi;
1045 if (pi[v] > 0)
1046 pi[v] -= delta;
1047
1048 // when some T-vertices have zero pi value, bloom their mates so
1049 // that matching can be further augmented
1050 if (label_T[v] != graph_traits< Graph >::null_vertex()
1051 && pi[v] == 0)
1052 put_T_label(v, label_T[v], outlet[v], pi[v]);
1053 }
1054
1055 // optimal solution reached, halt
1056 if (delta == delta1)
1057 return false;
1058
1059 // expand odd blossoms with zero dual variables and zero pi value of
1060 // their base vertices
1061 if (delta == delta2 && delta != delta3)
1062 {
1063 std::vector< blossom_ptr_t > new_top_blossoms;
1064 for (blossom_iterator_t bi = top_blossoms.begin();
1065 bi != top_blossoms.end();)
1066 {
1067 const blossom_ptr_t b = *bi;
1068 vertex_descriptor_t b_base = b->get_base();
1069 if (b->dual_var == 0
1070 && label_T[b_base] != graph_traits< Graph >::null_vertex()
1071 && pi[b_base] == 0)
1072 bi = expand_T_blossom(bi, new_top_blossoms);
1073 else
1074 ++bi;
1075 }
1076 top_blossoms.insert(top_blossoms.end(), new_top_blossoms.begin(),
1077 new_top_blossoms.end());
1078 }
1079
1080 while (true)
1081 {
1082 // find a zero-slack critical edge (v, w) of zero gamma values
1083 std::pair< edge_descriptor_t, bool > best_edge = null_edge;
1084 std::vector< vertex_descriptor_t > base_nodes;
1085 for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
1086 {
1087 if (*vi == base_vertex(*vi))
1088 base_nodes.push_back(*vi);
1089 }
1090 for (vertex_vec_iter_t i = base_nodes.begin();
1091 i != base_nodes.end(); ++i)
1092 {
1093 if (gamma[*i] == 0)
1094 {
1095 for (vertex_vec_iter_t j = base_nodes.begin();
1096 j != base_nodes.end(); ++j)
1097 {
1098 if (critical_edge[*i][*j] != null_edge
1099 && slack(critical_edge[*i][*j].first) == 0)
1100 best_edge = critical_edge[*i][*j];
1101 }
1102 }
1103 }
1104
1105 // if not found, continue finding other augment matching
1106 if (best_edge == null_edge)
1107 {
1108 bool augmented = augment_matching();
1109 return augmented || delta != delta1;
1110 }
1111 // if found, determine either augmenting or blossoming
1112 vertex_descriptor_t v = source(best_edge.first, g),
1113 w = target(best_edge.first, g);
1114 vertex_descriptor_t v_prime = base_vertex(v),
1115 w_prime = base_vertex(w), v_free_ancestor,
1116 w_free_ancestor;
1117 vertex_descriptor_t nca = nearest_common_ancestor(
1118 v_prime, w_prime, v_free_ancestor, w_free_ancestor);
1119 if (nca == graph_traits< Graph >::null_vertex())
1120 {
1121 augmenting(v, v_free_ancestor, w, w_free_ancestor);
1122 return true;
1123 }
1124 else
1125 blossoming(v, v_prime, w, w_prime, nca);
1126 }
1127
1128 return false;
1129 }
1130
get_current_matching(PropertyMap pm)1131 template < typename PropertyMap > void get_current_matching(PropertyMap pm)
1132 {
1133 vertex_iterator_t vi, vi_end;
1134 for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
1135 put(pm, *vi, mate[*vi]);
1136 }
1137
1138 private:
1139 const Graph& g;
1140 VertexIndexMap vm;
1141 const std::pair< edge_descriptor_t, bool > null_edge;
1142
1143 // storage for the property maps below
1144 std::vector< vertex_descriptor_t > mate_vector;
1145 std::vector< vertex_descriptor_t > label_S_vector, label_T_vector;
1146 std::vector< vertex_descriptor_t > outlet_vector;
1147 std::vector< vertex_descriptor_t > tau_idx_vector;
1148 std::vector< edge_property_t > dual_var_vector;
1149 std::vector< edge_property_t > pi_vector, gamma_vector, tau_vector;
1150 std::vector< blossom_ptr_t > in_blossom_vector;
1151 std::vector< std::pair< vertex_descriptor_t, vertex_descriptor_t > >
1152 old_label_vector;
1153 std::vector< vertex_to_edge_map_t > critical_edge_vector;
1154 std::vector< std::vector< std::pair< edge_descriptor_t, bool > > >
1155 critical_edge_vectors;
1156
1157 // iterator property maps
1158 vertex_to_vertex_map_t mate;
1159 vertex_to_vertex_map_t label_S; // v has an S-label -> v can be an even
1160 // vertex, label_S[v] is its mate
1161 vertex_to_vertex_map_t
1162 label_T; // v has a T-label -> v can be an odd vertex, label_T[v] is its
1163 // predecessor in aug_path
1164 vertex_to_vertex_map_t outlet;
1165 vertex_to_vertex_map_t tau_idx;
1166 vertex_to_weight_map_t dual_var;
1167 vertex_to_weight_map_t pi, gamma, tau;
1168 vertex_to_blossom_map_t
1169 in_blossom; // map any vertex v to the trivial blossom containing v
1170 vertex_to_pair_map_t old_label; // <old T-label, old outlet> before
1171 // relabeling or expanding T-blossoms
1172 vertex_pair_to_edge_map_t
1173 critical_edge; // an not matched edge (v, w) is critical if v and w
1174 // belongs to different S-blossoms
1175
1176 vertex_list_t aug_path;
1177 edge_list_t even_edges;
1178 std::vector< blossom_ptr_t > top_blossoms;
1179 };
1180
1181 template < typename Graph, typename MateMap, typename VertexIndexMap >
maximum_weighted_matching(const Graph & g,MateMap mate,VertexIndexMap vm)1182 void maximum_weighted_matching(const Graph& g, MateMap mate, VertexIndexMap vm)
1183 {
1184 empty_matching< Graph, MateMap >::find_matching(g, mate);
1185 weighted_augmenting_path_finder< Graph, MateMap, VertexIndexMap > augmentor(
1186 g, mate, vm);
1187
1188 // can have |V| times augmenting at most
1189 for (std::size_t t = 0; t < num_vertices(g); ++t)
1190 {
1191 bool augmented = false;
1192 while (!augmented)
1193 {
1194 augmented = augmentor.augment_matching();
1195 if (!augmented)
1196 {
1197 // halt if adjusting dual variables can't bring potential
1198 // augment
1199 if (!augmentor.adjust_dual())
1200 break;
1201 }
1202 }
1203 if (!augmented)
1204 break;
1205 }
1206
1207 augmentor.get_current_matching(mate);
1208 }
1209
1210 template < typename Graph, typename MateMap >
maximum_weighted_matching(const Graph & g,MateMap mate)1211 inline void maximum_weighted_matching(const Graph& g, MateMap mate)
1212 {
1213 maximum_weighted_matching(g, mate, get(vertex_index, g));
1214 }
1215
1216 // brute-force matcher searches all possible combinations of matched edges to
1217 // get the maximum weighted matching which can be used for testing on small
1218 // graphs (within dozens vertices)
1219 template < typename Graph, typename MateMap, typename VertexIndexMap >
1220 class brute_force_matching
1221 {
1222 public:
1223 typedef
1224 typename graph_traits< Graph >::vertex_descriptor vertex_descriptor_t;
1225 typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator_t;
1226 typedef
1227 typename std::vector< vertex_descriptor_t >::iterator vertex_vec_iter_t;
1228 typedef typename graph_traits< Graph >::edge_iterator edge_iterator_t;
1229 typedef boost::iterator_property_map< vertex_vec_iter_t, VertexIndexMap >
1230 vertex_to_vertex_map_t;
1231
brute_force_matching(const Graph & arg_g,MateMap arg_mate,VertexIndexMap arg_vm)1232 brute_force_matching(
1233 const Graph& arg_g, MateMap arg_mate, VertexIndexMap arg_vm)
1234 : g(arg_g)
1235 , vm(arg_vm)
1236 , mate_vector(num_vertices(g))
1237 , best_mate_vector(num_vertices(g))
1238 , mate(mate_vector.begin(), vm)
1239 , best_mate(best_mate_vector.begin(), vm)
1240 {
1241 vertex_iterator_t vi, vi_end;
1242 for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
1243 best_mate[*vi] = mate[*vi] = get(arg_mate, *vi);
1244 }
1245
find_matching(PropertyMap pm)1246 template < typename PropertyMap > void find_matching(PropertyMap pm)
1247 {
1248 edge_iterator_t ei;
1249 boost::tie(ei, ei_end) = edges(g);
1250 select_edge(ei);
1251
1252 vertex_iterator_t vi, vi_end;
1253 for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
1254 put(pm, *vi, best_mate[*vi]);
1255 }
1256
1257 private:
1258 const Graph& g;
1259 VertexIndexMap vm;
1260 std::vector< vertex_descriptor_t > mate_vector, best_mate_vector;
1261 vertex_to_vertex_map_t mate, best_mate;
1262 edge_iterator_t ei_end;
1263
select_edge(edge_iterator_t ei)1264 void select_edge(edge_iterator_t ei)
1265 {
1266 if (ei == ei_end)
1267 {
1268 if (matching_weight_sum(g, mate)
1269 > matching_weight_sum(g, best_mate))
1270 {
1271 vertex_iterator_t vi, vi_end;
1272 for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
1273 best_mate[*vi] = mate[*vi];
1274 }
1275 return;
1276 }
1277
1278 vertex_descriptor_t v, w;
1279 v = source(*ei, g);
1280 w = target(*ei, g);
1281
1282 select_edge(++ei);
1283
1284 if (mate[v] == graph_traits< Graph >::null_vertex()
1285 && mate[w] == graph_traits< Graph >::null_vertex())
1286 {
1287 mate[v] = w;
1288 mate[w] = v;
1289 select_edge(ei);
1290 mate[v] = mate[w] = graph_traits< Graph >::null_vertex();
1291 }
1292 }
1293 };
1294
1295 template < typename Graph, typename MateMap, typename VertexIndexMap >
brute_force_maximum_weighted_matching(const Graph & g,MateMap mate,VertexIndexMap vm)1296 void brute_force_maximum_weighted_matching(
1297 const Graph& g, MateMap mate, VertexIndexMap vm)
1298 {
1299 empty_matching< Graph, MateMap >::find_matching(g, mate);
1300 brute_force_matching< Graph, MateMap, VertexIndexMap > brute_force_matcher(
1301 g, mate, vm);
1302 brute_force_matcher.find_matching(mate);
1303 }
1304
1305 template < typename Graph, typename MateMap >
brute_force_maximum_weighted_matching(const Graph & g,MateMap mate)1306 inline void brute_force_maximum_weighted_matching(const Graph& g, MateMap mate)
1307 {
1308 brute_force_maximum_weighted_matching(g, mate, get(vertex_index, g));
1309 }
1310
1311 }
1312
1313 #endif
1314