//
// $COPYRIGHT$
//
//-*-c++-*-

#ifndef MMD_H
#define MMD_H

#include <vector>
#include <ggcl/mmd_aux.h>
#include <ggcl/ordered_stacks.h>

namespace matrix_ordering {


  /* Helper class for edge removal */
  template <class MarkerT, class NumberDecorator, 
	    class Stack, class ID>
  struct predicateRemoveEdge1 {
    predicateRemoveEdge1(MarkerT _marker, 
			 NumberDecorator _numbering,
			 Stack& n_e, ID _id)
      : marker(_marker), numbering(_numbering),
	neighbor_elements(&n_e), id(_id) {}
  
    template<class Edge>
    bool 
    operator()(Edge e) {
      const typename edge_traits<Edge>::vertex_type dist = target(e);
      const int d_id = id[dist];
      if (marker.is_tagged(d_id)) 
	return true;
    
      marker.mark_tag(d_id);
    
      if (numbering.is_numbered(d_id)) {
	neighbor_elements->push(d_id);
	return true;
      }
    
      return false;
    }
  
    MarkerT marker;
    NumberDecorator numbering;
    Stack* neighbor_elements;
    ID id;
  };


  /* Helper class for edge removal */
  template <class MarkerT, class ID>
  struct predicate_remove_tagged_edges {
    predicate_remove_tagged_edges(MarkerT _marker, ID _id)
      : marker(_marker), id(_id) {}
    template <class Edge>
  
    bool
    operator()(Edge e) {
      const typename edge_traits<Edge>::vertex_type dist = target(e);
      const int d_id = id[dist];
      if (marker.is_tagged(d_id))
	return true;
      return false;
    }
  
    MarkerT marker;
    ID id;
  };


  template <class DegreeLists, class NumberingT,
	    class MarkerT, class SuperNodeSize, 
	    class StackPool>
  class MMD {
  public:
    typedef typename DegreeLists::size_type size_type;

    MMD(DegreeLists d, NumberingT n, MarkerT m, 
	SuperNodeSize q, StackPool s, int _delta )
      : degreelists(d), numbering(n), marker(m), 
      qsize(q), stackpool(s), delta(_delta) {}

    template <class Graph, class ID>
    void run(Graph& G, ID id)  {
      typedef typename graph_traits<Graph>::vertex_type Vertex;

      degreelists.init(G.num_vertices());
      marker.init(G.num_vertices());

      /* 
       * Initialize degreelists.  Degreelists organizes the nodes according to
       * their degree.
       */
      typename graph_traits<Graph>::vertices_type::iterator iter;
      for (iter = vertices(G).begin(); iter != vertices(G).end(); ++iter) {
	const Vertex v = *iter;
	const size_type nid = id[v];
	const size_type met = adj(v).size(); //!!! Metric
	degreelists[met].push(nid);
      }
      
      /* 
       * 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()) {
	const size_type node = list_isolated.top();
	marker.mark_done(node);
	numbering(node);
	numbering.increment();
	list_isolated.pop();
      }
      
      size_type mdeg = 1;
      typename DegreeLists::stack list_mdeg = degreelists[mdeg];
      
      while (list_mdeg.empty()) {
	++mdeg;
	list_mdeg = degreelists[mdeg];
      }
      
      /* check if the whole eliminating process is done */
      while (!numbering.all_done()) {
    
	size_type mdlimit;
	if ( delta >= 0 )
	  mdlimit = mdeg + delta;
	else 
	  mdlimit = mdeg;

	typename StackPool::stack llist = stackpool.make_stack();
	
	/* multiple elimination */
	do {
	  
	  /* Find the next non-empty degree */
	  for (list_mdeg = degreelists[mdeg]; 
	       list_mdeg.empty() && mdeg <= mdlimit; 
	       ++mdeg, list_mdeg = degreelists[mdeg])
	    ;
	  if (mdeg > mdlimit)
	    break;
	  
	  const size_type node = list_mdeg.top();
	  list_mdeg.pop();
	  numbering(node);
	  
	  /* check if node is the last one */
	  if (numbering.all_done(qsize[node])) {
	    numbering.increment(qsize[node]);
	    break;
	  }
	  
	  marker.increment_tag();
	  marker.mark_tag(node);
	  
	  eliminate(G, G.vertex_of_id(node), id);
	  
	  numbering.increment(qsize[node]);
	  llist.push(node);
	} while (delta >= 0); /* multiple elimination */
	
	if (numbering.all_done()) 
	  break;
	
	update(G, llist, mdeg, id);
      }
    }
    
    template <class Graph, class Stack, class ID>
    void update(Graph& G, Stack llist, int& mdeg, ID id) {
      typedef typename graph_traits<Graph>::vertex_type Vertex;
      typedef typename vertex_traits<Vertex>::vertexlist_type AdjVec;
  
      int mdeg0;
      if (delta>=0) 
	mdeg0 = mdeg + delta + 1;
      else 
	mdeg0 = mdeg + 1;

      while (! llist.empty()) {	/*!!!*/
	int deg, deg0 = 0;

	marker.set_mtag(mdeg0);
	typename StackPool::stack q2list = stackpool.make_stack();
	typename StackPool::stack qxlist = stackpool.make_stack();

	Vertex current = G.vertex_of_id(llist.top());
	AdjVec work = adj(current);

	const typename AdjVec::iterator ie = work.end();
	for (typename AdjVec::iterator i = work.begin(); i != ie; ++i) {
	  Vertex inode = *i;
      
	  if (qsize[inode] != 0) {
	    deg0 += qsize[inode];
	
	    const int idi = id[inode];
	
	    marker.mark_mtag(idi);
	    if (degreelists.need_update(idi)) {
	      if (adj(inode).size() == 2)
		q2list.push(idi);
	      else
		qxlist.push(idi);
	    }
	  }
	}
    
	while (!q2list.empty()) { /*!!!*/
	  const int uid = q2list.top();
	  Vertex unode = G.vertex_of_id(uid);
      
	  /* if uid is outmatched by others, no need to update degree */
	  if (degreelists.outmatched_or_done(uid)) {
	    q2list.pop();
	    continue;
	  }
      
	  marker.increment_tag();
	  deg = deg0;
      
	  typename AdjVec::iterator nu = adj(unode).begin();
	  Vertex neighbor = *nu;
	  if (neighbor == unode) {
	    ++nu;
	    neighbor = *nu;
	  }
      
	  if (numbering.is_numbered(id[neighbor])) {
	
	    typename AdjVec::iterator i;
	    const typename AdjVec::iterator ie = adj(neighbor).end();
	    for (i = adj(neighbor).begin(); i != ie; ++i) {
	      const Vertex inode = *i;
	  
	      if (inode == unode || qsize[inode] == 0)
		continue;
	  
	      const int idi = id[inode];
	  
	      if (marker.is_tagged(idi)) {
		if (degreelists.need_update(idi)) {
		  if (adj(inode).size() == 2) { /* is indistinguishable */
		    qsize[unode] += qsize[inode];
		    qsize[inode] = 0;
		    numbering.indistinguishable(idi, id[unode]);
		    marker.mark_done(idi);
		    degreelists.mark(idi);
		  } else		              /*is outmatched */
		    degreelists.mark(idi);
		}
	      } else {
		marker.mark_tag(idi);
		deg += qsize[inode];
	      }
	    }
	  } else
	    deg += qsize[neighbor];
      
	  deg -= qsize[unode];
	  degreelists[deg].push(uid);
	  if (mdeg > deg) 
	    mdeg = deg;
	  q2list.pop();
	}
    
	while (!qxlist.empty()) {	/*!!!*/
	  const int uid = qxlist.top();
	  const Vertex unode = G.vertex_of_id(uid);
      
	  /* if uid is outmatched by others, no need to update degree */
	  if (degreelists.outmatched_or_done(uid)) {
	    qxlist.pop();
	    continue;
	  }
      
	  marker.increment_tag();
	  deg = deg0;

	  typename AdjVec::iterator i;
	  const typename AdjVec::iterator ie = adj(unode).end(); 
	  for (i = adj(unode).begin(); i != ie; ++i) { /*!!!*/
	
	    Vertex inode = *i;
	    const int idi = id[inode];
	
	    if (marker.is_tagged(idi)) 
	      continue;
	    marker.mark_tag(idi);
	
	    if (numbering.is_numbered(idi)) {
	  
	      typename AdjVec::iterator j;
	      const typename AdjVec::iterator je = adj(inode).end();
	      for (j = adj(inode).begin(); j != je; ++j) { /*!!!*/
		const Vertex jnode = *j;
		const int idj = id[jnode];
	    
		if (marker.is_not_tagged(idj)) {
		  marker.mark_tag(idj);
		  deg += qsize[jnode];
		}
	      }
	    } else
	      deg += qsize[inode];
	  }
      
	  deg -= qsize[unode];
	  degreelists[deg].push(uid);
      
	  if (mdeg > deg)
	    mdeg = deg;
      
	  qxlist.pop();
	}
    
	marker.set_tag_as_mtag();
	llist.pop();
      }

    }
    
    template <class Graph, class ID>
    void eliminate(Graph& G, typename graph_traits<Graph>::vertex_type node, 
		   ID id) {
      typedef typename graph_traits<Graph>::vertex_type Vertex;
      typedef typename vertex_traits<Vertex>::edgelist_type container;
      typedef typename vertex_traits<Vertex>::vertexlist_type::iterator 
	      AdjIter;
      container work = out_edges(node);
      
      typename StackPool::stack element_neighbor = stackpool.make_stack();
      
      /* Create two function objects for edge removal */
      predicateRemoveEdge1<MarkerT, NumberingT, 
	typename StackPool::stack, ID>
	p(marker, numbering, element_neighbor, id);
      predicate_remove_tagged_edges<MarkerT, ID> 
	p2(marker, id);
      
      /* 
       * Reconstruct the adjacent node list, push element neighbor in a List;
       * require asignment operator for edge mayfly type
       */
      typename container::iterator pos = 
	std::remove_copy_if(work.begin(), work.end(), work.begin(), p);
      work.erase(pos, work.end());
      
      while (!element_neighbor.empty()) {
	
	/* element absorb */
	int e_id = element_neighbor.top();
	Vertex element = G.vertex_of_id(e_id);
	
	const AdjIter i_end = adj(element).end();
	for (AdjIter i = adj(element).begin(); i != i_end; ++i) {
	  Vertex inode = *i;
	  const int idi = id[inode];
	  
	  if (!marker.is_tagged(idi) && !numbering.is_numbered(idi)) {
	    marker.mark_tag(idi);
	    work.push_back(idi);
	  }
	}
	element_neighbor.pop();
      }
      
      const AdjIter ve =  adj(node).end();
      for (AdjIter v = adj(node).begin(); v != ve; ++v) { /*!!!*/
	Vertex vnode = *v;
	const int idv = id[vnode];
	
	if (!degreelists.need_update(idv) 
	    && !degreelists.outmatched_or_done(idv))
	  degreelists.remove(idv);
	
	container vnode_out_edges = out_edges(vnode);
	typename container::iterator pos = 
	  std::remove_copy_if(vnode_out_edges.begin(), vnode_out_edges.end(), 
			      vnode_out_edges.begin(), p2);
	vnode_out_edges.erase(pos, vnode_out_edges.end());
	
	if (vnode_out_edges.empty()) { /* indistinguishable nodes */
	  qsize[node] += qsize[vnode];
	  qsize[vnode] = 0;
	  numbering.indistinguishable(idv, id[node]);
	  marker.mark_done(idv);
	  degreelists.mark(idv);
	} else {			     /* not indistinguishable nodes */
	  vnode_out_edges.push_back(id[node]);
	  degreelists.mark_need_update(idv);
	}
      }
      
    }
    
  protected:
    DegreeLists              degreelists;
    NumberingT               numbering;
    MarkerT                  marker;
    SuperNodeSize            qsize;
    StackPool                stackpool;
    int                      delta;
  };


  //: MMD algorithm
  // 
  //The implementation presently includes the enhancements for mass
  //elimination, incomplete degree update, multiple elimination, and
  //external degree.
  //
  //!category: ordering
  //!component: function
  //!definition: mmd.h
  //!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 IterI, class IterP, class IterQ>
  void 
  multiple_minimum_degree_ordering(Graph& G, IterI inverse_perm, 
				   IterP perm, IterQ qsize, int delta=0)
  {
    multiple_minimum_degree_ordering(G, inverse_perm, perm, 
				     qsize, delta, id_decorator());
  }

  //: MMD algorithm (with explicit ID decorator)
  // 
  //The implementation presently includes the enhancements for mass
  //elimination, incomplete degree update, multiple elimination, and
  //external degree.
  //
  //!category: ordering
  //!component: function
  //!definition: mmd.h
  //!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 IterI, class IterP, class IterQ, class ID>
  void 
  multiple_minimum_degree_ordering(Graph& G, IterI inverse_perm,
				   IterP perm, IterQ qsize, int delta, ID id)
  {
    typedef std::vector<int>                   Vector;
    typedef Vector::iterator                   Iter;
    typedef ordered_stacks<Iter, IterI, IterP> DegreeLists;

    int              n = G.num_vertices();

    Vector           head(n+1);
    Vector           workspace(n);
    Vector           mark(n);

    IterI            next = inverse_perm;
    IterP            prev = perm;

    DegreeLists      degreelists(head.begin(), next, prev);
    Marker<Iter>     marker(mark.begin(), n);
    Numbering<Iter>  numbering(next, n);
    Stacks<Iter>     stackpool(workspace.begin());

    MMD<DegreeLists, Numbering<Iter>, Marker<Iter>, IterQ, Stacks<Iter> >
      my_mmd(degreelists, numbering, marker, qsize, stackpool, delta);
    
    my_mmd.run(G, id);

    /* qsize[i] == 0,  -next[i]-1 is the parent of i
     *          >  0,  -next[i] is the inverse labeling
     */

    /* collect the permutation info */

    for (int i = 0; i < n; i++) {
      int size   = qsize[i];

      if ( size <= 0 ) {
	prev[i]  = next[i];          /* copy it for internal use */
	qsize[i] = next[i]+1;        /* record the possible representative */
      } else prev[i] = -next[i];     /* copy it for internal use */
    }
  
    for (int i = 1; i < n+1; i++) {
      if ( prev[i-1] > 0 )
	continue;                    /* done for this node */
      
      int parent   = i;
      while ( prev[parent-1] < 0 )
	parent     = -prev[parent-1];
      /*
        Above is to find the real supernode representative..
        If parent points to other node, parent is not the real 
	representative.
       */

      int root     = parent;
      int num      = prev[root-1] + 1;
      next[i-1]    = -num;
      prev[root-1] = num;
      /*
	Above is to number that node.  

	It does not change root labeling since its inverse labeling is
	in next[root]. Here root is the representative of a supernode.

	At the beginning, prev[root] is equal to inverse labeling of
	root. When i is found to be in supernode root, it assigns the
	inverse labeling of i to be prev[root] plus 1 and update
	prev[root]. Therefore, all non-representative nodes are
	numbered after the representative node.
      */

      parent = i;
      int next_node = - prev[parent-1];
      while (next_node > 0) {
	prev[parent-1]  = -root;     /* path compression in graph landuage */
	qnode[parent-1] = -root+1;   /* record the real supernode */
	parent          = next_node;
	next_node       = - prev[parent-1];
      }
    }
    /*
      Above is the path compression to speed up the process.
      Once a node found its representative, record it so that 
      next time it needs not to traverse all other nonrepresentative 
      nodes.
     */

    
    for (int i = 0; i < n; i++) {
      int num         = - next[i]-1;
      inverse_perm[i] = num;
      perm[num]       = i;
    }

    /* Conclusion:
     
       inverse_perm[i] is the inverse labeling,  i is in [0, n).
       perm[i] is the labeling, i is in [0, n).
       if qsize[i] > 0, i is the supernode and qsize[i] is the
                        number of nodes in this supernode.
       if qsize[i] < 0, i is in the supernode -qsize[i]. 
    */
  }
  

} /* namespace matrix_ordering */

#endif /*MMD_H*/