williamr@2: //=======================================================================
williamr@2: // Copyright 1997, 1998, 1999, 2000 University of Notre Dame.
williamr@2: // Copyright 2004, 2005 Trustees of Indiana University
williamr@2: // Authors: Andrew Lumsdaine, Lie-Quan Lee, Jeremy G. Siek,
williamr@2: //          Doug Gregor, D. Kevin McGrath
williamr@2: //
williamr@2: // Distributed under the Boost Software License, Version 1.0. (See
williamr@2: // accompanying file LICENSE_1_0.txt or copy at
williamr@2: // http://www.boost.org/LICENSE_1_0.txt)
williamr@2: //=======================================================================//
williamr@2: #ifndef BOOST_GRAPH_KING_HPP
williamr@2: #define BOOST_GRAPH_KING_HPP
williamr@2: 
williamr@2: #include <boost/config.hpp>
williamr@2: #include <boost/graph/detail/sparse_ordering.hpp>
williamr@2: 
williamr@2: /*
williamr@2:   King Algorithm for matrix reordering
williamr@2: */
williamr@2: 
williamr@2: namespace boost {
williamr@2:   namespace detail {
williamr@2:     template<typename OutputIterator, typename Buffer, typename Compare, 
williamr@2:              typename PseudoDegreeMap, typename VecMap, typename VertexIndexMap>
williamr@2:     class bfs_king_visitor:public default_bfs_visitor
williamr@2:     {
williamr@2:     public:
williamr@2:       bfs_king_visitor(OutputIterator *iter, Buffer *b, Compare compare, 
williamr@2:                        PseudoDegreeMap deg, std::vector<int> loc, VecMap color, 
williamr@2:                        VertexIndexMap vertices): 
williamr@2:         permutation(iter), Qptr(b), degree(deg), comp(compare), 
williamr@2:         Qlocation(loc), colors(color), vertex_map(vertices) { }
williamr@2:       
williamr@2:       template <typename Vertex, typename Graph>
williamr@2:       void finish_vertex(Vertex, Graph& g) {
williamr@2:         typename graph_traits<Graph>::out_edge_iterator ei, ei_end;
williamr@2:         Vertex v, w;
williamr@2: 
williamr@2:         typedef typename std::deque<Vertex>::iterator iterator;
williamr@2:         typedef typename std::deque<Vertex>::reverse_iterator reverse_iterator;
williamr@2: 
williamr@2:         reverse_iterator rend = Qptr->rend()-index_begin;
williamr@2:         reverse_iterator rbegin = Qptr->rbegin();
williamr@2: 
williamr@2: 
williamr@2:         //heap the vertices already there
williamr@2:         std::make_heap(rbegin, rend, boost::bind<bool>(comp, _2, _1));
williamr@2: 
williamr@2:         unsigned i = 0;
williamr@2:         
williamr@2:         for(i = index_begin; i != Qptr->size(); ++i){
williamr@2:           colors[get(vertex_map, (*Qptr)[i])] = 1;
williamr@2:           Qlocation[get(vertex_map, (*Qptr)[i])] = i;
williamr@2:         }
williamr@2: 
williamr@2:         i = 0;
williamr@2: 
williamr@2:         for( ; rbegin != rend; rend--){
williamr@2:           percolate_down<Vertex>(i);
williamr@2:           w = (*Qptr)[index_begin+i];
williamr@2:           for (tie(ei, ei_end) = out_edges(w, g); ei != ei_end; ++ei) {
williamr@2:             v = target(*ei, g);
williamr@2:             put(degree, v, get(degree, v) - 1);
williamr@2:     
williamr@2:             if (colors[get(vertex_map, v)] == 1) {
williamr@2:               percolate_up<Vertex>(get(vertex_map, v), i);            
williamr@2:             }
williamr@2:           }
williamr@2:           
williamr@2:           colors[get(vertex_map, w)] = 0;
williamr@2:           i++;
williamr@2:         }
williamr@2:       }
williamr@2:     
williamr@2:       template <typename Vertex, typename Graph>
williamr@2:       void examine_vertex(Vertex u, const Graph&) {
williamr@2:         
williamr@2:         *(*permutation)++ = u;
williamr@2:         index_begin = Qptr->size();
williamr@2:         
williamr@2:       }
williamr@2:     protected:
williamr@2: 
williamr@2: 
williamr@2:       //this function replaces pop_heap, and tracks state information
williamr@2:       template <typename Vertex>
williamr@2:       void percolate_down(int offset){
williamr@2:         typedef typename std::deque<Vertex>::reverse_iterator reverse_iterator;
williamr@2:         
williamr@2:         int heap_last = index_begin + offset;
williamr@2:         int heap_first = Qptr->size() - 1;
williamr@2:         
williamr@2:         //pop_heap functionality:
williamr@2:         //swap first, last
williamr@2:         std::swap((*Qptr)[heap_last], (*Qptr)[heap_first]);
williamr@2:         
williamr@2:         //swap in the location queue
williamr@2:         std::swap(Qlocation[heap_first], Qlocation[heap_last]);
williamr@2: 
williamr@2:         //set drifter, children
williamr@2:         int drifter = heap_first;
williamr@2:         int drifter_heap = Qptr->size() - drifter;
williamr@2: 
williamr@2:         int right_child_heap = drifter_heap * 2 + 1;
williamr@2:         int right_child = Qptr->size() - right_child_heap;
williamr@2: 
williamr@2:         int left_child_heap = drifter_heap * 2;
williamr@2:         int left_child = Qptr->size() - left_child_heap;
williamr@2: 
williamr@2:         //check that we are staying in the heap
williamr@2:         bool valid = (right_child < heap_last) ? false : true;
williamr@2:         
williamr@2:         //pick smallest child of drifter, and keep in mind there might only be left child
williamr@2:         int smallest_child = (valid && get(degree, (*Qptr)[left_child]) > get(degree,(*Qptr)[right_child])) ? 
williamr@2:           right_child : left_child;
williamr@2:         
williamr@2:         while(valid && smallest_child < heap_last && comp((*Qptr)[drifter], (*Qptr)[smallest_child])){
williamr@2:           
williamr@2:           //if smallest child smaller than drifter, swap them
williamr@2:           std::swap((*Qptr)[smallest_child], (*Qptr)[drifter]);
williamr@2:           std::swap(Qlocation[drifter], Qlocation[smallest_child]);
williamr@2: 
williamr@2:           //update the values, run again, as necessary
williamr@2:           drifter = smallest_child;
williamr@2:           drifter_heap = Qptr->size() - drifter;
williamr@2: 
williamr@2:           right_child_heap = drifter_heap * 2 + 1;
williamr@2:           right_child = Qptr->size() - right_child_heap;
williamr@2: 
williamr@2:           left_child_heap = drifter_heap * 2;
williamr@2:           left_child = Qptr->size() - left_child_heap;
williamr@2: 
williamr@2:           valid = (right_child < heap_last) ? false : true;
williamr@2: 
williamr@2:           smallest_child = (valid && get(degree, (*Qptr)[left_child]) > get(degree,(*Qptr)[right_child])) ? 
williamr@2:             right_child : left_child;
williamr@2:         }
williamr@2: 
williamr@2:       }
williamr@2: 
williamr@2: 
williamr@2:       
williamr@2:       // this is like percolate down, but we always compare against the
williamr@2:       // parent, as there is only a single choice
williamr@2:       template <typename Vertex>
williamr@2:       void percolate_up(int vertex, int offset){
williamr@2:         
williamr@2:         int child_location = Qlocation[vertex];
williamr@2:         int heap_child_location = Qptr->size() - child_location;
williamr@2:         int heap_parent_location = (int)(heap_child_location/2);
williamr@2:         unsigned parent_location = Qptr->size() - heap_parent_location; 
williamr@2: 
williamr@2:         bool valid = (heap_parent_location != 0 && child_location > index_begin + offset && 
williamr@2:                       parent_location < Qptr->size());
williamr@2: 
williamr@2:         while(valid && comp((*Qptr)[child_location], (*Qptr)[parent_location])){
williamr@2:           
williamr@2:           //swap in the heap
williamr@2:           std::swap((*Qptr)[child_location], (*Qptr)[parent_location]);
williamr@2:           
williamr@2:           //swap in the location queue
williamr@2:           std::swap(Qlocation[child_location], Qlocation[parent_location]);
williamr@2: 
williamr@2:           child_location = parent_location;
williamr@2:           heap_child_location = heap_parent_location;
williamr@2:           heap_parent_location = (int)(heap_child_location/2);
williamr@2:           parent_location = Qptr->size() - heap_parent_location; 
williamr@2:           valid = (heap_parent_location != 0 && child_location > index_begin + offset);
williamr@2:         }
williamr@2:       }
williamr@2:       
williamr@2:       OutputIterator *permutation;
williamr@2:       int index_begin;
williamr@2:       Buffer *Qptr;
williamr@2:       PseudoDegreeMap degree;
williamr@2:       Compare comp;
williamr@2:       std::vector<int> Qlocation;
williamr@2:       VecMap colors;
williamr@2:       VertexIndexMap vertex_map;
williamr@2:     };
williamr@2:   
williamr@2: 
williamr@2:   } // namespace detail  
williamr@2:   
williamr@2: 
williamr@2:   template<class Graph, class OutputIterator, class ColorMap, class DegreeMap,
williamr@2:            typename VertexIndexMap> 
williamr@2:   OutputIterator
williamr@2:   king_ordering(const Graph& g,
williamr@2:                 std::deque< typename graph_traits<Graph>::vertex_descriptor >
williamr@2:                   vertex_queue,
williamr@2:                 OutputIterator permutation, 
williamr@2:                 ColorMap color, DegreeMap degree,
williamr@2:                 VertexIndexMap index_map)
williamr@2:   {
williamr@2:     typedef typename property_traits<DegreeMap>::value_type ds_type;
williamr@2:     typedef typename property_traits<ColorMap>::value_type ColorValue;
williamr@2:     typedef color_traits<ColorValue> Color;
williamr@2:     typedef typename graph_traits<Graph>::vertex_descriptor Vertex;
williamr@2:     typedef iterator_property_map<typename std::vector<ds_type>::iterator, VertexIndexMap, ds_type, ds_type&> PseudoDegreeMap;
williamr@2:     typedef indirect_cmp<PseudoDegreeMap, std::less<ds_type> > Compare;
williamr@2:     typedef typename boost::sparse::sparse_ordering_queue<Vertex> queue;
williamr@2:     typedef typename detail::bfs_king_visitor<OutputIterator, queue, Compare,             
williamr@2:       PseudoDegreeMap, std::vector<int>, VertexIndexMap > Visitor;
williamr@2:     typedef typename graph_traits<Graph>::vertices_size_type
williamr@2:       vertices_size_type;
williamr@2:     std::vector<ds_type> pseudo_degree_vec(num_vertices(g));
williamr@2:     PseudoDegreeMap pseudo_degree(pseudo_degree_vec.begin(), index_map);
williamr@2:     
williamr@2:     typename graph_traits<Graph>::vertex_iterator ui, ui_end;    
williamr@2:     queue Q;
williamr@2:     // Copy degree to pseudo_degree
williamr@2:     // initialize the color map
williamr@2:     for (tie(ui, ui_end) = vertices(g); ui != ui_end; ++ui){
williamr@2:       put(pseudo_degree, *ui, get(degree, *ui));
williamr@2:       put(color, *ui, Color::white());
williamr@2:     }
williamr@2:     
williamr@2:     Compare comp(pseudo_degree);
williamr@2:     std::vector<int> colors(num_vertices(g));
williamr@2: 
williamr@2:     for(vertices_size_type i = 0; i < num_vertices(g); i++) 
williamr@2:       colors[i] = 0;
williamr@2: 
williamr@2:     std::vector<int> loc(num_vertices(g));
williamr@2: 
williamr@2:     //create the visitor
williamr@2:     Visitor vis(&permutation, &Q, comp, pseudo_degree, loc, colors, index_map);
williamr@2:     
williamr@2:     while( !vertex_queue.empty() ) {
williamr@2:       Vertex s = vertex_queue.front();
williamr@2:       vertex_queue.pop_front();
williamr@2:       
williamr@2:       //call BFS with visitor
williamr@2:       breadth_first_visit(g, s, Q, vis, color);
williamr@2:     }
williamr@2: 
williamr@2:     return permutation;
williamr@2:   }
williamr@2: 
williamr@2:   
williamr@2:   // This is the case where only a single starting vertex is supplied.
williamr@2:   template <class Graph, class OutputIterator,
williamr@2:             class ColorMap, class DegreeMap, typename VertexIndexMap>
williamr@2:   OutputIterator
williamr@2:   king_ordering(const Graph& g,
williamr@2:                 typename graph_traits<Graph>::vertex_descriptor s,
williamr@2:                 OutputIterator permutation, 
williamr@2:                 ColorMap color, DegreeMap degree, VertexIndexMap index_map)
williamr@2:   {
williamr@2: 
williamr@2:     std::deque< typename graph_traits<Graph>::vertex_descriptor > vertex_queue;
williamr@2:     vertex_queue.push_front( s );
williamr@2:     return king_ordering(g, vertex_queue, permutation, color, degree,
williamr@2:                          index_map);
williamr@2:   }
williamr@2: 
williamr@2:   
williamr@2:   template < class Graph, class OutputIterator, 
williamr@2:              class ColorMap, class DegreeMap, class VertexIndexMap>
williamr@2:   OutputIterator 
williamr@2:   king_ordering(const Graph& G, OutputIterator permutation, 
williamr@2:                 ColorMap color, DegreeMap degree, VertexIndexMap index_map)
williamr@2:   {
williamr@2:     if (vertices(G).first == vertices(G).second)
williamr@2:       return permutation;
williamr@2: 
williamr@2:     typedef typename boost::graph_traits<Graph>::vertex_descriptor Vertex;
williamr@2:     typedef typename boost::graph_traits<Graph>::vertex_iterator   VerIter;
williamr@2:     typedef typename property_traits<ColorMap>::value_type ColorValue;
williamr@2:     typedef color_traits<ColorValue> Color;
williamr@2: 
williamr@2:     std::deque<Vertex>      vertex_queue;
williamr@2: 
williamr@2:     // Mark everything white
williamr@2:     BGL_FORALL_VERTICES_T(v, G, Graph) put(color, v, Color::white());
williamr@2: 
williamr@2:     // Find one vertex from each connected component 
williamr@2:     BGL_FORALL_VERTICES_T(v, G, Graph) {
williamr@2:       if (get(color, v) == Color::white()) {
williamr@2:         depth_first_visit(G, v, dfs_visitor<>(), color);
williamr@2:         vertex_queue.push_back(v);
williamr@2:       }
williamr@2:     }
williamr@2: 
williamr@2:     // Find starting nodes for all vertices
williamr@2:     // TBD: How to do this with a directed graph?
williamr@2:     for (typename std::deque<Vertex>::iterator i = vertex_queue.begin();
williamr@2:          i != vertex_queue.end(); ++i)
williamr@2:       *i = find_starting_node(G, *i, color, degree);
williamr@2:     
williamr@2:     return king_ordering(G, vertex_queue, permutation, color, degree,
williamr@2:                          index_map);
williamr@2:   }
williamr@2: 
williamr@2:   template<typename Graph, typename OutputIterator, typename VertexIndexMap>
williamr@2:   OutputIterator 
williamr@2:   king_ordering(const Graph& G, OutputIterator permutation, 
williamr@2:                 VertexIndexMap index_map)
williamr@2:   {
williamr@2:     if (vertices(G).first == vertices(G).second)
williamr@2:       return permutation;
williamr@2: 
williamr@2:     typedef out_degree_property_map<Graph> DegreeMap;
williamr@2:     std::vector<default_color_type> colors(num_vertices(G));
williamr@2:     return king_ordering(G, permutation, 
williamr@2:                          make_iterator_property_map(&colors[0], index_map,
williamr@2:                                                     colors[0]),
williamr@2:                          make_out_degree_map(G), index_map);
williamr@2:   }
williamr@2: 
williamr@2:   template<typename Graph, typename OutputIterator>
williamr@2:   inline OutputIterator 
williamr@2:   king_ordering(const Graph& G, OutputIterator permutation)
williamr@2:   { return king_ordering(G, permutation, get(vertex_index, G)); }
williamr@2: 
williamr@2: } // namespace boost
williamr@2: 
williamr@2: 
williamr@2: #endif // BOOST_GRAPH_KING_HPP