max_cardinality_matching.hpp

support vector clustering for vc++

//=======================================================================
// 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]);