minimum_degree_ordering.hpp

support vector clustering for vc++

//-*-c++-*-
//=======================================================================
// Copyright 1997-2001 University of Notre Dame.
// Authors: Lie-Quan Lee, Jeremy Siek
//
// 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 MINIMUM_DEGREE_ORDERING_HPP
#define MINIMUM_DEGREE_ORDERING_HPP

#include <vector>
#include <cassert>
#include <boost/config.hpp>
#include <boost/pending/bucket_sorter.hpp>
#include <boost/detail/numeric_traits.hpp> // for integer_traits
#include <boost/graph/graph_traits.hpp>
#include <boost/property_map.hpp>

namespace boost {

  namespace detail {

    // 
    // Given a set of n integers (where the integer values range from
    // zero to n-1), we want to keep track of a collection of stacks
    // of integers. It so happens that an integer will appear in at
    // most one stack at a time, so the stacks form disjoint sets.
    // Because of these restrictions, we can use one big array to
    // store all the stacks, intertwined with one another.
    // No allocation/deallocation happens in the push()/pop() methods
    // so this is faster than using std::stack's.
    //
    template <class SignedInteger>
    class Stacks {
      typedef SignedInteger value_type;
      typedef typename std::vector<value_type>::size_type size_type;
    public:
      Stacks(size_type n) : data(n) {}
      
      //: stack 
      class stack {
        typedef typename std::vector<value_type>::iterator Iterator;
      public:
        stack(Iterator _data, const value_type& head)
          :  data(_data), current(head) {}

        // did not use default argument here to avoid internal compiler error
        // in g++.
        stack(Iterator _data)
          : data(_data), current(-(std::numeric_limits<value_type>::max)()) {}
        
        void pop() {
          assert(! empty());
          current = data[current];
        }
        void push(value_type v) {
          data[v] = current; 
          current = v;
        }
        bool empty() {
          return current == -(std::numeric_limits<value_type>::max)(); 
        }
        value_type& top() { return current; }
      private:
        Iterator data;
        value_type current;
      };

      // To return a stack object 
      stack make_stack()
        { return stack(data.begin()); }
    protected:
      std::vector<value_type> data;
    };


    // marker class, a generalization of coloring. 
    //
    // This class is to provide a generalization of coloring which has
    // complexity of amortized constant time to set all vertices' color
    // back to be untagged. It implemented by increasing a tag.
    //
    // The colors are:
    //   not tagged 
    //   tagged
    //   multiple_tagged
    //   done
    //
    template <class SignedInteger, class Vertex, class VertexIndexMap>
    class Marker {
      typedef SignedInteger value_type;
      typedef typename std::vector<value_type>::size_type size_type;
      
      static value_type done() 
      { return (std::numeric_limits<value_type>::max)()/2; }
    public:
      Marker(size_type _num, VertexIndexMap index_map) 
        : tag(1 - (std::numeric_limits<value_type>::max)()),
          data(_num, - (std::numeric_limits<value_type>::max)()),
          id(index_map) {}
      
      void mark_done(Vertex node) { data[get(id, node)] = done(); }
      
      bool is_done(Vertex node) { return data[get(id, node)] == done(); }
      
      void mark_tagged(Vertex node) { data[get(id, node)] = tag; }
      
      void mark_multiple_tagged(Vertex node) { data[get(id, node)] = multiple_tag; }
  
      bool is_tagged(Vertex node) const { return data[get(id, node)] >= tag; }

      bool is_not_tagged(Vertex node) const { return data[get(id, node)] < tag; }

      bool is_multiple_tagged(Vertex node) const 
        { return data[get(id, node)] >= multiple_tag; }

      void increment_tag() {
        const size_type num = data.size();
        ++tag;
        if ( tag >= done() ) {
          tag = 1 - (std::numeric_limits<value_type>::max)();
          for (size_type i = 0; i < num; ++i)
            if ( data[i] < done() ) 
              data[i] = - (std::numeric_limits<value_type>::max)();
        }
      }
      
      void set_multiple_tag(value_type mdeg0) 
      { 
        const size_type num = data.size();
        multiple_tag = tag + mdeg0; 
        
        if ( multiple_tag >= done() ) {
          tag = 1-(std::numeric_limits<value_type>::max)();
          
          for (size_type i=0; i<num; i++)
            if ( data[i] < done() ) 
              data[i] = -(std::numeric_limits<value_type>::max)();
          
          multiple_tag = tag + mdeg0; 
        }
      }
      
      void set_tag_as_multiple_tag() { tag = multiple_tag; }
      
    protected:
      value_type tag;
      value_type multiple_tag;
      std::vector<value_type> data;
      VertexIndexMap id;
    };
    
    template< class Iterator, class SignedInteger, 
       class Vertex, class VertexIndexMap, int offset = 1 >
    class Numbering {
      typedef SignedInteger number_type;
      number_type num; //start from 1 instead of zero
      Iterator   data;
      number_type max_num;
      VertexIndexMap id;
    public:
      Numbering(Iterator _data, number_type _max_num, VertexIndexMap id) 
        : num(1), data(_data), max_num(_max_num), id(id) {}
      void operator()(Vertex node) { data[get(id, node)] = -num; }
      bool all_done(number_type i = 0) const { return num + i > max_num; }
      void increment(number_type i = 1) { num += i; }
      bool is_numbered(Vertex node) const {
        return data[get(id, node)] < 0;
      }
      void indistinguishable(Vertex i, Vertex j) {
        data[get(id, i)] = - (get(id, j) + offset);
      }
    };

    template <class SignedInteger, class Vertex, class VertexIndexMap>
    class degreelists_marker {
    public:
      typedef SignedInteger value_type;
      typedef typename std::vector<value_type>::size_type size_type;
      degreelists_marker(size_type n, VertexIndexMap id)
        : marks(n, 0), id(id) {}
      void mark_need_update(Vertex i) { marks[get(id, i)] = 1;  }
      bool need_update(Vertex i) { return marks[get(id, i)] == 1; }
      bool outmatched_or_done (Vertex i) { return marks[get(id, i)] == -1; }
      void mark(Vertex i) { marks[get(id, i)] = -1; }
      void unmark(Vertex i) { marks[get(id, i)] = 0; }
    private:
      std::vector<value_type> marks;
      VertexIndexMap id;
    };

    // Helper function object for edge removal
    template <class Graph, class MarkerP, class NumberD, class Stack,
      class VertexIndexMap>
    class predicateRemoveEdge1 {
      typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
      typedef typename graph_traits<Graph>::edge_descriptor edge_t;
    public:
      predicateRemoveEdge1(Graph& _g, MarkerP& _marker, 
                           NumberD _numbering, Stack& n_e, VertexIndexMap id)
        : g(&_g), marker(&_marker), numbering(_numbering),
          neighbor_elements(&n_e), id(id) {}

      bool operator()(edge_t e) {
        vertex_t dist = target(e, *g);
        if ( marker->is_tagged(dist) )
          return true;
        marker->mark_tagged(dist);
        if (numbering.is_numbered(dist)) {
          neighbor_elements->push(get(id, dist));
          return true;
        }
        return false;
      }
    private:
      Graph*   g;
      MarkerP* marker;
      NumberD  numbering;
      Stack*   neighbor_elements;
      VertexIndexMap id;
    };

    // Helper function object for edge removal
    template <class Graph, class MarkerP>
    class predicate_remove_tagged_edges
    {
      typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
      typedef typename graph_traits<Graph>::edge_descriptor edge_t;
    public:
      predicate_remove_tagged_edges(Graph& _g, MarkerP& _marker)
        : g(&_g), marker(&_marker) {}

      bool operator()(edge_t e) {
        vertex_t dist = target(e, *g);
        if ( marker->is_tagged(dist) )
          return true;
        return false;
      }
    private:
      Graph*   g;
      MarkerP* marker;
    };

    template<class Graph, class DegreeMap, 
             class InversePermutationMap, 
             class PermutationMap,
             class SuperNodeMap, 
             class VertexIndexMap>
    class mmd_impl
    {
      // Typedefs
      typedef graph_traits<Graph> Traits;
      typedef typename Traits::vertices_size_type size_type;
      typedef typename detail::integer_traits<size_type>::difference_type 
        diff_t;
      typedef typename Traits::vertex_descriptor vertex_t;
      typedef typename Traits::adjacency_iterator adj_iter;
      typedef iterator_property_map<vertex_t*, 
        identity_property_map, vertex_t, vertex_t&> IndexVertexMap;
      typedef detail::Stacks<diff_t> Workspace;
      typedef bucket_sorter<size_type, vertex_t, DegreeMap, VertexIndexMap> 
        DegreeLists;
      typedef Numbering<InversePermutationMap, diff_t, vertex_t,VertexIndexMap>
        NumberingD;
      typedef degreelists_marker<diff_t, vertex_t, VertexIndexMap> 
        DegreeListsMarker;
      typedef Marker<diff_t, vertex_t, VertexIndexMap> MarkerP;

      // Data Members

      // input parameters
      Graph& G;
      int delta;
      DegreeMap degree;
      InversePermutationMap inverse_perm;
      PermutationMap perm;
      SuperNodeMap supernode_size;
      VertexIndexMap vertex_index_map;

      // internal data-structures
      std::vector<vertex_t> index_vertex_vec;
      size_type n;
      IndexVertexMap index_vertex_map;
      DegreeLists degreelists;
      NumberingD numbering;
      DegreeListsMarker degree_lists_marker;
      MarkerP marker;
      Workspace work_space;
    public:
      mmd_impl(Graph& g, size_type n_, int delta, DegreeMap degree, 
               InversePermutationMap inverse_perm, 
               PermutationMap perm,
               SuperNodeMap supernode_size, 
               VertexIndexMap id) 
        : G(g), delta(delta), degree(degree), 
        inverse_perm(inverse_perm), 
        perm(perm), 
        supernode_size(supernode_size), 
        vertex_index_map(id),
        index_vertex_vec(n_), 
        n(n_),
        degreelists(n_ + 1, n_, degree, id),
        numbering(inverse_perm, n_, vertex_index_map),
        degree_lists_marker(n_, vertex_index_map), 
        marker(n_, vertex_index_map),
        work_space(n_)
      {
        typename graph_traits<Graph>::vertex_iterator v, vend;
        size_type vid = 0;
        for (tie(v, vend) = vertices(G); v != vend; ++v, ++vid)
          index_vertex_vec[vid] = *v;
        index_vertex_map = IndexVertexMap(&index_vertex_vec[0]);

        // Initialize degreelists.  Degreelists organizes the nodes
        // according to their degree.
        for (tie(v, vend) = vertices(G); v != vend; ++v) {
          put(degree, *v, out_degree(*v, G));
          degreelists.push(*v);
        }
      }

      void do_mmd()
      {
        // Eliminate the isolated nodes -- these are simply the nodes
        // with no neighbors, which are accessible as a list (really, a