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