GenericTreeNode.h

/** \file GenericTreeNode.h
    This file defines the generic tree structures.
    @author Filippo Gioachin (previously Graeme Lufkin and Chao Huang)
    @version 2.0
*/

#ifndef GENERICTREENODE_H
#define GENERICTREENODE_H

#include "pup.h"

#include <map>
#include <vector>
#include <algorithm>
#include <set>
#include <list>

#include "OrientedBox.h"
#include "MultipoleMoments.h"
#include "keytype.h"
#include "GravityParticle.h"

namespace Tree {

  /** @brief This key is the identification of a node inside the global tree,
      and it is unique for the node. This is used to lookup nodes in any hash or
      cache table.

      The position of the leftmost bit which is 1 represents the depth of the
      node into the tree, all the bits at its right describe the path of this
      node into the tree, and the bits at its left are clearly 0 and unused.
   */
  // C++11 syntax
  // using NodeKey = KeyType;
  typedef KeyType NodeKey;
  static const int NodeKeyBits = 8*sizeof(NodeKey);

  /// This enumeration determines the different types of node a GenericTreeNode can be
  enum NodeType {
    Invalid = 1,
    Bucket,
    Internal,
    Boundary,
    NonLocal,
    Empty,
    Top,
    NonLocalBucket,
    Cached,
    CachedBucket,
    CachedEmpty
  };

  /// A simple counter for how many nodes are empty - statistics
  extern int numEmptyNodes;

class NodePool;

/// @brief Base class for tree nodes
  class GenericTreeNode {
#ifdef CHANGA_REFACTOR_WALKCHECK
    public:
    bool touched;
    char by;
#endif
  protected:
    NodeType myType;
    NodeKey key;

    GenericTreeNode() : myType(Invalid), key(0), parent(0), firstParticle(0),
	lastParticle(0), remoteIndex(0), iParticleTypes(0), nSPH(0) {
#if COSMO_STATS > 0
      used = false;
#endif
#if INTERLIST_VER > 0
      numBucketsBeneath=0;
      startBucket=-1;
#ifdef CUDA
      nodeArrayIndex = -1;
      bucketArrayIndex = -1;
#endif
#endif
#ifdef CHANGA_REFACTOR_WALKCHECK
      touched = false;
      by = 'I';
#endif
    }

  public:
#if COSMO_STATS > 0
    bool used;
#endif

    /// The moments for the gravity computation
    MultipoleMoments moments;
    /// The parent of this node, or null if none
    GenericTreeNode* parent;
    /// The axis-aligned bounding box of this node
    OrientedBox<cosmoType> boundingBox;
    /// The bounding box including search balls of this node
    OrientedBox<double> bndBoxBall;
    /// Mask of particle types contatained in this node
    unsigned int iParticleTypes;
    /// The number of SPH particles this node contains
    int64_t nSPH;
    /// An index for the first particle contained by this node, 0 means outside the node
    int firstParticle;
    /// An index to the last particle contained by this node, myNumParticles+1 means outside the node
    int lastParticle;
    /// An index to the *real* location of this node if this node is NonLocal, if
    /// it is Boundary or Internal or Bucket it is equal to thisIndex
    /// During Treebuid, it indicates whether remote moments are
    /// needed to calculate this nodes moment.
    int remoteIndex;
    /// Total number of particles contained (across all chares)
    unsigned int particleCount;
    /// Pointer to the first particle in this node
    GravityParticle *particlePointer;

    /// The greatest rung amoung the particles contained in this node (greater
    /// means faster). This information is limited to the nodes in the current
    /// TreePiece, and do not consider non-local data.
    int rungs;

#if INTERLIST_VER > 0
    /// @brief Number of buckets in this node
    int numBucketsBeneath;
    /// @brief index of first bucket in this node
    int startBucket;
#ifdef CUDA
    /// index in moments array sent to GPU
    int nodeArrayIndex;
    int bucketArrayIndex;
#endif
#endif
    /// center of smoothActive particles during smooth operation
    Vector3D<double> centerSm;
    /// Radius of bounding sphere of smoothActive particles
    double sizeSm;
    /// Maximum smoothing radius of smoothActive particles
    double fKeyMax;
    /// SMP rank of node owner
    int iRank;

    /// @brief Construct GenericTreeNode
    /// @param k NodeKey
    /// @param type NodeType
    /// @param first First particle index
    /// @param last  Last particle index
    /// @param p Parent node
    GenericTreeNode(NodeKey k, NodeType type, int first, int last, GenericTreeNode *p) : myType(type), key(k), parent(p), firstParticle(first), lastParticle(last), remoteIndex(0) {
#if INTERLIST_VER > 0
      numBucketsBeneath=0;
      startBucket=-1;
#ifdef CUDA
      nodeArrayIndex = -1;
      bucketArrayIndex = -1;
#endif
#endif
    }

    virtual ~GenericTreeNode() { }
    /// Recursively delete all nodes beneath this node.
    virtual void fullyDelete() = 0;

    /// return Tree::NodeType of node
    inline NodeType getType() const { return myType; }
    /// set Tree::NodeType of node
    inline void setType(NodeType t) { myType = t; }

    /// return unique Tree::NodeKey
    inline NodeKey getKey() const { return key; }

    /// return the number of children this node has
    virtual unsigned int numChildren() const = 0;
    /// return the pointers to the specified child of this node
    virtual GenericTreeNode* getChildren(int) = 0;
    /// set the specified child of this node to the passed pointer
    virtual void setChildren(int, GenericTreeNode*) = 0;
    /// return the keys for the specified child
    virtual NodeKey getChildKey(int) = 0;
    /// return the key for the parent
    virtual NodeKey getParentKey() = 0;
    /// return an integer with the number of the child reflecting the key
    virtual int whichChild(NodeKey childkey) = 0;
#if INTERLIST_VER > 0
    /// Is nodekey contained by this node
    virtual bool contains(NodeKey nodekey) = 0;
#endif

    /// Is the NodeType valid
    bool isValid(){
      return (myType != Invalid);
    }

    /// Is this a node in the cache
    bool isCached(){
      return (myType == Cached ||
              myType == CachedBucket ||
              myType == CachedEmpty);
    }

    /// Is this a node a bucket
    bool isBucket(){
      return (myType == Bucket ||
              myType == CachedBucket ||
              myType == NonLocalBucket);
    }

    /// \brief construct the children of the "this" node following the
    /// given logical criteria (Oct/Orb)
    virtual void makeOctChildren(GravityParticle *part, int totalPart, int level, NodePool *pool = NULL) = 0;
    virtual void makeOrbChildren(GravityParticle *part, int totalPart, int level, int rootsLevel, bool (*compFnPtr[])(GravityParticle, GravityParticle), bool spatial, NodePool *pool = NULL) = 0;

    /// get the top nodes corresponding to a particular number of
    /// chunks requested
    virtual void getChunks(int num, NodeKey *&ret) = 0;

    /// transform an internal node into a bucket
    inline void makeBucket(GravityParticle *part) {
      myType = Bucket;
      iRank = CkMyRank();
#if INTERLIST_VER > 0
      numBucketsBeneath = 1;
#endif
      calculateRadiusBox(moments, boundingBox);	/* set initial size */
      boundingBox.reset();
      bndBoxBall.reset();
      iParticleTypes = 0;
      nSPH = 0;
      rungs = 0;
      for (int i = firstParticle; i <= lastParticle; ++i) {
        moments += part[i];
        boundingBox.grow(part[i].position);
	if(TYPETest(&part[i], TYPE_GAS)) {
	    double fBallMax = part[i].fBallMax();
	    bndBoxBall.grow(part[i].position
			    + Vector3D<double>(fBallMax, fBallMax, fBallMax));
	    bndBoxBall.grow(part[i].position
			    - Vector3D<double>(fBallMax, fBallMax, fBallMax));
            nSPH++;
	    }
	iParticleTypes |= part[i].iType;
        if (part[i].rung > rungs) rungs = part[i].rung;
      }
      if(particleCount > 1)
	  calculateRadiusFarthestParticle(moments, &part[firstParticle],
					  &part[lastParticle+1]);
    }

    /// @brief initialize an empty node
    inline void makeEmpty() {
      myType = Empty;
      particleCount = 0;
#if INTERLIST_VER > 0
      numBucketsBeneath = 0;
#endif
      moments.clear();
      boundingBox.reset();
      bndBoxBall.reset();
      iParticleTypes = 0;
      nSPH = 0;
    }

    /// @brief print out a visualization of the tree for diagnostics
    void getGraphViz(std::ostream &out);

    /// @brief return the NodeKey of the lowest common ancestor.
    virtual NodeKey getLongestCommonPrefix(NodeKey k1, NodeKey k2)
    {
      CkAbort("getLongestCommonPrefix not implemented\n");
      return 0;
    }

    /// @brief depth of node corresponding to NodeKey
    virtual int getLevel(NodeKey k) = 0;
    /// @brief make a copy of the node
    virtual GenericTreeNode *clone() const = 0;

    /// @brief PUP node and children down to depth
    virtual void pup(PUP::er &p, int depth) = 0;
    /// @brief PUP just this node
    virtual void pup(PUP::er &p) {
      int iType;
      if(p.isUnpacking()) {
        p | iType;
  	    myType = (NodeType) iType;
#ifdef CUDA
  	    // so that newly shipped nodes are not mistakenly assumed
  	    // to be present on the GPU
  	    nodeArrayIndex = -1;
      bucketArrayIndex = -1;
#endif
      } else {
        iType = (int) myType;
        p | iType;
      }
      p | key;
      p | moments;
      p | boundingBox;
      p | bndBoxBall;
      p | iParticleTypes;
      p | nSPH;
      p | firstParticle;
      p | lastParticle;
      p | remoteIndex;
      p | particleCount;
#if INTERLIST_VER > 0
      p | numBucketsBeneath;
      p | startBucket;
#endif
      p | centerSm;
      p | sizeSm;
      p | fKeyMax;
#ifdef CHANGA_REFACTOR_WALKCHECK
      p | touched;
      p | by;
#endif
    }
  };

class BinaryTreeNode;

/// Utility to pool allocations of tree nodes
class NodePool {
    int next;
    int szPool;
    /// list of allocated pool blocks.
    std::list<BinaryTreeNode *> pools;
public:
    NodePool() {
        next = szPool = 128;  // Determines the size of the pools
    }
    ~NodePool();
    /// give out one node from the pool, allocating a new pool block if needed.
    BinaryTreeNode *alloc_one();
    BinaryTreeNode *alloc_one(NodeKey k, NodeType type, int first,
			      int nextlast, BinaryTreeNode *p);

    };

  /** A table of the nodes in my tree, indexed by their keys.
      @todo XXX: Make this lookup a hash table, so we get O(1) behavior instead of O(log N).
  */
  typedef std::map<NodeKey, GenericTreeNode *> NodeLookupType;

  /** @brief A TreeNode with two children */
  class BinaryTreeNode : public GenericTreeNode {
  protected:
  public:
    BinaryTreeNode* children[2];

    BinaryTreeNode() : GenericTreeNode() {
      children[0] = 0;
      children[1] = 0;
    }

  public:
    BinaryTreeNode(NodeKey k, NodeType type, int first, int nextlast, BinaryTreeNode *p) : GenericTreeNode(k, type, first, nextlast, p) {
      children[0] = 0;
      children[1] = 0;
    }

    void fullyDelete() {
      if (children[0] != NULL) children[0]->fullyDelete();
      delete children[0];
      if (children[1] != NULL) children[1]->fullyDelete();
      delete children[1];
    }

    unsigned int numChildren() const {
      return 2;
    }

    GenericTreeNode* getChildren(int i) {
      CkAssert(i>=0 && i<2);
      return children[i];
    }

    void setChildren(int i, GenericTreeNode* node) {
      CkAssert(i>=0 && i<2);
      children[i] = (BinaryTreeNode*)node;
    }

    NodeKey getChildKey(int i) {
      CkAssert(i>=0 && i<2);
      return (key<<1) + i;
    }

    NodeKey getParentKey() {
      return (key>>1);
    }

    int getLevel(NodeKey k) {
      int i = 0;
      k >>= 1;
      while (k!=0) {
	k >>= 1;
	++i;
      }
      return i;
    }

    NodeKey getLongestCommonPrefix(NodeKey k1, NodeKey k2){

      int l1 = getLevel(k1)+1;
      int l2 = getLevel(k2)+1;

      NodeKey a1 = k1 << (NodeKeyBits-l1);
      NodeKey a2 = k2 << (NodeKeyBits-l2);

      NodeKey a = a1 ^ a2;
      int i = 0;
      while( a < (NodeKey(1)<<(NodeKeyBits-1)) )
      {
        a = a << 1;
        i++;
      }
      return (k1 >> (l1-i)); // or k2 >> (l2-i)
    };


    int whichChild(NodeKey child) {
      int thisLevel = getLevel(key);
      int childLevel = getLevel(child);
      return ((child >> (childLevel-thisLevel-1)) ^ (key << 1));
    }

#if INTERLIST_VER > 0
    bool contains(NodeKey node) {
      int thisLevel = getLevel(key);
      int nodeLevel = getLevel(node);

      if(nodeLevel<thisLevel) return false;

      return ((node>>(nodeLevel-thisLevel))==key);
    }
#endif

    bool isLeftChild() const {
      return (dynamic_cast<BinaryTreeNode *>(parent) && dynamic_cast<BinaryTreeNode *>(parent)->children[0] == this);
    }

    bool isRightChild() const {
      return (dynamic_cast<BinaryTreeNode *>(parent) && dynamic_cast<BinaryTreeNode *>(parent)->children[1] == this);
    }

    BinaryTreeNode* getSibling() const {
      BinaryTreeNode* p = dynamic_cast<BinaryTreeNode *>(parent);
      if(p)
	return (p->children[0] == this ? p->children[1] : p->children[0]);
      else
	return 0;
    }

    /// Equally divide space into two child nodes.  The split
    /// direction is determined by level.
    /// For each child:
    ///   1. A key is assigned which encodes the tree branching path
    /// to the child
    ///	  2. indices to first and last particles are assigned.
    ///   3. "Type" of node is assigned (Internal, Boundary, etc.)
    ///        This generally depends on the nature of the sibling.
    ///   4. Pointer to particles is assigned.
    ///
    void makeOctChildren(GravityParticle *part, int totalPart, int level,
			 NodePool *pool=NULL) {
      if(pool) {
	  children[0] = pool->alloc_one();
	  children[1] = pool->alloc_one();
	  }
      else {
	  children[0] = new BinaryTreeNode();
	  children[1] = new BinaryTreeNode();
	  }
      children[0]->parent = this;
      children[1]->parent = this;
      children[0]->boundingBox = boundingBox;
      children[1]->boundingBox = boundingBox;
      double split;
      switch (level%3) {
      case 0:
	split = 0.5 * (boundingBox.greater_corner.x + boundingBox.lesser_corner.x);
	children[0]->boundingBox.greater_corner.x = split;
	children[1]->boundingBox.lesser_corner.x = split;
	break;
      case 1:
	split = 0.5 * (boundingBox.greater_corner.y + boundingBox.lesser_corner.y);
	children[0]->boundingBox.greater_corner.y = split;
	children[1]->boundingBox.lesser_corner.y = split;
	break;
      case 2:
	split = 0.5 * (boundingBox.greater_corner.z + boundingBox.lesser_corner.z);
	children[0]->boundingBox.greater_corner.z = split;
	children[1]->boundingBox.lesser_corner.z = split;
	break;
      }
      children[0]->key = key << 1;
      children[1]->key = (key << 1) + 1;
      children[0]->firstParticle = firstParticle;
      children[1]->lastParticle = lastParticle;

      // mask holds the bit that determines a particle's membership in
      // either the left (0) or right (1) child.
      SFC::Key mask = SFC::Key(1) << ((SFC::KeyBits-1) - level);
      SFC::Key leftBit = part[firstParticle].key & mask;
      SFC::Key rightBit = part[lastParticle].key & mask;

      if (leftBit == rightBit) {
	// we have only one child, the other is either NonLocal or Empty
	if (leftBit) {
	  // left child missing
	  if (firstParticle == 0) {
              // We are at the left domain boundary: the left child is
              // completely on another TreePiece
              children[0]->myType = NonLocal;
              children[1]->myType = Boundary;
	  } else {
	    children[0]->makeEmpty();
	    children[1]->myType = lastParticle==totalPart+1 ? Boundary : Internal;
	  }
	  children[0]->lastParticle = firstParticle-1;
          children[0]->particleCount = 0;
	  children[1]->firstParticle = firstParticle;
	  children[1]->particleCount = particleCount;
	} else {
	  // right child missing
	  if (lastParticle == totalPart+1) {
              // We are at the right domain boundary: the right child is
              // completely on another TreePiece
              children[1]->myType = NonLocal;
              children[0]->myType = Boundary;
	  } else {
	    children[1]->makeEmpty();
	    children[0]->myType = firstParticle==0 ? Boundary : Internal;
	  }
	  children[1]->firstParticle = lastParticle+1;
          children[1]->particleCount = 0;
	  children[0]->lastParticle = lastParticle;
	  children[0]->particleCount = particleCount;
	}
      } else if (leftBit < rightBit) {
	// both children are present
	if (firstParticle == 0 && leftBit != (part[1].key & mask)) {
          // the left child is NonLocal: the boundary particle (from
          // another TP) is in the left child, but the rest of the
          // particles are in the right child.
	  children[0]->myType = NonLocal;
	  children[0]->lastParticle = firstParticle;
	  children[1]->myType = lastParticle==totalPart+1 ? Boundary : Internal;
	  children[1]->firstParticle = firstParticle+1;
	  children[1]->particleCount = particleCount;
	} else if (lastParticle == totalPart+1 && rightBit != (part[totalPart].key & mask)) {
          // the right child is NonLocal: the boundary particle (from
          // another TP) is in the right child, but the rest of the
          // particles are in the left child.
	  children[1]->myType = NonLocal;
	  children[1]->firstParticle = lastParticle;
	  children[0]->myType = firstParticle==0 ? Boundary : Internal;
	  children[0]->lastParticle = lastParticle-1;
	  children[0]->particleCount = particleCount;
	} else {
          // the splitting point is somewhere in the middle
          // The following statement finds the first particle with
          // splitting bit (see the variable 'mask' above) set.
	  GravityParticle *splitParticle = std::lower_bound(&part[firstParticle],
			&part[lastParticle+1],
			(GravityParticle)(part[lastParticle].key
					  & ((~ (SFC::Key)0) << ((SFC::KeyBits-1)-level))));
	  children[0]->lastParticle = splitParticle - part - 1;
	  children[1]->firstParticle = splitParticle - part;
	  children[0]->particleCount = children[0]->lastParticle - firstParticle + 1;
	  children[1]->particleCount = lastParticle - children[1]->firstParticle + 1;
	  children[0]->myType = children[0]->particleCount==0 ? Empty : (firstParticle==0 ? Boundary : Internal);
	  children[1]->myType = children[1]->particleCount==0 ? Empty : (lastParticle==totalPart+1 ? Boundary : Internal);
          if (children[0]->myType==Empty) children[0]->makeEmpty();
          if (children[1]->myType==Empty) children[1]->makeEmpty();
	  if (firstParticle == 0) children[0]->particleCount --;
	  if (lastParticle == totalPart+1) children[1]->particleCount --;
	}
      } else {
	CkAbort("Ok, the particles should be ordered so, how the hell did we get here?!?");
      }
      children[0]->particlePointer = &part[children[0]->firstParticle];
      children[1]->particlePointer = &part[children[1]->firstParticle];
#if INTERLIST_VER > 0
      if(children[0]->myType == NonLocal) children[0]->numBucketsBeneath = 0;
      if(children[1]->myType == NonLocal) children[1]->numBucketsBeneath = 0;
      // empty nodes are makeEmpty()'ed, so that the numbucketsbeneath them are 0
#endif
    }

    // Constructs 2 children of this node based on ORB decomposition
    void makeOrbChildren(GravityParticle *part, int totalPart, int level,
		         int rootsLevel,
			 bool (*compFnPtr[])(GravityParticle, GravityParticle),
			 bool spatial, NodePool *pool=NULL) {

      if(pool) {
	  children[0] = pool->alloc_one();
	  children[1] = pool->alloc_one();
	  }
      else {
          children[0] = new BinaryTreeNode();
          children[1] = new BinaryTreeNode();
	  }
      children[0]->parent = this;
      children[1]->parent = this;
      children[0]->boundingBox = boundingBox;
      children[1]->boundingBox = boundingBox;

      children[0]->key = key << 1;
      children[1]->key = (key << 1) + 1;
      children[0]->firstParticle = firstParticle;
      children[1]->lastParticle = lastParticle;

      //CkPrintf("children keys:%lld,%lld\n",children[0]->key,children[1]->key);
      if(level<rootsLevel){
        //This branch is taken for levels above the TreePiece root level
        //TreePiece root level is the level from where onwards data is internal
        SFC::Key tmp = 1 << (level+1);
        tmp = children[0]->key - tmp;
        SFC::Key tmp2 = part[0].key >> (rootsLevel-level-1);

        if(level+1 != rootsLevel){
          if(tmp==tmp2){
            children[0]->myType = Boundary;
            children[1]->myType = NonLocal;
	          children[1]->firstParticle = lastParticle+1;
	          children[0]->lastParticle = lastParticle;
	          children[0]->particleCount = particleCount;
          }
          else{
            children[0]->myType = NonLocal;
            children[1]->myType = Boundary;
	          children[0]->lastParticle = firstParticle-1;
	          children[1]->firstParticle = firstParticle;
	          children[1]->particleCount = particleCount;
          }
        }
        else{ //Special case for TreePiece root level
          if(tmp==tmp2){
            children[0]->myType = Internal;
            children[1]->myType = NonLocal;
	          children[1]->firstParticle = lastParticle+1;
	          children[0]->lastParticle = lastParticle;
	          children[0]->particleCount = particleCount;
	          if(firstParticle==0)
              children[0]->firstParticle = firstParticle+1;
	          if(lastParticle==totalPart+1)
              children[0]->lastParticle = lastParticle-1;
          }
          else{
            children[0]->myType = NonLocal;
            children[1]->myType = Internal;
	          children[0]->lastParticle = firstParticle-1;
	          children[1]->firstParticle = firstParticle;
	          children[1]->particleCount = particleCount;
	          if(firstParticle==0)
              children[1]->firstParticle = firstParticle+1;
	          if(lastParticle==totalPart+1)
              children[1]->lastParticle = lastParticle-1;
          }
        }
      }
      else{ //Below the TreePiece root level
        float len=0.0,len2=0.0;
        int dim;
	
	if(spatial) { // "squeeze" box before division
	    boundingBox.reset();
	    for (int i = firstParticle; i <= lastParticle; ++i)
		boundingBox.grow(part[i].position);
	    }

        len=boundingBox.greater_corner.x-boundingBox.lesser_corner.x;
        dim=0;
        CkAssert(len>=0.0);
        len2=boundingBox.greater_corner.y-boundingBox.lesser_corner.y;
        CkAssert(len2>=0.0);
        if(len2>len) { len = len2; dim=1; }
        len2=boundingBox.greater_corner.z-boundingBox.lesser_corner.z;
        CkAssert(len2>=0.0);
        if(len2>len) { len = len2; dim=2; }

        //Sort the particles in longest dimension
        std::sort(&part[firstParticle],&part[lastParticle+1],compFnPtr[dim]);

        //Middle particle is the splitting point to get 2 children

        //Compress the bounding box of this node
        boundingBox.greater_corner[dim] = part[lastParticle].position[dim];
        boundingBox.lesser_corner[dim] = part[firstParticle].position[dim];

	if(!spatial) {
	    if((lastParticle-firstParticle+1)%2==0){
	      children[0]->lastParticle = firstParticle + (lastParticle-firstParticle-1)/2;
	      children[1]->firstParticle = firstParticle + (lastParticle-firstParticle+1)/2;
	      children[0]->particleCount = (lastParticle-firstParticle+1)/2;
	      children[1]->particleCount = (lastParticle-firstParticle+1)/2;
	    }
	    else{
	      children[0]->lastParticle = firstParticle + (lastParticle-firstParticle)/2;
	      children[1]->firstParticle = firstParticle + (lastParticle-firstParticle+2)/2;
	      children[0]->particleCount = (lastParticle-firstParticle+2)/2;
	      children[1]->particleCount = (lastParticle-firstParticle)/2;
	    }
	}
	else {
	    GravityParticle dummy;  // dummy for splitting
	    GravityParticle* divStart = &part[firstParticle];
	    Vector3D<double> divide(0.0,0.0,0.0);
	    divide[dim] = part[firstParticle].position[dim] + 0.5*len;
	    dummy.position = divide;
	    GravityParticle* divEnd
		= std::upper_bound(&part[firstParticle],&part[lastParticle+1],
				   dummy,compFnPtr[dim]);
	    children[0]->lastParticle = firstParticle + (divEnd - divStart) - 1;
	    children[1]->firstParticle = firstParticle + (divEnd - divStart);
	    children[0]->particleCount = divEnd - divStart;
	    children[1]->particleCount = 1 + (lastParticle-firstParticle)
		- (divEnd - divStart);
	    CkAssert(children[0]->particleCount > 0);
	    CkAssert(children[1]->particleCount > 0);
	}

        children[0]->myType = Internal;
        children[1]->myType = Internal;

        //Compress the bounding boxes of both children too
        children[0]->boundingBox.lesser_corner = boundingBox.lesser_corner;
        children[0]->boundingBox.greater_corner = boundingBox.greater_corner;
        children[0]->boundingBox.greater_corner[dim] = part[children[0]->lastParticle].position[dim];

        children[1]->boundingBox.lesser_corner = boundingBox.lesser_corner;
        children[1]->boundingBox.lesser_corner[dim] = part[children[1]->firstParticle].position[dim];
        children[1]->boundingBox.greater_corner = boundingBox.greater_corner;

      }
      children[0]->particlePointer = &part[children[0]->firstParticle];
      children[1]->particlePointer = &part[children[1]->firstParticle];

    }

    // implemented in the .cpp
    GenericTreeNode *clone() const;

/// @brief Get a number of top level NodeKeys which together make a
/// complete tree.
/// @param num Number of NodeKeys to generate
/// @param ret Array in which to store generated NodeKeys.
///
    void getChunks(int num, NodeKey *&ret) {
      int i = 0;
      int base = num;
      while (base > 1) {
	base >>= 1;
	i++;
      }
      base = 1 << i;
      // base now contains the greatest power of two less or equal to num
      int additional = num - base;
      if (ret==NULL) ret = new NodeKey[num];
      for (int j=additional, k=additional*2; j<base; ++j, ++k) ret[k] = j + base;
      base <<= 1;
      for (int j=0; j<additional*2; ++j) ret[j] = j + base;

    }

    int countDepth(int depth) {
      int count = 1;
      if (depth != 0) {
        if (children[0] != NULL) count += children[0]->countDepth(depth-1);
        if (children[1] != NULL) count += children[1]->countDepth(depth-1);
      }
      return count;
    }

    // extraSpace is used by the CacheInterface to store a pointer to
    // the message so it can be freed when we are finished.
    int packNodes(BinaryTreeNode *buffer, int depth, int extraSpace=0) {
      //CkPrintf("Entering packNodes: this=%p, buffer=%p, depth=%d\n",this,buffer,depth);
      //*buffer = *this;
      memcpy((void *)buffer, (void *)this, sizeof(*this));
      buffer->parent = NULL;
      buffer->particlePointer = NULL;
#if INTERLIST_VER > 0 && defined CUDA
      buffer->nodeArrayIndex = -1;
      buffer->bucketArrayIndex = -1;
#endif
      int used = 1;
      if (depth != 0) {
        if (children[0] != NULL) {
          BinaryTreeNode *nextBuf = (BinaryTreeNode *) (((char*)buffer) + used * ALIGN_DEFAULT(sizeof(BinaryTreeNode)+extraSpace));
          buffer->children[0] = (BinaryTreeNode*)(((char*)nextBuf) - ((char*)buffer));
          //CkPrintf("Entering child 0: offset %ld\n",buffer->children[0]);
          used += children[0]->packNodes(nextBuf, depth-1, extraSpace);
        } else {
          //CkPrintf("Excluding child 0\n");
          buffer->children[0] = NULL;
        }
        if (children[1] != NULL) {
          BinaryTreeNode *nextBuf = (BinaryTreeNode *) (((char*)buffer) + used * ALIGN_DEFAULT(sizeof(BinaryTreeNode)+extraSpace));
          buffer->children[1] = (BinaryTreeNode*)(((char*)nextBuf) - ((char*)buffer));
          //CkPrintf("Entering child 1: offset %ld\n",buffer->children[1]);
          used += children[1]->packNodes(nextBuf, depth-1, extraSpace);
        } else {
          //CkPrintf("Excluding child 1\n");
          buffer->children[1] = NULL;
        }
      } else {
        //CkPrintf("Depth reached\n");
        buffer->children[0] = buffer->children[1] = NULL;
      }
      //CkAssert((long int)buffer->children[0] < 10000);
      //CkAssert((long int)buffer->children[1] < 10000);
      //CkPrintf("Returning used = %d\n",used);
      return used;
    }

    void unpackNodes() {
      if (children[0] != NULL) {
        //CkAssert((long int)children[0] < 10000);
        children[0] = (BinaryTreeNode*)(((long int)children[0]) + ((char*)this));
        children[0]->parent = this;
        children[0]->unpackNodes();
      }
      if (children[1] != NULL) {
        //CkAssert((long int)children[1] < 10000);
        children[1] = (BinaryTreeNode*)(((long int)children[1]) + ((char*)this));
        children[1]->parent = this;
        children[1]->unpackNodes();
      }
    }

    void pup(PUP::er &p) { pup(p, -1); }
    void pup(PUP::er &p, int depth);/* {
      //CkPrintf("Pupper of BinaryTreeNode(%d) called for %s (%d)\n",depth,p.isPacking()?"Packing":p.isUnpacking()?"Unpacking":"Sizing",p.isSizing()?((PUP::sizer*)&p)->size():((PUP::mem*)&p)->size());
      GenericTreeNode::pup(p);
      int isNull;
      for (int i=0; i<2; ++i) {
	isNull = (children[i]==NULL || depth==0) ? 0 : 1;
	p | isNull;
	CkAssert(isNull==0 || isNull==1);
	if (isNull != 0 && depth != 0) {
	  if (p.isUnpacking()) children[i] = new BinaryTreeNode();
	  children[i]->pup(p, depth-1);
	  if (p.isUnpacking()) children[i]->parent = this;
	}
      }
      }*/
  };

inline
NodePool::~NodePool() {
        while(!pools.empty()) {
            delete [] pools.front();
            pools.pop_front();
            }
        }

inline BinaryTreeNode *
NodePool::alloc_one() {
        if(next >= szPool) { // need new pool block
            BinaryTreeNode *nodes = new BinaryTreeNode[szPool];
            next = 0;
            pools.push_front(nodes);
            }
        return &pools.front()[next++];
        }

inline BinaryTreeNode *
NodePool::alloc_one(NodeKey k, NodeType type, int first, int nextlast,
		    BinaryTreeNode *p) {
	BinaryTreeNode *one = alloc_one();
	return new (one) BinaryTreeNode(k, type, first, nextlast, p);
	}

/// Class for Oct tree where each node has 8 direct children.
  class OctTreeNode : public GenericTreeNode {
  protected:
  public:
    /// the 8 different children
    OctTreeNode* children[8];

    OctTreeNode() : GenericTreeNode() {
      for (int i = 0; i < 8; i++) {
        children[i] = 0;
      }
    }

  public:
    OctTreeNode(NodeKey k, NodeType type, int first, int last, BinaryTreeNode *p) : GenericTreeNode(k, type, first, last, p) {
      for (int i = 0; i < 8; i++) {
        children[i] = 0;
      }
    }

    void fullyDelete() {
      for (int i = 0; i < numChildren(); i++) {
        if (children[i] != NULL) {
          children[i]->fullyDelete();
          delete children[i];
        }
      }
    }

    virtual unsigned int numChildren() const {
      return 8;
    }

    virtual GenericTreeNode* getChildren(int i) {
      CkAssert(i>=0 && i<8);
      return children[i];
    }

    virtual void setChildren(int i, GenericTreeNode* node) {
      CkAssert(i>=0 && i<8);
      children[i] = (OctTreeNode*)node;
    }

    NodeKey getChildKey(int i) {
      CkAssert(i>=0 && i<8);
      return (key<<3) + i;
    }

    NodeKey getParentKey() {
      return (key>>3);
    }

    int whichChild(NodeKey child) {
      return (child ^ (key<<3));
    }

    int getLevel(NodeKey k) {
      int i = 0;
      k >>= 3;
      while (k!=0) {
	k >>= 3;
	++i;
      }
      return i;
    }

#if INTERLIST_VER > 0
    bool contains(NodeKey node) {

      return true;
    }
#endif

    void makeOctChildren(GravityParticle *part, int totalPart, int level,
			 NodePool *pool = NULL) {}

    void makeOrbChildren(GravityParticle *part, int totalPart, int level, int rootsLevel, bool (*compFnPtr[])(GravityParticle, GravityParticle), bool spatial,
			 NodePool *pool = NULL) {}

    GenericTreeNode *clone() const {
      OctTreeNode *tmp = new OctTreeNode();
      *tmp = *this;
      return tmp;
    }

    /*
    int getNumChunks(int num) {
      int i = 0;
      while (num > 1) {
	num >>= 3;
	i++;
      }
      return 1 << (3*i);
    }
    */

    void getChunks(int num, NodeKey *&ret) {
      return;
    }

    void pup(PUP::er &p);
    void pup(PUP::er &p, int depth);
  };

  /** Possible types of trees: the first word means the physical structure used
      to implement the tree (Binary or Oct), the second work means the logic
      used to decompose the particles of the space amoung the nodes (Oct or
      ORB). Notice that the SFC domain decomposition builds an oct-tree!
   */
  enum GenericTrees {
    Binary_Oct,
    Oct_Oct,
    Binary_ORB
  };

  /** Added the weight balancer routine*/

  class compare{ //Defines the comparison operator on the map used in balancer
  public:
    compare(){}

    bool operator()(NodeKey key1, NodeKey key2) const {

      NodeKey tmp = NodeKey(1);
      int len1=0, len2=0;
      int cnt=1;

      while(tmp<=key1 || tmp<=key2){
        tmp<<=1;
        cnt++;
        if(len1==0 && tmp>key1){
          len1=cnt-1;
        }
        if(len2==0 && tmp>key2){
          len2=cnt-1;
        }
      }

      if(len1==len2){
        return key1<key2;
      }
      else if(len1>len2){
        key1>>=(len1-len2);
        return key1<key2;
      }
      else{
        key2>>=(len2-len1);
        return key1<key2;
      }
    }
  };

  /// @brief PUP a NodeType
  inline void operator|(PUP::er &p, NodeType &nt) {
    int nti;
    if (p.isUnpacking()) {
      p | nti;
      nt = (NodeType)nti;
    } else {
      nti = (int)nt;
      p | nti;
    }
  }

  /// @brief PUP a tree type.
  inline void operator|(PUP::er &p, GenericTrees &gt) {
    int gti;
    if (p.isUnpacking()) {
      p | gti;
      gt = (GenericTrees)gti;
    } else {
      gti = (int)gt;
      p | gti;
    }
  }

} //close namespace Tree

#endif //GENERICTREENODE_H