//-*-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

         // stack) at location 0.  Since these don't affect any other

         // nodes, we can eliminate them without doing degree updates.

         typename DegreeLists::stack list_isolated = degreelists[0];

         while (!list_isolated.empty()) {

           vertex_t node = list_isolated.top();

           marker.mark_done(node);

           numbering(node);

           numbering.increment();

           list_isolated.pop();

         }

         size_type min_degree = 1;

         typename DegreeLists::stack list_min_degree = degreelists[min_degree];

         while (list_min_degree.empty()) {

           ++min_degree;

           list_min_degree = degreelists[min_degree];

         }

         // check if the whole eliminating process is done

         while (!numbering.all_done()) {

           size_type min_degree_limit = min_degree + delta; // WARNING

           typename Workspace::stack llist = work_space.make_stack();

           // multiple elimination

           while (delta >= 0) {

             // Find the next non-empty degree

             for (list_min_degree = degreelists[min_degree]; 

                  list_min_degree.empty() && min_degree <= min_degree_limit; 

                  ++min_degree, list_min_degree = degreelists[min_degree])

               ;

             if (min_degree > min_degree_limit)

               break;

             const vertex_t node = list_min_degree.top();

             const size_type node_id = get(vertex_index_map, node);

             list_min_degree.pop();

             numbering(node);

             // check if node is the last one

             if (numbering.all_done(supernode_size[node])) {

               numbering.increment(supernode_size[node]);

               break;

             }

             marker.increment_tag();

             marker.mark_tagged(node);

             this->eliminate(node);

             numbering.increment(supernode_size[node]);

             llist.push(node_id);

           } // multiple elimination

           if (numbering.all_done()) 

             break;

           this->update( llist, min_degree);

         }

       } // do_mmd()

       void eliminate(vertex_t node)

       {

         typename Workspace::stack element_neighbor = work_space.make_stack();

         // Create two function objects for edge removal

         typedef typename Workspace::stack WorkStack;

         predicateRemoveEdge1<Graph, MarkerP, NumberingD, 

                              WorkStack, VertexIndexMap>

           p(G, marker, numbering, element_neighbor, vertex_index_map);

         predicate_remove_tagged_edges<Graph, MarkerP> p2(G, marker);

         // Reconstruct the adjacent node list, push element neighbor in a List.

         remove_out_edge_if(node, p, G);

         //during removal element neighbors are collected.

         while (!element_neighbor.empty()) {

           // element absorb

           size_type e_id = element_neighbor.top();

           vertex_t element = get(index_vertex_map, e_id);

           adj_iter i, i_end;

           for (tie(i, i_end) = adjacent_vertices(element, G); i != i_end; ++i){

             vertex_t i_node = *i;

             if (!marker.is_tagged(i_node) && !numbering.is_numbered(i_node)) {

               marker.mark_tagged(i_node);

               add_edge(node, i_node, G);

             }

           }

           element_neighbor.pop();

         }

         adj_iter v, ve;

         for (tie(v, ve) = adjacent_vertices(node, G); v != ve; ++v) {

           vertex_t v_node = *v;

           if (!degree_lists_marker.need_update(v_node) 

               && !degree_lists_marker.outmatched_or_done(v_node)) {

             degreelists.remove(v_node);

           }

           //update out edges of v_node

           remove_out_edge_if(v_node, p2, G);

           if ( out_degree(v_node, G) == 0 ) { // indistinguishable nodes

             supernode_size[node] += supernode_size[v_node];

             supernode_size[v_node] = 0;

             numbering.indistinguishable(v_node, node);

             marker.mark_done(v_node);

             degree_lists_marker.mark(v_node);

           } else {                           // not indistinguishable nodes

             add_edge(v_node, node, G);

             degree_lists_marker.mark_need_update(v_node);

           }

         }

       } // eliminate()

       template <class Stack>

       void update(Stack llist, size_type& min_degree)

       {

         size_type min_degree0 = min_degree + delta + 1;

         while (! llist.empty()) {

           size_type deg, deg0 = 0;

           marker.set_multiple_tag(min_degree0);

           typename Workspace::stack q2list = work_space.make_stack();

           typename Workspace::stack qxlist = work_space.make_stack();

           vertex_t current = get(index_vertex_map, llist.top());

           adj_iter i, ie;

           for (tie(i,ie) = adjacent_vertices(current, G); i != ie; ++i) {

             vertex_t i_node = *i;

             const size_type i_id = get(vertex_index_map, i_node);

             if (supernode_size[i_node] != 0) {

               deg0 += supernode_size[i_node];

               marker.mark_multiple_tagged(i_node);

               if (degree_lists_marker.need_update(i_node)) {

                 if (out_degree(i_node, G) == 2)

                   q2list.push(i_id);

                 else

                   qxlist.push(i_id);

               }

             }

           }

           while (!q2list.empty()) {

             const size_type u_id = q2list.top();

             vertex_t u_node = get(index_vertex_map, u_id);

             // if u_id is outmatched by others, no need to update degree

             if (degree_lists_marker.outmatched_or_done(u_node)) {

               q2list.pop();

               continue;

             }

             marker.increment_tag();

             deg = deg0;

             adj_iter nu = adjacent_vertices(u_node, G).first;

             vertex_t neighbor = *nu;

             if (neighbor == u_node) {

               ++nu;

               neighbor = *nu;

             }

             if (numbering.is_numbered(neighbor)) {

               adj_iter i, ie;

               for (tie(i,ie) = adjacent_vertices(neighbor, G);

                    i != ie; ++i) {

                 const vertex_t i_node = *i;

                 if (i_node == u_node || supernode_size[i_node] == 0)

                   continue;

                 if (marker.is_tagged(i_node)) {

                   if (degree_lists_marker.need_update(i_node)) {

                     if ( out_degree(i_node, G) == 2 ) { // is indistinguishable

                       supernode_size[u_node] += supernode_size[i_node];

                       supernode_size[i_node] = 0;

                       numbering.indistinguishable(i_node, u_node);

                       marker.mark_done(i_node);

                       degree_lists_marker.mark(i_node);

                     } else                        // is outmatched

                       degree_lists_marker.mark(i_node);

                   }

                 } else {

                   marker.mark_tagged(i_node);

                   deg += supernode_size[i_node];

                 }

               }

             } else

               deg += supernode_size[neighbor];

             deg -= supernode_size[u_node];

             degree[u_node] = deg; //update degree

             degreelists[deg].push(u_node);

             //u_id has been pushed back into degreelists

             degree_lists_marker.unmark(u_node);

             if (min_degree > deg) 

               min_degree = deg;

             q2list.pop();

           } // while (!q2list.empty())

           while (!qxlist.empty()) {

             const size_type u_id = qxlist.top();

             const vertex_t u_node = get(index_vertex_map, u_id);

             // if u_id is outmatched by others, no need to update degree

             if (degree_lists_marker.outmatched_or_done(u_node)) {

               qxlist.pop();

               continue;

             }

             marker.increment_tag();

             deg = deg0;

             adj_iter i, ie;

             for (tie(i, ie) = adjacent_vertices(u_node, G); i != ie; ++i) {

               vertex_t i_node = *i;

               if (marker.is_tagged(i_node)) 

                 continue;

               marker.mark_tagged(i_node);

               if (numbering.is_numbered(i_node)) {

                 adj_iter j, je;

                 for (tie(j, je) = adjacent_vertices(i_node, G); j != je; ++j) {

                   const vertex_t j_node = *j;

                   if (marker.is_not_tagged(j_node)) {

                     marker.mark_tagged(j_node);

                     deg += supernode_size[j_node];

                   }

                 }

               } else

                 deg += supernode_size[i_node];

             } // for adjacent vertices of u_node

             deg -= supernode_size[u_node];

             degree[u_node] = deg;

             degreelists[deg].push(u_node);

             // u_id has been pushed back into degreelists

             degree_lists_marker.unmark(u_node); 

             if (min_degree > deg)

               min_degree = deg;

             qxlist.pop();

           } // while (!qxlist.empty()) {

           marker.set_tag_as_multiple_tag();

           llist.pop();

         } // while (! llist.empty())

       } // update()

       void build_permutation(InversePermutationMap next,

                              PermutationMap prev) 

       {

         // collect the permutation info

         size_type i;

         for (i = 0; i < n; ++i) {

           diff_t size = supernode_size[get(index_vertex_map, i)];

           if ( size <= 0 ) {

             prev[i] = next[i];

             supernode_size[get(index_vertex_map, i)]

               = next[i] + 1;  // record the supernode info

           } else

             prev[i] = - next[i];

         }

         for (i = 1; i < n + 1; ++i) {

           if ( prev[i-1] > 0 )

             continue;

           diff_t parent = i;

           while ( prev[parent - 1] < 0 ) {

             parent = - prev[parent - 1];

           }

           diff_t root = parent;

           diff_t num = prev[root - 1] + 1;

           next[i-1] = - num;

           prev[root-1] = num;

           parent = i;

           diff_t next_node = - prev[parent - 1];

           while (next_node > 0) {

             prev[parent-1] = - root;

             parent = next_node;

             next_node = - prev[parent - 1];

           }

         }

         for (i = 0; i < n; i++) {

           diff_t num = - next[i] - 1;

           next[i] = num;

           prev[num] = i;

         }

       } // build_permutation()

     };

   } //namespace detail

   // MMD algorithm

   // 

   //The implementation presently includes the enhancements for mass

   //elimination, incomplete degree update, multiple elimination, and

   //external degree.

   //

   //Important Note: This implementation requires the BGL graph to be

   //directed.  Therefore, nonzero entry (i, j) in a symmetrical matrix

   //A coresponds to two directed edges (i->j and j->i).

   //

   //see Alan George and Joseph W. H. Liu, The Evolution of the Minimum

   //Degree Ordering Algorithm, SIAM Review, 31, 1989, Page 1-19

   template<class Graph, class DegreeMap, 

            class InversePermutationMap, 

            class PermutationMap,

            class SuperNodeMap, class VertexIndexMap>

   void minimum_degree_ordering

     (Graph& G, 

      DegreeMap degree, 

      InversePermutationMap inverse_perm, 

      PermutationMap perm, 

      SuperNodeMap supernode_size, 

      int delta, 

      VertexIndexMap vertex_index_map)

   {

     detail::mmd_impl<Graph,DegreeMap,InversePermutationMap,

       PermutationMap, SuperNodeMap, VertexIndexMap> 

       impl(G, num_vertices(G), delta, degree, inverse_perm, 

            perm, supernode_size, vertex_index_map);

     impl.do_mmd();

     impl.build_permutation(inverse_perm, perm);

   }

 } // namespace boost

 #endif // MINIMUM_DEGREE_ORDERING_HPP