max_cardinality_matching.hpp

support vector clustering for vc++

      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