Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:37:33 

0001 //=======================================================================
0002 // Copyright (c) 2018 Yi Ji
0003 //
0004 // Distributed under the Boost Software License, Version 1.0.
0005 // (See accompanying file LICENSE_1_0.txt or copy at
0006 // http://www.boost.org/LICENSE_1_0.txt)
0007 //
0008 //=======================================================================
0009 
0010 #ifndef BOOST_GRAPH_MAXIMUM_WEIGHTED_MATCHING_HPP
0011 #define BOOST_GRAPH_MAXIMUM_WEIGHTED_MATCHING_HPP
0012 
0013 #include <algorithm> // for std::iter_swap
0014 #include <boost/shared_ptr.hpp>
0015 #include <boost/make_shared.hpp>
0016 #include <boost/graph/max_cardinality_matching.hpp>
0017 
0018 namespace boost
0019 {
0020 template < typename Graph, typename MateMap, typename VertexIndexMap >
0021 typename property_traits<
0022     typename property_map< Graph, edge_weight_t >::type >::value_type
0023 matching_weight_sum(const Graph& g, MateMap mate, VertexIndexMap vm)
0024 {
0025     typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator_t;
0026     typedef
0027         typename graph_traits< Graph >::vertex_descriptor vertex_descriptor_t;
0028     typedef typename property_traits< typename property_map< Graph,
0029         edge_weight_t >::type >::value_type edge_property_t;
0030 
0031     edge_property_t weight_sum = 0;
0032     vertex_iterator_t vi, vi_end;
0033 
0034     for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
0035     {
0036         vertex_descriptor_t v = *vi;
0037         if (get(mate, v) != graph_traits< Graph >::null_vertex()
0038             && get(vm, v) < get(vm, get(mate, v)))
0039             weight_sum += get(edge_weight, g, edge(v, mate[v], g).first);
0040     }
0041     return weight_sum;
0042 }
0043 
0044 template < typename Graph, typename MateMap >
0045 inline typename property_traits<
0046     typename property_map< Graph, edge_weight_t >::type >::value_type
0047 matching_weight_sum(const Graph& g, MateMap mate)
0048 {
0049     return matching_weight_sum(g, mate, get(vertex_index, g));
0050 }
0051 
0052 template < typename Graph, typename MateMap, typename VertexIndexMap >
0053 class weighted_augmenting_path_finder
0054 {
0055 public:
0056     template < typename T > struct map_vertex_to_
0057     {
0058         typedef boost::iterator_property_map<
0059             typename std::vector< T >::iterator, VertexIndexMap >
0060             type;
0061     };
0062     typedef typename graph::detail::VERTEX_STATE vertex_state_t;
0063     typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator_t;
0064     typedef
0065         typename graph_traits< Graph >::vertex_descriptor vertex_descriptor_t;
0066     typedef typename std::vector< vertex_descriptor_t >::const_iterator
0067         vertex_vec_iter_t;
0068     typedef
0069         typename graph_traits< Graph >::out_edge_iterator out_edge_iterator_t;
0070     typedef typename graph_traits< Graph >::edge_descriptor edge_descriptor_t;
0071     typedef typename graph_traits< Graph >::edge_iterator edge_iterator_t;
0072     typedef typename property_traits< typename property_map< Graph,
0073         edge_weight_t >::type >::value_type edge_property_t;
0074     typedef std::deque< vertex_descriptor_t > vertex_list_t;
0075     typedef std::vector< edge_descriptor_t > edge_list_t;
0076     typedef typename map_vertex_to_< vertex_descriptor_t >::type
0077         vertex_to_vertex_map_t;
0078     typedef
0079         typename map_vertex_to_< edge_property_t >::type vertex_to_weight_map_t;
0080     typedef typename map_vertex_to_< bool >::type vertex_to_bool_map_t;
0081     typedef typename map_vertex_to_< std::pair< vertex_descriptor_t,
0082         vertex_descriptor_t > >::type vertex_to_pair_map_t;
0083     typedef
0084         typename map_vertex_to_< std::pair< edge_descriptor_t, bool > >::type
0085             vertex_to_edge_map_t;
0086     typedef typename map_vertex_to_< vertex_to_edge_map_t >::type
0087         vertex_pair_to_edge_map_t;
0088 
0089     class blossom
0090     {
0091     public:
0092         typedef boost::shared_ptr< blossom > blossom_ptr_t;
0093         std::vector< blossom_ptr_t > sub_blossoms;
0094         edge_property_t dual_var;
0095         blossom_ptr_t father;
0096 
0097         blossom() : dual_var(0), father(blossom_ptr_t()) {}
0098 
0099         // get the base vertex of a blossom by recursively getting
0100         // its base sub-blossom, which is always the first one in
0101         // sub_blossoms because of how we create and maintain blossoms
0102         virtual vertex_descriptor_t get_base() const
0103         {
0104             const blossom* b = this;
0105             while (!b->sub_blossoms.empty())
0106                 b = b->sub_blossoms[0].get();
0107             return b->get_base();
0108         }
0109 
0110         // set a sub-blossom as a blossom's base by exchanging it
0111         // with its first sub-blossom
0112         void set_base(const blossom_ptr_t& sub)
0113         {
0114             for (blossom_iterator_t bi = sub_blossoms.begin();
0115                  bi != sub_blossoms.end(); ++bi)
0116             {
0117                 if (sub.get() == bi->get())
0118                 {
0119                     std::iter_swap(sub_blossoms.begin(), bi);
0120                     break;
0121                 }
0122             }
0123         }
0124 
0125         // get all vertices inside recursively
0126         virtual std::vector< vertex_descriptor_t > vertices() const
0127         {
0128             std::vector< vertex_descriptor_t > all_vertices;
0129             for (typename std::vector< blossom_ptr_t >::const_iterator bi
0130                  = sub_blossoms.begin();
0131                  bi != sub_blossoms.end(); ++bi)
0132             {
0133                 std::vector< vertex_descriptor_t > some_vertices
0134                     = (*bi)->vertices();
0135                 all_vertices.insert(all_vertices.end(), some_vertices.begin(),
0136                     some_vertices.end());
0137             }
0138             return all_vertices;
0139         }
0140     };
0141 
0142     // a trivial_blossom only has one vertex and no sub-blossom;
0143     // for each vertex v, in_blossom[v] is the trivial_blossom that contains it
0144     // directly
0145     class trivial_blossom : public blossom
0146     {
0147     public:
0148         trivial_blossom(vertex_descriptor_t v) : trivial_vertex(v) {}
0149         virtual vertex_descriptor_t get_base() const { return trivial_vertex; }
0150 
0151         virtual std::vector< vertex_descriptor_t > vertices() const
0152         {
0153             std::vector< vertex_descriptor_t > all_vertices;
0154             all_vertices.push_back(trivial_vertex);
0155             return all_vertices;
0156         }
0157 
0158     private:
0159         vertex_descriptor_t trivial_vertex;
0160     };
0161 
0162     typedef boost::shared_ptr< blossom > blossom_ptr_t;
0163     typedef typename std::vector< blossom_ptr_t >::iterator blossom_iterator_t;
0164     typedef
0165         typename map_vertex_to_< blossom_ptr_t >::type vertex_to_blossom_map_t;
0166 
0167     weighted_augmenting_path_finder(
0168         const Graph& arg_g, MateMap arg_mate, VertexIndexMap arg_vm)
0169     : g(arg_g)
0170     , vm(arg_vm)
0171     , null_edge(std::pair< edge_descriptor_t, bool >(
0172           num_edges(g) == 0 ? edge_descriptor_t() : *edges(g).first, false))
0173     , mate_vector(num_vertices(g))
0174     , label_S_vector(num_vertices(g), graph_traits< Graph >::null_vertex())
0175     , label_T_vector(num_vertices(g), graph_traits< Graph >::null_vertex())
0176     , outlet_vector(num_vertices(g), graph_traits< Graph >::null_vertex())
0177     , tau_idx_vector(num_vertices(g), graph_traits< Graph >::null_vertex())
0178     , dual_var_vector(std::vector< edge_property_t >(
0179           num_vertices(g), std::numeric_limits< edge_property_t >::min()))
0180     , pi_vector(std::vector< edge_property_t >(
0181           num_vertices(g), std::numeric_limits< edge_property_t >::max()))
0182     , gamma_vector(std::vector< edge_property_t >(
0183           num_vertices(g), std::numeric_limits< edge_property_t >::max()))
0184     , tau_vector(std::vector< edge_property_t >(
0185           num_vertices(g), std::numeric_limits< edge_property_t >::max()))
0186     , in_blossom_vector(num_vertices(g))
0187     , old_label_vector(num_vertices(g))
0188     , critical_edge_vectors(num_vertices(g),
0189           std::vector< std::pair< edge_descriptor_t, bool > >(
0190               num_vertices(g), null_edge))
0191     ,
0192 
0193         mate(mate_vector.begin(), vm)
0194     , label_S(label_S_vector.begin(), vm)
0195     , label_T(label_T_vector.begin(), vm)
0196     , outlet(outlet_vector.begin(), vm)
0197     , tau_idx(tau_idx_vector.begin(), vm)
0198     , dual_var(dual_var_vector.begin(), vm)
0199     , pi(pi_vector.begin(), vm)
0200     , gamma(gamma_vector.begin(), vm)
0201     , tau(tau_vector.begin(), vm)
0202     , in_blossom(in_blossom_vector.begin(), vm)
0203     , old_label(old_label_vector.begin(), vm)
0204     {
0205         vertex_iterator_t vi, vi_end;
0206         edge_iterator_t ei, ei_end;
0207 
0208         edge_property_t max_weight
0209             = std::numeric_limits< edge_property_t >::min();
0210         for (boost::tie(ei, ei_end) = edges(g); ei != ei_end; ++ei)
0211             max_weight = std::max(max_weight, get(edge_weight, g, *ei));
0212 
0213         typename std::vector<
0214             std::vector< std::pair< edge_descriptor_t, bool > > >::iterator vei;
0215 
0216         for (boost::tie(vi, vi_end) = vertices(g),
0217                             vei = critical_edge_vectors.begin();
0218              vi != vi_end; ++vi, ++vei)
0219         {
0220             vertex_descriptor_t u = *vi;
0221             mate[u] = get(arg_mate, u);
0222             dual_var[u] = 2 * max_weight;
0223             in_blossom[u] = boost::make_shared< trivial_blossom >(u);
0224             outlet[u] = u;
0225             critical_edge_vector.push_back(
0226                 vertex_to_edge_map_t(vei->begin(), vm));
0227         }
0228 
0229         critical_edge
0230             = vertex_pair_to_edge_map_t(critical_edge_vector.begin(), vm);
0231 
0232         init();
0233     }
0234 
0235     // return the top blossom where v is contained inside
0236     blossom_ptr_t in_top_blossom(vertex_descriptor_t v) const
0237     {
0238         blossom_ptr_t b = in_blossom[v];
0239         while (b->father != blossom_ptr_t())
0240             b = b->father;
0241         return b;
0242     }
0243 
0244     // check if vertex v is in blossom b
0245     bool is_in_blossom(blossom_ptr_t b, vertex_descriptor_t v) const
0246     {
0247         if (v == graph_traits< Graph >::null_vertex())
0248             return false;
0249         blossom_ptr_t vb = in_blossom[v]->father;
0250         while (vb != blossom_ptr_t())
0251         {
0252             if (vb.get() == b.get())
0253                 return true;
0254             vb = vb->father;
0255         }
0256         return false;
0257     }
0258 
0259     // return the base vertex of the top blossom that contains v
0260     inline vertex_descriptor_t base_vertex(vertex_descriptor_t v) const
0261     {
0262         return in_top_blossom(v)->get_base();
0263     }
0264 
0265     // add an existed top blossom of base vertex v into new top
0266     // blossom b as its sub-blossom
0267     void add_sub_blossom(blossom_ptr_t b, vertex_descriptor_t v)
0268     {
0269         blossom_ptr_t sub = in_top_blossom(v);
0270         sub->father = b;
0271         b->sub_blossoms.push_back(sub);
0272         if (sub->sub_blossoms.empty())
0273             return;
0274         for (blossom_iterator_t bi = top_blossoms.begin();
0275              bi != top_blossoms.end(); ++bi)
0276         {
0277             if (bi->get() == sub.get())
0278             {
0279                 top_blossoms.erase(bi);
0280                 break;
0281             }
0282         }
0283     }
0284 
0285     // when a top blossom is created or its base vertex getting an S-label,
0286     // add all edges incident to this blossom into even_edges
0287     void bloom(blossom_ptr_t b)
0288     {
0289         std::vector< vertex_descriptor_t > vertices_of_b = b->vertices();
0290         vertex_vec_iter_t vi;
0291         for (vi = vertices_of_b.begin(); vi != vertices_of_b.end(); ++vi)
0292         {
0293             out_edge_iterator_t oei, oei_end;
0294             for (boost::tie(oei, oei_end) = out_edges(*vi, g); oei != oei_end;
0295                  ++oei)
0296             {
0297                 if (target(*oei, g) != *vi && mate[*vi] != target(*oei, g))
0298                     even_edges.push_back(*oei);
0299             }
0300         }
0301     }
0302 
0303     // assigning a T-label to a non S-vertex, along with outlet and updating pi
0304     // value if updated pi[v] equals zero, augment the matching from its mate
0305     // vertex
0306     void put_T_label(vertex_descriptor_t v, vertex_descriptor_t T_label,
0307         vertex_descriptor_t outlet_v, edge_property_t pi_v)
0308     {
0309         if (label_S[v] != graph_traits< Graph >::null_vertex())
0310             return;
0311 
0312         label_T[v] = T_label;
0313         outlet[v] = outlet_v;
0314         pi[v] = pi_v;
0315 
0316         vertex_descriptor_t v_mate = mate[v];
0317         if (pi[v] == 0)
0318         {
0319             label_T[v_mate] = graph_traits< Graph >::null_vertex();
0320             label_S[v_mate] = v;
0321             bloom(in_top_blossom(v_mate));
0322         }
0323     }
0324 
0325     // get the missing T-label for a to-be-expanded base vertex
0326     // the missing T-label is the last vertex of the path from outlet[v] to v
0327     std::pair< vertex_descriptor_t, vertex_descriptor_t > missing_label(
0328         vertex_descriptor_t b_base)
0329     {
0330         vertex_descriptor_t missing_outlet = outlet[b_base];
0331 
0332         if (outlet[b_base] == b_base)
0333             return std::make_pair(
0334                 graph_traits< Graph >::null_vertex(), missing_outlet);
0335 
0336         vertex_iterator_t vi, vi_end;
0337         for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
0338             old_label[*vi] = std::make_pair(label_T[*vi], outlet[*vi]);
0339 
0340         std::pair< vertex_descriptor_t, vertex_state_t > child(
0341             outlet[b_base], graph::detail::V_EVEN);
0342         blossom_ptr_t b = in_blossom[child.first];
0343         for (; b->father->father != blossom_ptr_t(); b = b->father)
0344             ;
0345         child.first = b->get_base();
0346 
0347         if (child.first == b_base)
0348             return std::make_pair(
0349                 graph_traits< Graph >::null_vertex(), missing_outlet);
0350 
0351         while (true)
0352         {
0353             std::pair< vertex_descriptor_t, vertex_state_t > child_parent
0354                 = parent(child, true);
0355 
0356             for (b = in_blossom[child_parent.first];
0357                  b->father->father != blossom_ptr_t(); b = b->father)
0358                 ;
0359             missing_outlet = child_parent.first;
0360             child_parent.first = b->get_base();
0361 
0362             if (child_parent.first == b_base)
0363                 break;
0364             else
0365                 child = child_parent;
0366         }
0367         return std::make_pair(child.first, missing_outlet);
0368     }
0369 
0370     // expand a top blossom, put all its non-trivial sub-blossoms into
0371     // top_blossoms
0372     blossom_iterator_t expand_blossom(
0373         blossom_iterator_t bi, std::vector< blossom_ptr_t >& new_ones)
0374     {
0375         blossom_ptr_t b = *bi;
0376         for (blossom_iterator_t i = b->sub_blossoms.begin();
0377              i != b->sub_blossoms.end(); ++i)
0378         {
0379             blossom_ptr_t sub_blossom = *i;
0380             vertex_descriptor_t sub_base = sub_blossom->get_base();
0381             label_S[sub_base] = label_T[sub_base]
0382                 = graph_traits< Graph >::null_vertex();
0383             outlet[sub_base] = sub_base;
0384             sub_blossom->father = blossom_ptr_t();
0385             // new top blossoms cannot be pushed back into top_blossoms
0386             // immediately, because push_back() may cause reallocation and then
0387             // invalid iterators
0388             if (!sub_blossom->sub_blossoms.empty())
0389                 new_ones.push_back(sub_blossom);
0390         }
0391         return top_blossoms.erase(bi);
0392     }
0393 
0394     // when expanding a T-blossom with base v, it requires more operations:
0395     // supply the missing T-labels for new base vertices by picking the minimum
0396     // tau from vertices of each corresponding new top-blossoms; when label_T[v]
0397     // is null or we have a smaller tau from missing_label(v), replace T-label
0398     // and outlet of v (but don't bloom v)
0399     blossom_iterator_t expand_T_blossom(
0400         blossom_iterator_t bi, std::vector< blossom_ptr_t >& new_ones)
0401     {
0402         blossom_ptr_t b = *bi;
0403 
0404         vertex_descriptor_t b_base = b->get_base();
0405         std::pair< vertex_descriptor_t, vertex_descriptor_t > T_and_outlet
0406             = missing_label(b_base);
0407 
0408         blossom_iterator_t next_bi = expand_blossom(bi, new_ones);
0409 
0410         for (blossom_iterator_t i = b->sub_blossoms.begin();
0411              i != b->sub_blossoms.end(); ++i)
0412         {
0413             blossom_ptr_t sub_blossom = *i;
0414             vertex_descriptor_t sub_base = sub_blossom->get_base();
0415             vertex_descriptor_t min_tau_v
0416                 = graph_traits< Graph >::null_vertex();
0417             edge_property_t min_tau
0418                 = std::numeric_limits< edge_property_t >::max();
0419 
0420             std::vector< vertex_descriptor_t > sub_vertices
0421                 = sub_blossom->vertices();
0422             for (vertex_vec_iter_t v = sub_vertices.begin();
0423                  v != sub_vertices.end(); ++v)
0424             {
0425                 if (tau[*v] < min_tau)
0426                 {
0427                     min_tau = tau[*v];
0428                     min_tau_v = *v;
0429                 }
0430             }
0431 
0432             if (min_tau < std::numeric_limits< edge_property_t >::max())
0433                 put_T_label(
0434                     sub_base, tau_idx[min_tau_v], min_tau_v, tau[min_tau_v]);
0435         }
0436 
0437         if (label_T[b_base] == graph_traits< Graph >::null_vertex()
0438             || tau[old_label[b_base].second] < pi[b_base])
0439             boost::tie(label_T[b_base], outlet[b_base]) = T_and_outlet;
0440 
0441         return next_bi;
0442     }
0443 
0444     // when vertices v and w are matched to each other by augmenting,
0445     // we must set v/w as base vertex of any blossom who contains v/w and
0446     // is a sub-blossom of their lowest (smallest) common blossom
0447     void adjust_blossom(vertex_descriptor_t v, vertex_descriptor_t w)
0448     {
0449         blossom_ptr_t vb = in_blossom[v], wb = in_blossom[w],
0450                       lowest_common_blossom;
0451         std::vector< blossom_ptr_t > v_ancestors, w_ancestors;
0452 
0453         while (vb->father != blossom_ptr_t())
0454         {
0455             v_ancestors.push_back(vb->father);
0456             vb = vb->father;
0457         }
0458         while (wb->father != blossom_ptr_t())
0459         {
0460             w_ancestors.push_back(wb->father);
0461             wb = wb->father;
0462         }
0463 
0464         typename std::vector< blossom_ptr_t >::reverse_iterator i, j;
0465         i = v_ancestors.rbegin();
0466         j = w_ancestors.rbegin();
0467         while (i != v_ancestors.rend() && j != w_ancestors.rend()
0468             && i->get() == j->get())
0469         {
0470             lowest_common_blossom = *i;
0471             ++i;
0472             ++j;
0473         }
0474 
0475         vb = in_blossom[v];
0476         wb = in_blossom[w];
0477         while (vb->father != lowest_common_blossom)
0478         {
0479             vb->father->set_base(vb);
0480             vb = vb->father;
0481         }
0482         while (wb->father != lowest_common_blossom)
0483         {
0484             wb->father->set_base(wb);
0485             wb = wb->father;
0486         }
0487     }
0488 
0489     // every edge weight is multiplied by 4 to ensure integer weights
0490     // throughout the algorithm if all input weights are integers
0491     inline edge_property_t slack(const edge_descriptor_t& e) const
0492     {
0493         vertex_descriptor_t v, w;
0494         v = source(e, g);
0495         w = target(e, g);
0496         return dual_var[v] + dual_var[w] - 4 * get(edge_weight, g, e);
0497     }
0498 
0499     // backtrace one step on vertex v along the augmenting path
0500     // by its labels and its vertex state;
0501     // boolean parameter "use_old" means whether we are updating labels,
0502     // if we are, then we use old labels to backtrace and also we
0503     // don't jump to its base vertex when we reach an odd vertex
0504     std::pair< vertex_descriptor_t, vertex_state_t > parent(
0505         std::pair< vertex_descriptor_t, vertex_state_t > v,
0506         bool use_old = false) const
0507     {
0508         if (v.second == graph::detail::V_EVEN)
0509         {
0510             // a paranoid check: label_S shoule be the same as mate in
0511             // backtracing
0512             if (label_S[v.first] == graph_traits< Graph >::null_vertex())
0513                 label_S[v.first] = mate[v.first];
0514             return std::make_pair(label_S[v.first], graph::detail::V_ODD);
0515         }
0516         else if (v.second == graph::detail::V_ODD)
0517         {
0518             vertex_descriptor_t w = use_old ? old_label[v.first].first
0519                                             : base_vertex(label_T[v.first]);
0520             return std::make_pair(w, graph::detail::V_EVEN);
0521         }
0522         return std::make_pair(v.first, graph::detail::V_UNREACHED);
0523     }
0524 
0525     // backtrace from vertices v and w to their free (unmatched) ancesters,
0526     // return the nearest common ancestor (null_vertex if none) of v and w
0527     vertex_descriptor_t nearest_common_ancestor(vertex_descriptor_t v,
0528         vertex_descriptor_t w, vertex_descriptor_t& v_free_ancestor,
0529         vertex_descriptor_t& w_free_ancestor) const
0530     {
0531         std::pair< vertex_descriptor_t, vertex_state_t > v_up(
0532             v, graph::detail::V_EVEN);
0533         std::pair< vertex_descriptor_t, vertex_state_t > w_up(
0534             w, graph::detail::V_EVEN);
0535         vertex_descriptor_t nca;
0536         nca = w_free_ancestor = v_free_ancestor
0537             = graph_traits< Graph >::null_vertex();
0538 
0539         std::vector< bool > ancestor_of_w_vector(num_vertices(g), false);
0540         std::vector< bool > ancestor_of_v_vector(num_vertices(g), false);
0541         vertex_to_bool_map_t ancestor_of_w(ancestor_of_w_vector.begin(), vm);
0542         vertex_to_bool_map_t ancestor_of_v(ancestor_of_v_vector.begin(), vm);
0543 
0544         while (nca == graph_traits< Graph >::null_vertex()
0545             && (v_free_ancestor == graph_traits< Graph >::null_vertex()
0546                 || w_free_ancestor == graph_traits< Graph >::null_vertex()))
0547         {
0548             ancestor_of_w[w_up.first] = true;
0549             ancestor_of_v[v_up.first] = true;
0550 
0551             if (w_free_ancestor == graph_traits< Graph >::null_vertex())
0552                 w_up = parent(w_up);
0553             if (v_free_ancestor == graph_traits< Graph >::null_vertex())
0554                 v_up = parent(v_up);
0555 
0556             if (mate[v_up.first] == graph_traits< Graph >::null_vertex())
0557                 v_free_ancestor = v_up.first;
0558             if (mate[w_up.first] == graph_traits< Graph >::null_vertex())
0559                 w_free_ancestor = w_up.first;
0560 
0561             if (ancestor_of_w[v_up.first] == true || v_up.first == w_up.first)
0562                 nca = v_up.first;
0563             else if (ancestor_of_v[w_up.first] == true)
0564                 nca = w_up.first;
0565             else if (v_free_ancestor == w_free_ancestor
0566                 && v_free_ancestor != graph_traits< Graph >::null_vertex())
0567                 nca = v_up.first;
0568         }
0569 
0570         return nca;
0571     }
0572 
0573     // when a new top blossom b is created by connecting (v, w), we add
0574     // sub-blossoms into b along backtracing from v_prime and w_prime to
0575     // stop_vertex (the base vertex); also, we set labels and outlet for each
0576     // base vertex we pass by
0577     void make_blossom(blossom_ptr_t b, vertex_descriptor_t w_prime,
0578         vertex_descriptor_t v_prime, vertex_descriptor_t stop_vertex)
0579     {
0580         std::pair< vertex_descriptor_t, vertex_state_t > u(
0581             v_prime, graph::detail::V_ODD);
0582         std::pair< vertex_descriptor_t, vertex_state_t > u_up(
0583             w_prime, graph::detail::V_EVEN);
0584 
0585         for (; u_up.first != stop_vertex; u = u_up, u_up = parent(u))
0586         {
0587             if (u_up.second == graph::detail::V_EVEN)
0588             {
0589                 if (!in_top_blossom(u_up.first)->sub_blossoms.empty())
0590                     outlet[u_up.first] = label_T[u.first];
0591                 label_T[u_up.first] = outlet[u.first];
0592             }
0593             else if (u_up.second == graph::detail::V_ODD)
0594                 label_S[u_up.first] = u.first;
0595 
0596             add_sub_blossom(b, u_up.first);
0597         }
0598     }
0599 
0600     // the design of recursively expanding augmenting path in
0601     // (reversed_)retrieve_augmenting_path functions is inspired by same
0602     // functions in max_cardinality_matching.hpp; except that in weighted
0603     // matching, we use "outlet" vertices instead of "bridge" vertex pairs: if
0604     // blossom b is the smallest non-trivial blossom that contains its base
0605     // vertex v, then v and outlet[v] are where augmenting path enters and
0606     // leaves b
0607     void retrieve_augmenting_path(
0608         vertex_descriptor_t v, vertex_descriptor_t w, vertex_state_t v_state)
0609     {
0610         if (v == w)
0611             aug_path.push_back(v);
0612         else if (v_state == graph::detail::V_EVEN)
0613         {
0614             aug_path.push_back(v);
0615             retrieve_augmenting_path(label_S[v], w, graph::detail::V_ODD);
0616         }
0617         else if (v_state == graph::detail::V_ODD)
0618         {
0619             if (outlet[v] == v)
0620                 aug_path.push_back(v);
0621             else
0622                 reversed_retrieve_augmenting_path(
0623                     outlet[v], v, graph::detail::V_EVEN);
0624             retrieve_augmenting_path(label_T[v], w, graph::detail::V_EVEN);
0625         }
0626     }
0627 
0628     void reversed_retrieve_augmenting_path(
0629         vertex_descriptor_t v, vertex_descriptor_t w, vertex_state_t v_state)
0630     {
0631         if (v == w)
0632             aug_path.push_back(v);
0633         else if (v_state == graph::detail::V_EVEN)
0634         {
0635             reversed_retrieve_augmenting_path(
0636                 label_S[v], w, graph::detail::V_ODD);
0637             aug_path.push_back(v);
0638         }
0639         else if (v_state == graph::detail::V_ODD)
0640         {
0641             reversed_retrieve_augmenting_path(
0642                 label_T[v], w, graph::detail::V_EVEN);
0643             if (outlet[v] != v)
0644                 retrieve_augmenting_path(outlet[v], v, graph::detail::V_EVEN);
0645             else
0646                 aug_path.push_back(v);
0647         }
0648     }
0649 
0650     // correct labels for vertices in the augmenting path
0651     void relabel(vertex_descriptor_t v)
0652     {
0653         blossom_ptr_t b = in_blossom[v]->father;
0654 
0655         if (!is_in_blossom(b, mate[v]))
0656         { // if v is a new base vertex
0657             std::pair< vertex_descriptor_t, vertex_state_t > u(
0658                 v, graph::detail::V_EVEN);
0659             while (label_S[u.first] != u.first
0660                 && is_in_blossom(b, label_S[u.first]))
0661                 u = parent(u, true);
0662 
0663             vertex_descriptor_t old_base = u.first;
0664             if (label_S[old_base] != old_base)
0665             { // if old base is not exposed
0666                 label_T[v] = label_S[old_base];
0667                 outlet[v] = old_base;
0668             }
0669             else
0670             { // if old base is exposed then new label_T[v] is not in b,
0671                 // we must (i) make b2 the smallest blossom containing v but not
0672                 // as base vertex (ii) backtrace from b2's new base vertex to b
0673                 label_T[v] = graph_traits< Graph >::null_vertex();
0674                 for (b = b->father; b != blossom_ptr_t() && b->get_base() == v;
0675                      b = b->father)
0676                     ;
0677                 if (b != blossom_ptr_t())
0678                 {
0679                     u = std::make_pair(b->get_base(), graph::detail::V_ODD);
0680                     while (!is_in_blossom(
0681                         in_blossom[v]->father, old_label[u.first].first))
0682                         u = parent(u, true);
0683                     label_T[v] = u.first;
0684                     outlet[v] = old_label[u.first].first;
0685                 }
0686             }
0687         }
0688         else if (label_S[v] == v || !is_in_blossom(b, label_S[v]))
0689         { // if v is an old base vertex
0690             // let u be the new base vertex; backtrace from u's old T-label
0691             std::pair< vertex_descriptor_t, vertex_state_t > u(
0692                 b->get_base(), graph::detail::V_ODD);
0693             while (
0694                 old_label[u.first].first != graph_traits< Graph >::null_vertex()
0695                 && old_label[u.first].first != v)
0696                 u = parent(u, true);
0697             label_T[v] = old_label[u.first].second;
0698             outlet[v] = v;
0699         }
0700         else // if v is neither a new nor an old base vertex
0701             label_T[v] = label_S[v];
0702     }
0703 
0704     void augmenting(vertex_descriptor_t v, vertex_descriptor_t v_free_ancestor,
0705         vertex_descriptor_t w, vertex_descriptor_t w_free_ancestor)
0706     {
0707         vertex_iterator_t vi, vi_end;
0708 
0709         // retrieve the augmenting path and put it in aug_path
0710         reversed_retrieve_augmenting_path(
0711             v, v_free_ancestor, graph::detail::V_EVEN);
0712         retrieve_augmenting_path(w, w_free_ancestor, graph::detail::V_EVEN);
0713 
0714         // augment the matching along aug_path
0715         vertex_descriptor_t a, b;
0716         vertex_list_t reversed_aug_path;
0717         while (!aug_path.empty())
0718         {
0719             a = aug_path.front();
0720             aug_path.pop_front();
0721             reversed_aug_path.push_back(a);
0722             b = aug_path.front();
0723             aug_path.pop_front();
0724             reversed_aug_path.push_back(b);
0725 
0726             mate[a] = b;
0727             mate[b] = a;
0728 
0729             // reset base vertex for every blossom in augment path
0730             adjust_blossom(a, b);
0731         }
0732 
0733         for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
0734             old_label[*vi] = std::make_pair(label_T[*vi], outlet[*vi]);
0735 
0736         // correct labels for in-blossom vertices along aug_path
0737         while (!reversed_aug_path.empty())
0738         {
0739             a = reversed_aug_path.front();
0740             reversed_aug_path.pop_front();
0741 
0742             if (in_blossom[a]->father != blossom_ptr_t())
0743                 relabel(a);
0744         }
0745 
0746         for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
0747         {
0748             vertex_descriptor_t u = *vi;
0749             if (mate[u] != graph_traits< Graph >::null_vertex())
0750                 label_S[u] = mate[u];
0751         }
0752 
0753         // expand blossoms with zero dual variables
0754         std::vector< blossom_ptr_t > new_top_blossoms;
0755         for (blossom_iterator_t bi = top_blossoms.begin();
0756              bi != top_blossoms.end();)
0757         {
0758             if ((*bi)->dual_var <= 0)
0759                 bi = expand_blossom(bi, new_top_blossoms);
0760             else
0761                 ++bi;
0762         }
0763         top_blossoms.insert(top_blossoms.end(), new_top_blossoms.begin(),
0764             new_top_blossoms.end());
0765         init();
0766     }
0767 
0768     // create a new blossom and set labels for vertices inside
0769     void blossoming(vertex_descriptor_t v, vertex_descriptor_t v_prime,
0770         vertex_descriptor_t w, vertex_descriptor_t w_prime,
0771         vertex_descriptor_t nca)
0772     {
0773         vertex_iterator_t vi, vi_end;
0774 
0775         std::vector< bool > is_old_base_vector(num_vertices(g));
0776         vertex_to_bool_map_t is_old_base(is_old_base_vector.begin(), vm);
0777         for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
0778         {
0779             if (*vi == base_vertex(*vi))
0780                 is_old_base[*vi] = true;
0781         }
0782 
0783         blossom_ptr_t b = boost::make_shared< blossom >();
0784         add_sub_blossom(b, nca);
0785 
0786         label_T[w_prime] = v;
0787         label_T[v_prime] = w;
0788         outlet[w_prime] = w;
0789         outlet[v_prime] = v;
0790 
0791         make_blossom(b, w_prime, v_prime, nca);
0792         make_blossom(b, v_prime, w_prime, nca);
0793 
0794         label_T[nca] = graph_traits< Graph >::null_vertex();
0795         outlet[nca] = nca;
0796 
0797         top_blossoms.push_back(b);
0798         bloom(b);
0799 
0800         // set gamma[b_base] = min_slack{critical_edge(b_base, other_base)}
0801         // where each critical edge is updated before, by
0802         // argmin{slack(old_bases_in_b, other_base)};
0803         vertex_vec_iter_t i, j;
0804         std::vector< vertex_descriptor_t > b_vertices = b->vertices(),
0805                                            old_base_in_b, other_base;
0806         vertex_descriptor_t b_base = b->get_base();
0807         for (i = b_vertices.begin(); i != b_vertices.end(); ++i)
0808         {
0809             if (is_old_base[*i])
0810                 old_base_in_b.push_back(*i);
0811         }
0812         for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
0813         {
0814             if (*vi != b_base && *vi == base_vertex(*vi))
0815                 other_base.push_back(*vi);
0816         }
0817         for (i = other_base.begin(); i != other_base.end(); ++i)
0818         {
0819             edge_property_t min_slack
0820                 = std::numeric_limits< edge_property_t >::max();
0821             std::pair< edge_descriptor_t, bool > b_vi = null_edge;
0822             for (j = old_base_in_b.begin(); j != old_base_in_b.end(); ++j)
0823             {
0824                 if (critical_edge[*j][*i] != null_edge
0825                     && min_slack > slack(critical_edge[*j][*i].first))
0826                 {
0827                     min_slack = slack(critical_edge[*j][*i].first);
0828                     b_vi = critical_edge[*j][*i];
0829                 }
0830             }
0831             critical_edge[b_base][*i] = critical_edge[*i][b_base] = b_vi;
0832         }
0833         gamma[b_base] = std::numeric_limits< edge_property_t >::max();
0834         for (i = other_base.begin(); i != other_base.end(); ++i)
0835         {
0836             if (critical_edge[b_base][*i] != null_edge)
0837                 gamma[b_base] = std::min(
0838                     gamma[b_base], slack(critical_edge[b_base][*i].first));
0839         }
0840     }
0841 
0842     void init()
0843     {
0844         even_edges.clear();
0845 
0846         vertex_iterator_t vi, vi_end;
0847         typename std::vector<
0848             std::vector< std::pair< edge_descriptor_t, bool > > >::iterator vei;
0849 
0850         for (boost::tie(vi, vi_end) = vertices(g),
0851                             vei = critical_edge_vectors.begin();
0852              vi != vi_end; ++vi, ++vei)
0853         {
0854             vertex_descriptor_t u = *vi;
0855             out_edge_iterator_t ei, ei_end;
0856 
0857             gamma[u] = tau[u] = pi[u]
0858                 = std::numeric_limits< edge_property_t >::max();
0859             std::fill(vei->begin(), vei->end(), null_edge);
0860 
0861             if (base_vertex(u) != u)
0862                 continue;
0863 
0864             label_S[u] = label_T[u] = graph_traits< Graph >::null_vertex();
0865             outlet[u] = u;
0866 
0867             if (mate[u] == graph_traits< Graph >::null_vertex())
0868             {
0869                 label_S[u] = u;
0870                 bloom(in_top_blossom(u));
0871             }
0872         }
0873     }
0874 
0875     bool augment_matching()
0876     {
0877         vertex_descriptor_t v, w, w_free_ancestor, v_free_ancestor;
0878         v = w = w_free_ancestor = v_free_ancestor
0879             = graph_traits< Graph >::null_vertex();
0880         bool found_alternating_path = false;
0881 
0882         // note that we only use edges of zero slack value for augmenting
0883         while (!even_edges.empty() && !found_alternating_path)
0884         {
0885             // search for augmenting paths depth-first
0886             edge_descriptor_t current_edge = even_edges.back();
0887             even_edges.pop_back();
0888 
0889             v = source(current_edge, g);
0890             w = target(current_edge, g);
0891 
0892             vertex_descriptor_t v_prime = base_vertex(v);
0893             vertex_descriptor_t w_prime = base_vertex(w);
0894 
0895             // w_prime == v_prime implies that we get an edge that has been
0896             // shrunk into a blossom
0897             if (v_prime == w_prime)
0898                 continue;
0899 
0900             // a paranoid check
0901             if (label_S[v_prime] == graph_traits< Graph >::null_vertex())
0902             {
0903                 std::swap(v_prime, w_prime);
0904                 std::swap(v, w);
0905             }
0906 
0907             // w_prime may be unlabeled or have a T-label; replace the existed
0908             // T-label if the edge slack is smaller than current pi[w_prime] and
0909             // update it. Note that a T-label is "deserved" only when pi equals
0910             // zero. also update tau and tau_idx so that tau_idx becomes T-label
0911             // when a T-blossom is expanded
0912             if (label_S[w_prime] == graph_traits< Graph >::null_vertex())
0913             {
0914                 if (slack(current_edge) < pi[w_prime])
0915                     put_T_label(w_prime, v, w, slack(current_edge));
0916                 if (slack(current_edge) < tau[w])
0917                 {
0918                     if (in_blossom[w]->father == blossom_ptr_t()
0919                         || label_T[w_prime] == v
0920                         || label_T[w_prime]
0921                             == graph_traits< Graph >::null_vertex()
0922                         || nearest_common_ancestor(v_prime, label_T[w_prime],
0923                                v_free_ancestor, w_free_ancestor)
0924                             == graph_traits< Graph >::null_vertex())
0925                     {
0926                         tau[w] = slack(current_edge);
0927                         tau_idx[w] = v;
0928                     }
0929                 }
0930             }
0931 
0932             else
0933             {
0934                 if (slack(current_edge) > 0)
0935                 {
0936                     // update gamma and critical_edges when we have a smaller
0937                     // edge slack
0938                     gamma[v_prime]
0939                         = std::min(gamma[v_prime], slack(current_edge));
0940                     gamma[w_prime]
0941                         = std::min(gamma[w_prime], slack(current_edge));
0942                     if (critical_edge[v_prime][w_prime] == null_edge
0943                         || slack(critical_edge[v_prime][w_prime].first)
0944                             > slack(current_edge))
0945                     {
0946                         critical_edge[v_prime][w_prime]
0947                             = std::pair< edge_descriptor_t, bool >(
0948                                 current_edge, true);
0949                         critical_edge[w_prime][v_prime]
0950                             = std::pair< edge_descriptor_t, bool >(
0951                                 current_edge, true);
0952                     }
0953                     continue;
0954                 }
0955                 else if (slack(current_edge) == 0)
0956                 {
0957                     // if nca is null_vertex then we have an augmenting path;
0958                     // otherwise we have a new top blossom with nca as its base
0959                     // vertex
0960                     vertex_descriptor_t nca = nearest_common_ancestor(
0961                         v_prime, w_prime, v_free_ancestor, w_free_ancestor);
0962 
0963                     if (nca == graph_traits< Graph >::null_vertex())
0964                         found_alternating_path
0965                             = true; // to break out of the loop
0966                     else
0967                         blossoming(v, v_prime, w, w_prime, nca);
0968                 }
0969             }
0970         }
0971 
0972         if (!found_alternating_path)
0973             return false;
0974 
0975         augmenting(v, v_free_ancestor, w, w_free_ancestor);
0976         return true;
0977     }
0978 
0979     // slack the vertex and blossom dual variables when there is no augmenting
0980     // path found according to the primal-dual method
0981     bool adjust_dual()
0982     {
0983         edge_property_t delta1, delta2, delta3, delta4, delta;
0984         delta1 = delta2 = delta3 = delta4
0985             = std::numeric_limits< edge_property_t >::max();
0986 
0987         vertex_iterator_t vi, vi_end;
0988 
0989         for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
0990         {
0991             delta1 = std::min(delta1, dual_var[*vi]);
0992             delta4 = pi[*vi] > 0 ? std::min(delta4, pi[*vi]) : delta4;
0993             if (*vi == base_vertex(*vi))
0994                 delta3 = std::min(delta3, gamma[*vi] / 2);
0995         }
0996 
0997         for (blossom_iterator_t bi = top_blossoms.begin();
0998              bi != top_blossoms.end(); ++bi)
0999         {
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 
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 >
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 >
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 
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 
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 
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 >
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 >
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
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!