//=======================================================================
// 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/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(boost::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)
  {
    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( boost::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&, MateMap, VertexIndexMap)
    { }

    inline bool augment_matching() { return false; }

    template <typename PropertyMap>
    void get_current_matching(PropertyMap) {}
  };


  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(boost::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(boost::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(boost::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( boost::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(boost::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(boost::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(boost::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(boost::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( boost::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(boost::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(boost::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(boost::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, Graph&) 
      {
        m_parity = false; 
      }
      
      template <class Vertex, class Graph>
      void discover_vertex(Vertex, Graph&) 
      {
        m_parity = !m_parity;
        m_parity ? ++m_count : --m_count;
      }
      
    protected:
      SizeType& m_count;
      
    private:
      bool m_parity;
      
    };

  }//namespace detail


  template <typename Graph, typename MateMap, 
            typename VertexIndexMap = dummy_property_map>
  struct no_matching_verifier
  {
    inline static bool 
    verify_matching(const Graph&, MateMap, VertexIndexMap) 
    { 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(boost::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