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