//=======================================================================

 // Copyright (c) 2005 Aaron Windsor

 //

 // Distributed under the Boost Software License, Version 1.0. 

 // (See accompanying file LICENSE_1_0.txt or copy at

 // http://www.boost.org/LICENSE_1_0.txt)

 //

 //=======================================================================

 #ifndef BOOST_GRAPH_MAXIMUM_CARDINALITY_MATCHING_HPP

 #define BOOST_GRAPH_MAXIMUM_CARDINALITY_MATCHING_HPP

 #include <vector>

 #include <list>

 #include <deque>

 #include <algorithm>                     // for std::sort and std::stable_sort

 #include <utility>                       // for std::pair

 #include <boost/property_map.hpp>

 #include <boost/utility.hpp>             // for boost::tie

 #include <boost/graph/graph_traits.hpp>  

 #include <boost/graph/visitors.hpp>

 #include <boost/graph/depth_first_search.hpp>

 #include <boost/graph/filtered_graph.hpp>

 #include <boost/pending/disjoint_sets.hpp>

 #include <boost/assert.hpp>

 namespace boost

 {

   namespace graph { namespace detail {

     enum { V_EVEN, V_ODD, V_UNREACHED };

   } } // end namespace graph::detail

   template <typename Graph, typename MateMap, typename VertexIndexMap>

   typename graph_traits<Graph>::vertices_size_type 

   matching_size(const Graph& g, MateMap mate, VertexIndexMap vm)

   {

     typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator_t;

     typedef typename graph_traits<Graph>::vertex_descriptor

       vertex_descriptor_t;

     typedef typename graph_traits<Graph>::vertices_size_type v_size_t;

     v_size_t size_of_matching = 0;

     vertex_iterator_t vi, vi_end;

     for(tie(vi,vi_end) = vertices(g); vi != vi_end; ++vi)

       {

     vertex_descriptor_t v = *vi;

     if (get(mate,v) != graph_traits<Graph>::null_vertex() 

             && get(vm,v) < get(vm,get(mate,v)))

       ++size_of_matching;

       }

     return size_of_matching;

   }

   template <typename Graph, typename MateMap>

   inline typename graph_traits<Graph>::vertices_size_type

   matching_size(const Graph& g, MateMap mate)

   {

     return matching_size(g, mate, get(vertex_index,g));

   }

   template <typename Graph, typename MateMap, typename VertexIndexMap>

   bool is_a_matching(const Graph& g, MateMap mate, VertexIndexMap vm)

   {

     typedef typename graph_traits<Graph>::vertex_descriptor

       vertex_descriptor_t;

     typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator_t;

     vertex_iterator_t vi, vi_end;

     for( tie(vi,vi_end) = vertices(g); vi != vi_end; ++vi)

       {

     vertex_descriptor_t v = *vi;

     if (get(mate,v) != graph_traits<Graph>::null_vertex() 

             && v != get(mate,get(mate,v)))

       return false;

       }    

     return true;

   }

   template <typename Graph, typename MateMap>

   inline bool is_a_matching(const Graph& g, MateMap mate)

   {

     return is_a_matching(g, mate, get(vertex_index,g));

   }

   //***************************************************************************

   //***************************************************************************

   //               Maximum Cardinality Matching Functors 

   //***************************************************************************

   //***************************************************************************

   template <typename Graph, typename MateMap, 

             typename VertexIndexMap = dummy_property_map>

   struct no_augmenting_path_finder

   {

     no_augmenting_path_finder(const Graph& g, MateMap mate, VertexIndexMap vm)

     { }

     inline bool augment_matching() { return false; }

     template <typename PropertyMap>

     void get_current_matching(PropertyMap p) {}

   };

   template <typename Graph, typename MateMap, typename VertexIndexMap>

   class edmonds_augmenting_path_finder

   {

     // This implementation of Edmonds' matching algorithm closely

     // follows Tarjan's description of the algorithm in "Data

     // Structures and Network Algorithms."

   public:

     //generates the type of an iterator property map from vertices to type X

     template <typename X>

     struct map_vertex_to_ 

     { 

       typedef boost::iterator_property_map<typename std::vector<X>::iterator,

                                            VertexIndexMap> type; 

     };

     typedef typename graph_traits<Graph>::vertex_descriptor

       vertex_descriptor_t;

     typedef typename std::pair< vertex_descriptor_t, vertex_descriptor_t >

       vertex_pair_t;

     typedef typename graph_traits<Graph>::edge_descriptor edge_descriptor_t; 

     typedef typename graph_traits<Graph>::vertices_size_type v_size_t;

     typedef typename graph_traits<Graph>::edges_size_type e_size_t;

     typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator_t;

     typedef typename graph_traits<Graph>::out_edge_iterator 

       out_edge_iterator_t;

     typedef typename std::deque<vertex_descriptor_t> vertex_list_t;

     typedef typename std::vector<edge_descriptor_t> edge_list_t;

     typedef typename map_vertex_to_<vertex_descriptor_t>::type 

       vertex_to_vertex_map_t;

     typedef typename map_vertex_to_<int>::type vertex_to_int_map_t;

     typedef typename map_vertex_to_<vertex_pair_t>::type 

       vertex_to_vertex_pair_map_t;

     typedef typename map_vertex_to_<v_size_t>::type vertex_to_vsize_map_t;

     typedef typename map_vertex_to_<e_size_t>::type vertex_to_esize_map_t;

     edmonds_augmenting_path_finder(const Graph& arg_g, MateMap arg_mate, 

                                    VertexIndexMap arg_vm) : 

       g(arg_g),

       vm(arg_vm),

       n_vertices(num_vertices(arg_g)),

       mate_vector(n_vertices),

       ancestor_of_v_vector(n_vertices),

       ancestor_of_w_vector(n_vertices),

       vertex_state_vector(n_vertices),

       origin_vector(n_vertices),

       pred_vector(n_vertices),

       bridge_vector(n_vertices),

       ds_parent_vector(n_vertices),

       ds_rank_vector(n_vertices),

       mate(mate_vector.begin(), vm),

       ancestor_of_v(ancestor_of_v_vector.begin(), vm),

       ancestor_of_w(ancestor_of_w_vector.begin(), vm),

       vertex_state(vertex_state_vector.begin(), vm),

       origin(origin_vector.begin(), vm),

       pred(pred_vector.begin(), vm),

       bridge(bridge_vector.begin(), vm),

       ds_parent_map(ds_parent_vector.begin(), vm),

       ds_rank_map(ds_rank_vector.begin(), vm),

       ds(ds_rank_map, ds_parent_map)

     {

       vertex_iterator_t vi, vi_end;

       for(tie(vi,vi_end) = vertices(g); vi != vi_end; ++vi)

     mate[*vi] = get(arg_mate, *vi);

     }

     bool augment_matching()

     {

       //As an optimization, some of these values can be saved from one

       //iteration to the next instead of being re-initialized each

       //iteration, allowing for "lazy blossom expansion." This is not

       //currently implemented.

       e_size_t timestamp = 0;

       even_edges.clear();

       vertex_iterator_t vi, vi_end;

       for(tie(vi,vi_end) = vertices(g); vi != vi_end; ++vi)

     {

       vertex_descriptor_t u = *vi;

       origin[u] = u;

       pred[u] = u;

       ancestor_of_v[u] = 0;

       ancestor_of_w[u] = 0;

       ds.make_set(u);

       if (mate[u] == graph_traits<Graph>::null_vertex())

         {

           vertex_state[u] = graph::detail::V_EVEN;

           out_edge_iterator_t ei, ei_end;

           for(tie(ei,ei_end) = out_edges(u,g); ei != ei_end; ++ei)

         even_edges.push_back( *ei );

         }

       else

         vertex_state[u] = graph::detail::V_UNREACHED;      

     }

       //end initializations

       vertex_descriptor_t v,w,w_free_ancestor,v_free_ancestor;

       w_free_ancestor = graph_traits<Graph>::null_vertex();

       v_free_ancestor = graph_traits<Graph>::null_vertex(); 

       bool found_alternating_path = false;

       while(!even_edges.empty() && !found_alternating_path)

     {

       // since we push even edges onto the back of the list as

       // they're discovered, taking them off the back will search

       // for augmenting paths depth-first.

       edge_descriptor_t current_edge = even_edges.back();

       even_edges.pop_back();

       v = source(current_edge,g);

       w = target(current_edge,g);

       vertex_descriptor_t v_prime = origin[ds.find_set(v)];

       vertex_descriptor_t w_prime = origin[ds.find_set(w)];

       // because of the way we put all of the edges on the queue,

       // v_prime should be labeled V_EVEN; the following is a

       // little paranoid but it could happen...

       if (vertex_state[v_prime] != graph::detail::V_EVEN)

         {

           std::swap(v_prime,w_prime);

           std::swap(v,w);

         }

       if (vertex_state[w_prime] == graph::detail::V_UNREACHED)

         {

           vertex_state[w_prime] = graph::detail::V_ODD;

           vertex_state[mate[w_prime]] = graph::detail::V_EVEN;

           out_edge_iterator_t ei, ei_end;

           for( tie(ei,ei_end) = out_edges(mate[w_prime], g); ei != ei_end; ++ei)

         even_edges.push_back(*ei);

           pred[w_prime] = v;

         }

           //w_prime == v_prime can happen below if we get an edge that has been

           //shrunk into a blossom

       else if (vertex_state[w_prime] == graph::detail::V_EVEN && w_prime != v_prime) 

         {                                                             

           vertex_descriptor_t w_up = w_prime;

           vertex_descriptor_t v_up = v_prime;

           vertex_descriptor_t nearest_common_ancestor 

                 = graph_traits<Graph>::null_vertex();

           w_free_ancestor = graph_traits<Graph>::null_vertex();

           v_free_ancestor = graph_traits<Graph>::null_vertex();

           // We now need to distinguish between the case that

           // w_prime and v_prime share an ancestor under the

           // "parent" relation, in which case we've found a

           // blossom and should shrink it, or the case that

           // w_prime and v_prime both have distinct ancestors that

           // are free, in which case we've found an alternating

           // path between those two ancestors.

           ++timestamp;

           while (nearest_common_ancestor == graph_traits<Graph>::null_vertex() && 

              (v_free_ancestor == graph_traits<Graph>::null_vertex() || 

               w_free_ancestor == graph_traits<Graph>::null_vertex()

               )

              )

         {

           ancestor_of_w[w_up] = timestamp;

           ancestor_of_v[v_up] = timestamp;

           if (w_free_ancestor == graph_traits<Graph>::null_vertex())

             w_up = parent(w_up);

           if (v_free_ancestor == graph_traits<Graph>::null_vertex())

             v_up = parent(v_up);

           if (mate[v_up] == graph_traits<Graph>::null_vertex())

             v_free_ancestor = v_up;

           if (mate[w_up] == graph_traits<Graph>::null_vertex())

             w_free_ancestor = w_up;

           if (ancestor_of_w[v_up] == timestamp)

             nearest_common_ancestor = v_up;

           else if (ancestor_of_v[w_up] == timestamp)

             nearest_common_ancestor = w_up;

           else if (v_free_ancestor == w_free_ancestor && 

                v_free_ancestor != graph_traits<Graph>::null_vertex())

             nearest_common_ancestor = v_up;

         }

           if (nearest_common_ancestor == graph_traits<Graph>::null_vertex())

         found_alternating_path = true; //to break out of the loop

           else

         {

           //shrink the blossom

           link_and_set_bridges(w_prime, nearest_common_ancestor, std::make_pair(w,v));

           link_and_set_bridges(v_prime, nearest_common_ancestor, std::make_pair(v,w));

         }

         }      

     }

       if (!found_alternating_path)

     return false;

       // retrieve the augmenting path and put it in aug_path

       reversed_retrieve_augmenting_path(v, v_free_ancestor);

       retrieve_augmenting_path(w, w_free_ancestor);

       // augment the matching along aug_path

       vertex_descriptor_t a,b;

       while (!aug_path.empty())

     {

       a = aug_path.front();

       aug_path.pop_front();

       b = aug_path.front();

       aug_path.pop_front();

       mate[a] = b;

       mate[b] = a;

     }

       return true;

     }

     template <typename PropertyMap>

     void get_current_matching(PropertyMap pm)

     {

       vertex_iterator_t vi,vi_end;

       for(tie(vi,vi_end) = vertices(g); vi != vi_end; ++vi)

     put(pm, *vi, mate[*vi]);

     }

     template <typename PropertyMap>

     void get_vertex_state_map(PropertyMap pm)

     {

       vertex_iterator_t vi,vi_end;

       for(tie(vi,vi_end) = vertices(g); vi != vi_end; ++vi)

     put(pm, *vi, vertex_state[origin[ds.find_set(*vi)]]);

     }

   private:    

     vertex_descriptor_t parent(vertex_descriptor_t x)

     {

       if (vertex_state[x] == graph::detail::V_EVEN 

           && mate[x] != graph_traits<Graph>::null_vertex())

     return mate[x];

       else if (vertex_state[x] == graph::detail::V_ODD)

     return origin[ds.find_set(pred[x])];

       else

     return x;

     }

     void link_and_set_bridges(vertex_descriptor_t x, 

                               vertex_descriptor_t stop_vertex, 

                   vertex_pair_t the_bridge)

     {

       for(vertex_descriptor_t v = x; v != stop_vertex; v = parent(v))

     {

       ds.union_set(v, stop_vertex);

       origin[ds.find_set(stop_vertex)] = stop_vertex;

       if (vertex_state[v] == graph::detail::V_ODD)

         {

           bridge[v] = the_bridge;

           out_edge_iterator_t oei, oei_end;

           for(tie(oei, oei_end) = out_edges(v,g); oei != oei_end; ++oei)

         even_edges.push_back(*oei);

         }

     }

     }

     // Since none of the STL containers support both constant-time

     // concatenation and reversal, the process of expanding an

     // augmenting path once we know one exists is a little more

     // complicated than it has to be. If we know the path is from v to

     // w, then the augmenting path is recursively defined as:

     //

     // path(v,w) = [v], if v = w

     //           = concat([v, mate[v]], path(pred[mate[v]], w), 

     //                if v != w and vertex_state[v] == graph::detail::V_EVEN

     //           = concat([v], reverse(path(x,mate[v])), path(y,w)),

     //                if v != w, vertex_state[v] == graph::detail::V_ODD, and bridge[v] = (x,y)

     //

     // These next two mutually recursive functions implement this definition.

     void retrieve_augmenting_path(vertex_descriptor_t v, vertex_descriptor_t w)  

     {

       if (v == w)

     aug_path.push_back(v);

       else if (vertex_state[v] == graph::detail::V_EVEN)

     {

       aug_path.push_back(v);

       aug_path.push_back(mate[v]);

       retrieve_augmenting_path(pred[mate[v]], w);

     }

       else //vertex_state[v] == graph::detail::V_ODD 

     {

       aug_path.push_back(v);

       reversed_retrieve_augmenting_path(bridge[v].first, mate[v]);

       retrieve_augmenting_path(bridge[v].second, w);

     }

     }

     void reversed_retrieve_augmenting_path(vertex_descriptor_t v,

                                            vertex_descriptor_t w)  

     {

       if (v == w)

     aug_path.push_back(v);

       else if (vertex_state[v] == graph::detail::V_EVEN)

     {

       reversed_retrieve_augmenting_path(pred[mate[v]], w);

       aug_path.push_back(mate[v]);

       aug_path.push_back(v);

     }

       else //vertex_state[v] == graph::detail::V_ODD 

     {

       reversed_retrieve_augmenting_path(bridge[v].second, w);

       retrieve_augmenting_path(bridge[v].first, mate[v]);

       aug_path.push_back(v);

     }

     }

     //private data members

     const Graph& g;

     VertexIndexMap vm;

     v_size_t n_vertices;

     //storage for the property maps below

     std::vector<vertex_descriptor_t> mate_vector;

     std::vector<e_size_t> ancestor_of_v_vector;

     std::vector<e_size_t> ancestor_of_w_vector;

     std::vector<int> vertex_state_vector;

     std::vector<vertex_descriptor_t> origin_vector;

     std::vector<vertex_descriptor_t> pred_vector;

     std::vector<vertex_pair_t> bridge_vector;

     std::vector<vertex_descriptor_t> ds_parent_vector;

     std::vector<v_size_t> ds_rank_vector;

     //iterator property maps

     vertex_to_vertex_map_t mate;

     vertex_to_esize_map_t ancestor_of_v;

     vertex_to_esize_map_t ancestor_of_w;

     vertex_to_int_map_t vertex_state;

     vertex_to_vertex_map_t origin;

     vertex_to_vertex_map_t pred;

     vertex_to_vertex_pair_map_t bridge;

     vertex_to_vertex_map_t ds_parent_map;

     vertex_to_vsize_map_t ds_rank_map;

     vertex_list_t aug_path;

     edge_list_t even_edges;

     disjoint_sets< vertex_to_vsize_map_t, vertex_to_vertex_map_t > ds;

   };

   //***************************************************************************

   //***************************************************************************

   //               Initial Matching Functors

   //***************************************************************************

   //***************************************************************************

   template <typename Graph, typename MateMap>

   struct greedy_matching

   {

     typedef typename graph_traits< Graph >::vertex_descriptor vertex_descriptor_t;

     typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator_t;

     typedef typename graph_traits< Graph >::edge_descriptor edge_descriptor_t; 

     typedef typename graph_traits< Graph >::edge_iterator edge_iterator_t;

     static void find_matching(const Graph& g, MateMap mate)

     {

       vertex_iterator_t vi, vi_end;

       for(tie(vi,vi_end) = vertices(g); vi != vi_end; ++vi)

     put(mate, *vi, graph_traits<Graph>::null_vertex());

       edge_iterator_t ei, ei_end;

       for( tie(ei, ei_end) = edges(g); ei != ei_end; ++ei)

     {

       edge_descriptor_t e = *ei;

       vertex_descriptor_t u = source(e,g);

       vertex_descriptor_t v = target(e,g);

       if (get(mate,u) == get(mate,v))  

         //only way equality can hold is if

             //   mate[u] == mate[v] == null_vertex

         {

           put(mate,u,v);

           put(mate,v,u);

         }

     }    

     }

   };

   template <typename Graph, typename MateMap>

   struct extra_greedy_matching

   {

     // The "extra greedy matching" is formed by repeating the

     // following procedure as many times as possible: Choose the

     // unmatched vertex v of minimum non-zero degree.  Choose the

     // neighbor w of v which is unmatched and has minimum degree over

     // all of v's neighbors. Add (u,v) to the matching. Ties for

     // either choice are broken arbitrarily. This procedure takes time

     // O(m log n), where m is the number of edges in the graph and n

     // is the number of vertices.

     typedef typename graph_traits< Graph >::vertex_descriptor

       vertex_descriptor_t;

     typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator_t;

     typedef typename graph_traits< Graph >::edge_descriptor edge_descriptor_t; 

     typedef typename graph_traits< Graph >::edge_iterator edge_iterator_t;

     typedef std::pair<vertex_descriptor_t, vertex_descriptor_t> vertex_pair_t;

     struct select_first

     {

       inline static vertex_descriptor_t select_vertex(const vertex_pair_t p) 

       {return p.first;}

     };

     struct select_second

     {

       inline static vertex_descriptor_t select_vertex(const vertex_pair_t p) 

       {return p.second;}

     };

     template <class PairSelector>

     class less_than_by_degree

     {

     public:

       less_than_by_degree(const Graph& g): m_g(g) {}

       bool operator() (const vertex_pair_t x, const vertex_pair_t y)

       {

     return 

       out_degree(PairSelector::select_vertex(x), m_g) 

       < out_degree(PairSelector::select_vertex(y), m_g);

       }

     private:

       const Graph& m_g;

     };

     static void find_matching(const Graph& g, MateMap mate)

     {

       typedef std::vector<std::pair<vertex_descriptor_t, vertex_descriptor_t> >

         directed_edges_vector_t;

       directed_edges_vector_t edge_list;

       vertex_iterator_t vi, vi_end;

       for(tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)

     put(mate, *vi, graph_traits<Graph>::null_vertex());

       edge_iterator_t ei, ei_end;

       for(tie(ei, ei_end) = edges(g); ei != ei_end; ++ei)

     {

       edge_descriptor_t e = *ei;

       vertex_descriptor_t u = source(e,g);

       vertex_descriptor_t v = target(e,g);

       edge_list.push_back(std::make_pair(u,v));

       edge_list.push_back(std::make_pair(v,u));

     }

       //sort the edges by the degree of the target, then (using a

       //stable sort) by degree of the source

       std::sort(edge_list.begin(), edge_list.end(), 

                 less_than_by_degree<select_second>(g));

       std::stable_sort(edge_list.begin(), edge_list.end(), 

                        less_than_by_degree<select_first>(g));

       //construct the extra greedy matching

       for(typename directed_edges_vector_t::const_iterator itr = edge_list.begin(); itr != edge_list.end(); ++itr)

     {

       if (get(mate,itr->first) == get(mate,itr->second)) 

         //only way equality can hold is if mate[itr->first] == mate[itr->second] == null_vertex

         {

           put(mate, itr->first, itr->second);

           put(mate, itr->second, itr->first);

         }

     }    

     }

   };

   template <typename Graph, typename MateMap>

   struct empty_matching

   { 

     typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator_t;

     static void find_matching(const Graph& g, MateMap mate)

     {

       vertex_iterator_t vi, vi_end;

       for(tie(vi,vi_end) = vertices(g); vi != vi_end; ++vi)

     put(mate, *vi, graph_traits<Graph>::null_vertex());

     }

   };

   //***************************************************************************

   //***************************************************************************

   //               Matching Verifiers

   //***************************************************************************

   //***************************************************************************

   namespace detail

   {

     template <typename SizeType>

     class odd_components_counter : public dfs_visitor<>

     // This depth-first search visitor will count the number of connected 

     // components with an odd number of vertices. It's used by 

     // maximum_matching_verifier.

     {

     public:

       odd_components_counter(SizeType& c_count):

     m_count(c_count)

       {

     m_count = 0;

       }

       template <class Vertex, class Graph>

       void start_vertex(Vertex v, Graph&) 

       {  

     addend = -1; 

       }

       template <class Vertex, class Graph>

       void discover_vertex(Vertex u, Graph&) 

       {

     addend *= -1;

     m_count += addend;

       }

     protected:

       SizeType& m_count;

     private:

       SizeType addend;

     };

   }//namespace detail

   template <typename Graph, typename MateMap, 

             typename VertexIndexMap = dummy_property_map>

   struct no_matching_verifier

   {

     inline static bool 

     verify_matching(const Graph& g, MateMap mate, VertexIndexMap vm) 

     { return true;}

   };

   template <typename Graph, typename MateMap, typename VertexIndexMap>

   struct maximum_cardinality_matching_verifier

   {

     template <typename X>

     struct map_vertex_to_

     { 

       typedef boost::iterator_property_map<typename std::vector<X>::iterator,

                                            VertexIndexMap> type; 

     };

     typedef typename graph_traits<Graph>::vertex_descriptor 

       vertex_descriptor_t;

     typedef typename graph_traits<Graph>::vertices_size_type v_size_t;

     typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator_t;

     typedef typename map_vertex_to_<int>::type vertex_to_int_map_t;

     typedef typename map_vertex_to_<vertex_descriptor_t>::type 

       vertex_to_vertex_map_t;

     template <typename VertexStateMap>

     struct non_odd_vertex {

       //this predicate is used to create a filtered graph that

       //excludes vertices labeled "graph::detail::V_ODD"

       non_odd_vertex() : vertex_state(0) { }

       non_odd_vertex(VertexStateMap* arg_vertex_state) 

         : vertex_state(arg_vertex_state) { }

       template <typename Vertex>

       bool operator()(const Vertex& v) const {

     BOOST_ASSERT(vertex_state);

     return get(*vertex_state, v) != graph::detail::V_ODD;

       }

       VertexStateMap* vertex_state;

     };

     static bool verify_matching(const Graph& g, MateMap mate, VertexIndexMap vm)

     {

       //For any graph G, let o(G) be the number of connected

       //components in G of odd size. For a subset S of G's vertex set

       //V(G), let (G - S) represent the subgraph of G induced by

       //removing all vertices in S from G. Let M(G) be the size of the

       //maximum cardinality matching in G. Then the Tutte-Berge

       //formula guarantees that

       //

       //           2 * M(G) = min ( |V(G)| + |U| + o(G - U) )

       //

       //where the minimum is taken over all subsets U of

       //V(G). Edmonds' algorithm finds a set U that achieves the

       //minimum in the above formula, namely the vertices labeled

       //"ODD." This function runs one iteration of Edmonds' algorithm

       //to find U, then verifies that the size of the matching given

       //by mate satisfies the Tutte-Berge formula.

       //first, make sure it's a valid matching

       if (!is_a_matching(g,mate,vm))

       return false;

       //We'll try to augment the matching once. This serves two

       //purposes: first, if we find some augmenting path, the matching

       //is obviously non-maximum. Second, running edmonds' algorithm

       //on a graph with no augmenting path will create the

       //Edmonds-Gallai decomposition that we need as a certificate of

       //maximality - we can get it by looking at the vertex_state map

       //that results.

       edmonds_augmenting_path_finder<Graph,MateMap,VertexIndexMap>

         augmentor(g,mate,vm);

       if (augmentor.augment_matching())

       return false;

       std::vector<int> vertex_state_vector(num_vertices(g));

       vertex_to_int_map_t vertex_state(vertex_state_vector.begin(), vm);

       augmentor.get_vertex_state_map(vertex_state);

       //count the number of graph::detail::V_ODD vertices

       v_size_t num_odd_vertices = 0;

       vertex_iterator_t vi, vi_end;

       for(tie(vi,vi_end) = vertices(g); vi != vi_end; ++vi)

     if (vertex_state[*vi] == graph::detail::V_ODD)

       ++num_odd_vertices;

       //count the number of connected components with odd cardinality

       //in the graph without graph::detail::V_ODD vertices

       non_odd_vertex<vertex_to_int_map_t> filter(&vertex_state);

       filtered_graph<Graph, keep_all, non_odd_vertex<vertex_to_int_map_t> > fg(g, keep_all(), filter);

       v_size_t num_odd_components;

       detail::odd_components_counter<v_size_t> occ(num_odd_components);

       depth_first_search(fg, visitor(occ).vertex_index_map(vm));

       if (2 * matching_size(g,mate,vm) == num_vertices(g) + num_odd_vertices - num_odd_components)

     return true;

       else

     return false;

     }

   };

   template <typename Graph, 

         typename MateMap,

         typename VertexIndexMap,

         template <typename, typename, typename> class AugmentingPathFinder, 

         template <typename, typename> class InitialMatchingFinder,

         template <typename, typename, typename> class MatchingVerifier>

   bool matching(const Graph& g, MateMap mate, VertexIndexMap vm)

   {

     InitialMatchingFinder<Graph,MateMap>::find_matching(g,mate);

     AugmentingPathFinder<Graph,MateMap,VertexIndexMap> augmentor(g,mate,vm);

     bool not_maximum_yet = true;

     while(not_maximum_yet)

       {

     not_maximum_yet = augmentor.augment_matching();

       }

     augmentor.get_current_matching(mate);

     return MatchingVerifier<Graph,MateMap,VertexIndexMap>::verify_matching(g,mate,vm);    

   }

   template <typename Graph, typename MateMap, typename VertexIndexMap>

   inline bool checked_edmonds_maximum_cardinality_matching(const Graph& g, MateMap mate, VertexIndexMap vm)

   {

     return matching 

       < Graph, MateMap, VertexIndexMap,

         edmonds_augmenting_path_finder, extra_greedy_matching, maximum_cardinality_matching_verifier>

       (g, mate, vm);

   }

   template <typename Graph, typename MateMap>

   inline bool checked_edmonds_maximum_cardinality_matching(const Graph& g, MateMap mate)

   {

     return checked_edmonds_maximum_cardinality_matching(g, mate, get(vertex_index,g));

   }

   template <typename Graph, typename MateMap, typename VertexIndexMap>

   inline void edmonds_maximum_cardinality_matching(const Graph& g, MateMap mate, VertexIndexMap vm)

   {

     matching < Graph, MateMap, VertexIndexMap,

                edmonds_augmenting_path_finder, extra_greedy_matching, no_matching_verifier>

       (g, mate, vm);

   }

   template <typename Graph, typename MateMap>

   inline void edmonds_maximum_cardinality_matching(const Graph& g, MateMap mate)

   {

     edmonds_maximum_cardinality_matching(g, mate, get(vertex_index,g));

   }

 }//namespace boost

 #endif //BOOST_GRAPH_MAXIMUM_CARDINALITY_MATCHING_HPP