tree.hpp

#ifndef TREE_HPP
#define TREE_HPP


#include <assert.h>
#include <typeinfo>
#include <type_traits>
#include <algorithm>
#include <functional>
#include <set>
#include <vector>
#include <deque>
#include <iostream>
#include <random>
#include <cstdint>




#include <hmlp.h>
#include <hmlp_base.hpp>
#include <primitives/combinatorics.hpp>
using namespace std;
using namespace hmlp;

#define REPORT_ANN_STATUS 0

//#define DEBUG_TREE 1


bool has_uneven_split = false;




namespace hmlp
{

class MortonHelper
{
  public:

    typedef pair<size_t, size_t> Recursor;

    static Recursor Root() 
    { 
      return Recursor( 0, 0 ); 
    };

    static Recursor RecurLeft( Recursor r ) 
    { 
      return Recursor( ( r.first << 1 ) + 0, r.second + 1 ); 
    };

    static Recursor RecurRight( Recursor r ) 
    { 
      return Recursor( ( r.first << 1 ) + 1, r.second + 1 ); 
    };

    static size_t MortonID( Recursor r )
    {
      size_t shift = Shift( r.second );
      return ( r.first << shift ) + r.second;
    };

    static size_t SiblingMortonID( Recursor r )
    {
      size_t shift = Shift( r.second );
      if ( r.first % 2 )
        return ( ( r.first - 1 ) << shift ) + r.second;
      else
        return ( ( r.first + 1 ) << shift ) + r.second;
    };

    static int Morton2Rank( size_t it, int size )
    {
      size_t itdepth = Depth( it );
      size_t mpidepth = 0;
      while ( size >>= 1 ) mpidepth ++;
      if ( itdepth > mpidepth ) itdepth = mpidepth;
      size_t itshift = Shift( itdepth );
      return ( it >> itshift ) << ( mpidepth - itdepth );
    }; 
    static void Morton2Offsets( Recursor r, size_t depth, vector<size_t> &offsets )
    {
      if ( r.second == depth ) 
      {
        offsets.push_back( r.first );
      }
      else
      {
        Morton2Offsets( RecurLeft( r ), depth, offsets );
        Morton2Offsets( RecurRight( r ), depth, offsets );
      }
    }; 
    static vector<size_t> Morton2Offsets( size_t me, size_t depth )
    {
      vector<size_t> offsets;
      size_t mydepth = Depth( me );
      assert( mydepth <= depth );
      Recursor r( me >> Shift( mydepth ), mydepth );
      Morton2Offsets( r, depth, offsets );
      return offsets;
    }; 
    static bool IsMyParent( size_t me, size_t it )
    {
      size_t itlevel = Depth( it );
      size_t mylevel = Depth( me );
      size_t itshift = Shift( itlevel );
      bool is_my_parent = !( ( me ^ it ) >> itshift ) && ( itlevel <= mylevel );
    #ifdef DEBUG_TREE
      hmlp_print_binary( me );
      hmlp_print_binary( it );
      hmlp_print_binary( ( me ^ it ) >> itshift );
      printf( "ismyparent %d itlevel %lu mylevel %lu shift %lu fixed shift %d\n",
          is_my_parent, itlevel, mylevel, itshift, 1 << LEVELOFFSET );
    #endif
      return is_my_parent;
    }; 
    template<typename TQUERY>
    static bool ContainAny( size_t target, TQUERY &querys )
    {
      for ( auto & q : querys ) 
        if ( IsMyParent( q, target ) ) return true;
      return false;
    }; 
  private:

    static size_t Depth( size_t it ) 
    {
      size_t filter = ( 1 << LEVELOFFSET ) - 1;
      return it & filter;
    }; 
    static size_t Shift( size_t depth )
    {
      return ( 1 << LEVELOFFSET ) - depth + LEVELOFFSET;
    }; 
    const static int LEVELOFFSET = 4;

}; 
template<typename T>
bool less_first( const pair<T, size_t> &a, const pair<T, size_t> &b )
{
  return ( a.first < b.first );
};
template<typename T>
bool less_second( const pair<T, size_t> &a, const pair<T, size_t> &b )
{
  return ( a.second < b.second );
};
template<typename T>
bool equal_second( const pair<T, size_t> &a, const pair<T, size_t> &b )
{
  return ( a.second == b.second );
};
  
  

template<typename T>
void MergeNeighbors( size_t k, pair<T, size_t> *A, 
    pair<T, size_t> *B, vector<pair<T, size_t>> &aux )
{
  if ( aux.size() != 2 * k ) aux.resize( 2 * k );

  for ( size_t i = 0; i < k; i++ ) aux[     i ] = A[ i ];
  for ( size_t i = 0; i < k; i++ ) aux[ k + i ] = B[ i ];

  sort( aux.begin(), aux.end(), less_second<T> );
  auto it = unique( aux.begin(), aux.end(), equal_second<T> );
  sort( aux.begin(), it, less_first<T> );

  for ( size_t i = 0; i < k; i++ ) A[ i ] = aux[ i ];
}; 
template<typename T>
void MergeNeighbors( size_t k, size_t n,
  vector<pair<T, size_t>> &A, vector<pair<T, size_t>> &B )
{
  assert( A.size() >= n * k && B.size() >= n * k );
   #pragma omp parallel
  {
    vector<pair<T, size_t> > aux( 2 * k );
    #pragma omp for
    for( size_t i = 0; i < n; i++ ) 
    {
      MergeNeighbors( k, &(A[ i * k ]), &(B[ i * k ]), aux );
    }
  }
}; 
namespace tree
{
template<typename NODE>
class IndexPermuteTask : public Task
{
  public:

    NODE *arg;

    void Set( NODE *user_arg )
    {
      name = string( "Permutation" );
      arg = user_arg;
      // Need an accurate cost model.
      cost = 1.0;
    };

    void DependencyAnalysis()
    {
      arg->DependencyAnalysis( RW, this );
      if ( !arg->isleaf )
      {
        arg->lchild->DependencyAnalysis( R, this );
        arg->rchild->DependencyAnalysis( R, this );
      }
      this->TryEnqueue();
    };


    void Execute( Worker* user_worker )
    {
      auto &gids = arg->gids; 
      auto *lchild = arg->lchild;
      auto *rchild = arg->rchild;
      
      if ( !arg->isleaf )
      {
        auto &lgids = lchild->gids;
        auto &rgids = rchild->gids;
        gids = lgids;
        gids.insert( gids.end(), rgids.begin(), rgids.end() );
      }
    };

}; 
template<typename NODE>
class SplitTask : public Task
{
  public:

    NODE *arg = NULL;

    void Set( NODE *user_arg )
    {
      name = string( "Split" );
      arg = user_arg;
      // Need an accurate cost model.
      cost = 1.0;
    };

      void DependencyAnalysis() { arg->DependOnParent( this ); };

    void Execute( Worker* user_worker ) { arg->Split(); };

}; 
// *  @brief This is the default ball tree splitter. Given coordinates,
// *         compute the direction from the two most far away points.
// *         Project all points to this line and split into two groups
// *         using a median select.
// *
// *  @para
// *
// *  @TODO  Need to explit the parallelism.
// */ 
//template<int N_SPLIT, typename T>
//struct centersplit
//{
//  /** closure */
//  Data<T> *Coordinate = NULL;
//
//  inline vector<vector<size_t> > operator()
//  ( 
//    vector<size_t>& gids
//  ) const 
//  {
//    assert( N_SPLIT == 2 );
//
//    Data<T> &X = *Coordinate;
//    size_t d = X.row();
//    size_t n = gids.size();
//
//    T rcx0 = 0.0, rx01 = 0.0;
//    size_t x0, x1;
//    vector<vector<size_t> > split( N_SPLIT );
//
//
//    vector<T> centroid = combinatorics::Mean( d, n, X, gids );
//    vector<T> direction( d );
//    vector<T> projection( n, 0.0 );
//
//    //printf( "After Mean\n" );
//
//    // Compute the farest x0 point from the centroid
//    for ( size_t i = 0; i < n; i ++ )
//    {
//      T rcx = 0.0;
//      for ( size_t p = 0; p < d; p ++ )
//      {
//        T tmp = X( p, gids[ i ] ) - centroid[ p ];
//
//
//        rcx += tmp * tmp;
//      }
//      if ( rcx > rcx0 ) 
//      {
//        rcx0 = rcx;
//        x0 = i;
//      }
//    }
//
//    //printf( "After Farest\n" );
//    //for ( int p = 0; p < d; p ++ )
//    //{
//    //}
//    //printf( "\n" );
//
//    // Compute the farest point x1 from x0
//    for ( size_t i = 0; i < n; i ++ )
//    {
//      T rxx = 0.0;
//      for ( size_t p = 0; p < d; p ++ )
//      {
//          T tmp = X( p, gids[ i ] ) - X( p, gids[ x0 ] );
//        rxx += tmp * tmp;
//      }
//      if ( rxx > rx01 )
//      {
//        rx01 = rxx;
//        x1 = i;
//      }
//    }
//
//
//
//    // Compute direction
//    for ( size_t p = 0; p < d; p ++ )
//    {
//      direction[ p ] = X( p, gids[ x1 ] ) - X( p, gids[ x0 ] );
//    }
//
//    //printf( "After Direction\n" );
//    //for ( int p = 0; p < d; p ++ )
//    //{
//    //  printf( "%5.2lf ", direction[ p ] );
//    //}
//    //printf( "\n" );
//    //exit( 1 );
//
//
//
//    // Compute projection
//    projection.resize( n, 0.0 );
//    for ( size_t i = 0; i < n; i ++ )
//      for ( size_t p = 0; p < d; p ++ )
//        projection[ i ] += X( p, gids[ i ] ) * direction[ p ];
//
//    //printf( "After Projetion\n" );
//    //for ( int p = 0; p < d; p ++ )
//    //{
//    //  printf( "%5.2lf ", projec[ p ] );
//    //}
//    //printf( "\n" );
//
//
//
//    /** Parallel median search */
//    T median;
//    
//    if ( 1 )
//    {
//      median = combinatorics::Select( n, n / 2, projection );
//    }
//    else
//    {
//      auto proj_copy = projection;
//      sort( proj_copy.begin(), proj_copy.end() );
//      median = proj_copy[ n / 2 ];
//    }
//
//
//
//    split[ 0 ].reserve( n / 2 + 1 );
//    split[ 1 ].reserve( n / 2 + 1 );
//
//    /** TODO: Can be parallelized */
//    vector<size_t> middle;
//    for ( size_t i = 0; i < n; i ++ )
//    {
//      if      ( projection[ i ] < median ) split[ 0 ].push_back( i );
//      else if ( projection[ i ] > median ) split[ 1 ].push_back( i );
//      else                                 middle.push_back( i );
//    }
//
//    for ( size_t i = 0; i < middle.size(); i ++ )
//    {
//      if ( split[ 0 ].size() <= split[ 1 ].size() ) split[ 0 ].push_back( middle[ i ] );
//      else                                          split[ 1 ].push_back( middle[ i ] );
//    }
//
//
//    //printf( "split median %lf left %d right %d\n", 
//    //    median,
//    //    (int)split[ 0 ].size(), (int)split[ 1 ].size() );
//
//    //if ( split[ 0 ].size() > 0.6 * n ||
//    //     split[ 1 ].size() > 0.6 * n )
//    //{
//    //  for ( int i = 0; i < n; i ++ )
//    //  {
//    //    printf( "%E ", projection[ i ] );
//    //  } 
//    //  printf( "\n" );
//    //}
//
//
//    return split; 
//  };
//};
//
//
// *  @brief This is the splitter used in the randomized tree. Given
// *         coordinates, project all points onto a random direction
// *         and split into two groups using a median select.
// *
// *  @para
// *
// *  @TODO  Need to explit the parallelism.
// */ 
//template<int N_SPLIT, typename T>
//struct randomsplit
//{
//  /** Closure */
//  Data<T> *Coordinate = NULL;
//
//  inline vector<vector<size_t> > operator()
//  ( 
//    vector<size_t>& gids
//  ) const 
//  {
//    assert( N_SPLIT == 2 );
//
//    Data<T> &X = *Coordinate;
//    size_t d = X.row();
//    size_t n = gids.size();
//
//    vector<vector<size_t> > split( N_SPLIT );
//
//    vector<T> direction( d );
//    vector<T> projection( n, 0.0 );
//
//    // Compute random direction
//    static default_random_engine generator;
//    normal_distribution<T> distribution;
//    for ( int p = 0; p < d; p ++ )
//    {
//      direction[ p ] = distribution( generator );
//    }
//
//    // Compute projection
//    projection.resize( n, 0.0 );
//    for ( size_t i = 0; i < n; i ++ )
//      for ( size_t p = 0; p < d; p ++ )
//        projection[ i ] += X( p, gids[ i ] ) * direction[ p ];
//
//
//    // Parallel median search
//    // T median = Select( n, n / 2, projection );
//    auto proj_copy = projection;
//    sort( proj_copy.begin(), proj_copy.end() );
//    T median = proj_copy[ n / 2 ];
//
//    split[ 0 ].reserve( n / 2 + 1 );
//    split[ 1 ].reserve( n / 2 + 1 );
//
//    /** TODO: Can be parallelized */
//    vector<size_t> middle;
//    for ( size_t i = 0; i < n; i ++ )
//    {
//      if      ( projection[ i ] < median ) split[ 0 ].push_back( i );
//      else if ( projection[ i ] > median ) split[ 1 ].push_back( i );
//      else                                 middle.push_back( i );
//    }
//
//    for ( size_t i = 0; i < middle.size(); i ++ )
//    {
//      if ( split[ 0 ].size() <= split[ 1 ].size() ) split[ 0 ].push_back( middle[ i ] );
//      else                                          split[ 1 ].push_back( middle[ i ] );
//    }
//
//
//    //printf( "split median %lf left %d right %d\n", 
//    //    median,
//    //    (int)split[ 0 ].size(), (int)split[ 1 ].size() );
//
//    //if ( split[ 0 ].size() > 0.6 * n ||
//    //     split[ 1 ].size() > 0.6 * n )
//    //{
//    //  for ( int i = 0; i < n; i ++ )
//    //  {
//    //    printf( "%E ", projection[ i ] );
//    //  } 
//    //  printf( "\n" );
//    //}
//
//
//    return split; 
//  };
//};


//template<typename SETUP, int N_CHILDREN, typename NODEDATA>
template<typename SETUP, typename NODEDATA>
class Node : public ReadWrite
{
  public:

    typedef typename SETUP::T T;
    static const int N_CHILDREN = 2;

    Node( SETUP* setup, size_t n, size_t l, 
        Node *parent, unordered_map<size_t, Node*> *morton2node, Lock *treelock )
    {
      this->setup = setup;
      this->n = n;
      this->l = l;
      this->morton = 0;
      this->treelist_id = 0;
      this->gids.resize( n );
      this->isleaf = false;
      this->parent = parent;
      this->lchild = NULL;
      this->rchild = NULL;
      this->morton2node = morton2node;
      this->treelock = treelock;
      for ( int i = 0; i < N_CHILDREN; i++ ) kids[ i ] = NULL;
    };

    Node( SETUP *setup, int n, int l, vector<size_t> gids,
      Node *parent, unordered_map<size_t, Node*> *morton2node, Lock *treelock )
    {
      this->setup = setup;
      this->n = n;
      this->l = l;
      this->morton = 0;
      this->treelist_id = 0;
      this->gids = gids;
      this->isleaf = false;
      this->parent = parent;
      this->lchild = NULL;
      this->rchild = NULL;
      this->morton2node = morton2node;
      this->treelock = treelock;
      for ( int i = 0; i < N_CHILDREN; i++ ) kids[ i ] = NULL;
    };
  

    Node( size_t morton ) { this->morton = morton; };

    ~Node() {};

    void Resize( int n )
    {
      this->n = n;
      gids.resize( n );
    };


    void Split()
    {
      try
      {
        if ( isleaf ) return;

        int m = setup->m;
        int max_depth = setup->max_depth;

        double beg = omp_get_wtime();
        auto split = setup->splitter( gids );
        double splitter_time = omp_get_wtime() - beg;
        //printf( "splitter %5.3lfs\n", splitter_time );

        if ( std::abs( (int)split[ 0 ].size() - (int)split[ 1 ].size() ) > 1 )
        {
          if ( !has_uneven_split )
          {
            printf( "\n\nWARNING! uneven split. Using random split instead %lu %lu\n\n",
                split[ 0 ].size(), split[ 1 ].size() );
            has_uneven_split = true;
          }
          //printf( "split[ 0 ].size() %lu split[ 1 ].size() %lu\n", 
          //  split[ 0 ].size(), split[ 1 ].size() );
          split[ 0 ].resize( gids.size() / 2 );
          split[ 1 ].resize( gids.size() - ( gids.size() / 2 ) );
          //#pragma omp parallel for
          for ( size_t i = 0; i < gids.size(); i ++ )
          {
            if ( i < gids.size() / 2 ) split[ 0 ][ i ] = i;
            else                       split[ 1 ][ i - ( gids.size() / 2 ) ] = i;
          }
        }

        for ( size_t i = 0; i < N_CHILDREN; i ++ )
        {
          int nchild = split[ i ].size();

          kids[ i ]->Resize( nchild );
          for ( int j = 0; j < nchild; j ++ )
          {
            kids[ i ]->gids[ j ] = gids[ split[ i ][ j ] ];
          }
        }
      }
      catch ( const exception & e )
      {
        cout << e.what() << endl;
      }
    }; 
    bool ContainAny( vector<size_t> &queries )
    {
      if ( !setup->morton.size() )
      {
        printf( "Morton id was not initialized.\n" );
        exit( 1 );
      }
      for ( size_t i = 0; i < queries.size(); i ++ )
      {
        if ( MortonHelper::IsMyParent( setup->morton[ queries[ i ] ], morton ) ) 
        {
#ifdef DEBUG_TREE
          printf( "\n" );
          hmlp_print_binary( setup->morton[ queries[ i ] ] );
          hmlp_print_binary( morton );
          printf( "\n" );
#endif
          return true;
        }
      }
      return false;

    }; 
    bool ContainAny( set<Node*> &querys )
    {
      if ( !setup->morton.size() )
      {
        printf( "Morton id was not initialized.\n" );
        exit( 1 );
      }
      for ( auto it = querys.begin(); it != querys.end(); it ++ )
      {
        if ( MortonHelper::IsMyParent( (*it)->morton, morton ) ) 
        {
          return true;
        }
      }
      return false;

    }; 
    void Print()
    {
      printf( "l %lu offset %lu n %lu\n", this->l, this->offset, this->n );
      hmlp_print_binary( this->morton );
    };


    void DependOnChildren( Task *task )
    {
      if ( this->lchild ) this->lchild->DependencyAnalysis( R, task );
      if ( this->rchild ) this->rchild->DependencyAnalysis( R, task );
      this->DependencyAnalysis( RW, task );
      task->TryEnqueue();
    };

    void DependOnParent( Task *task )
    {
      this->DependencyAnalysis( R, task );
      if ( this->lchild ) this->lchild->DependencyAnalysis( RW, task );
      if ( this->rchild ) this->rchild->DependencyAnalysis( RW, task );
      task->TryEnqueue();
    };

    void DependOnNoOne( Task *task )
    {
      this->DependencyAnalysis( RW, task );
      task->TryEnqueue();
    };


    SETUP *setup = NULL;

    NODEDATA data;

    size_t n;

    size_t l;

    size_t morton = 0;
    size_t offset = 0;

    size_t treelist_id; 

    vector<size_t> gids;

    set<size_t> FarIDs;
    set<Node*>  FarNodes;
    set<size_t> FarNodeMortonIDs;

    set<size_t> NearIDs;
    set<Node*>  NearNodes;
    set<size_t> NearNodeMortonIDs;

    set<size_t> NNFarIDs;
    set<Node*>  NNFarNodes;
    set<Node*>  ProposedNNFarNodes;
    set<size_t> NNFarNodeMortonIDs;

    set<size_t> NNNearIDs;
    set<Node*>  NNNearNodes;
    set<Node*>  ProposedNNNearNodes;
    set<size_t> NNNearNodeMortonIDs;

    //set<size_t> HSSNear;
    //set<size_t> HSSFar;
    //set<size_t> FMMNear;
    //set<size_t> FMMFar;


    vector<map<size_t, Data<T>>> DistFar;
    vector<map<size_t, Data<T>>> DistNear;




    Lock *treelock = NULL;

    Node *kids[ N_CHILDREN ];
    Node *lchild  = NULL; 
    Node *rchild  = NULL;
    Node *sibling = NULL;
    Node *parent  = NULL;
    unordered_map<size_t, Node*> *morton2node = NULL;

    bool isleaf;

  private:

}; 
template<typename SPLITTER, typename DATATYPE>
class Setup
{
  public:

    typedef DATATYPE T;

    Setup() {};

    ~Setup() {};

    size_t m = 0;
    
    size_t max_depth = 15;

    //Data<T> *X = NULL;

    Data<pair<T, size_t>> *NN = NULL;

    vector<size_t> morton;

    SPLITTER splitter;



      vector<size_t> ContainAny( vector<size_t> &queries, size_t target )
    {
         vector<size_t> validation( queries.size(), 0 );

      if ( !morton.size() )
      {
        printf( "Morton id was not initialized.\n" );
        exit( 1 );
      }

      for ( size_t i = 0; i < queries.size(); i ++ )
      {
        //auto it = this->setup->morton.find( queries[ i ] );

            //if ( it != this->setup->morton.end() )
            //{
        //  if ( tree::IsMyParent( *it, this->morton ) ) validation[ i ] = 1;
            //}


       //if ( tree::IsMyParent( morton[ queries[ i ] ], target ) ) 
       if ( MortonHelper::IsMyParent( morton[ queries[ i ] ], target ) ) 
             validation[ i ] = 1;

      }
      return validation;

    }; 
}; 
template<class SETUP, class NODEDATA>
class Tree
{
  public:

    typedef typename SETUP::T T;
    typedef Node<SETUP, NODEDATA> NODE;

    static const int N_CHILDREN = 2;



    SETUP setup;

    size_t n = 0;

    size_t m = 0;

    size_t depth = 0;


    Lock lock;

    vector<NODE*> treelist;

    unordered_map<size_t, NODE*> morton2node;

    Tree() {};

    ~Tree()
    {
      //printf( "~Tree() shared treelist.size() %lu treequeue.size() %lu\n",
      //    treelist.size(), treequeue.size() );
      for ( int i = 0; i < treelist.size(); i ++ )
      {
        if ( treelist[ i ] ) delete treelist[ i ];
      }
      morton2node.clear();
      //printf( "end ~Tree() shared\n" );
    };

    void Offset( NODE *node, size_t offset )
    {
      if ( node )
      {
        node->offset = offset;
        if ( node->lchild )
        {
          Offset( node->lchild, offset + 0 );
          Offset( node->rchild, offset + node->lchild->gids.size() );
        }
      }
    }; 
    void RecursiveMorton( NODE *node, MortonHelper::Recursor r ) 
    {
      if ( !node ) return;
      node->morton = MortonHelper::MortonID( r );
      RecursiveMorton( node->lchild, MortonHelper::RecurLeft( r ) );
      RecursiveMorton( node->rchild, MortonHelper::RecurRight( r ) );
      if ( !node->lchild )
      {
        for ( auto it : node->gids ) setup.morton[ it ] = node->morton;
      }
    };


    void AllocateNodes( NODE *root )
    {
         int glb_depth = std::ceil( std::log2( (double)n / m ) );
         if ( glb_depth > setup.max_depth ) glb_depth = setup.max_depth;
         depth = glb_depth - root->l;

         //printf( "local AllocateNodes n %lu m %lu glb_depth %d loc_depth %lu\n", 
         //    n, m, glb_depth, depth );

      for ( auto node_ptr : treelist ) delete node_ptr;
      treelist.clear();
      morton2node.clear();
      treelist.reserve( 1 << ( depth + 1 ) );
      deque<NODE*> treequeue;
      treequeue.push_back( root );


      while ( auto *node = treequeue.front() )
      {
        node->treelist_id = treelist.size();
        if ( node->l < glb_depth )
        {
          for ( int i = 0; i < N_CHILDREN; i ++ )
          {
            node->kids[ i ] = new NODE( &setup, 0, node->l + 1, node, &morton2node, &lock );
            treequeue.push_back( node->kids[ i ] );
          }
               node->lchild = node->kids[ 0 ];
               node->rchild = node->kids[ 1 ];
          if ( node->lchild ) node->lchild->sibling = node->rchild;
          if ( node->rchild ) node->rchild->sibling = node->lchild;
        }
        else
        {
               node->isleaf = true;
          treequeue.push_back( NULL );
        }
        treelist.push_back( node );
        treequeue.pop_front();
      }

    }; 
    void TreePartition()
    {
      double beg, alloc_time, split_time, morton_time, permute_time;

      this->n = setup.ProblemSize();
      this->m = setup.LeafNodeSize();

      global_indices.clear();
      for ( size_t i = 0; i < n; i ++ ) global_indices.push_back( i );

      has_uneven_split = false;

      beg = omp_get_wtime();
      AllocateNodes( new NODE( &setup, n, 0, global_indices, NULL, &morton2node, &lock ) );
      alloc_time = omp_get_wtime() - beg;

      beg = omp_get_wtime();
      SplitTask<NODE> splittask;
      TraverseDown( splittask );
      ExecuteAllTasks();
      split_time = omp_get_wtime() - beg;


      setup.morton.resize( n );
      RecursiveMorton( treelist[ 0 ], MortonHelper::Root() );


      Offset( treelist[ 0 ], 0 );

      morton2node.clear();
      for ( size_t i = 0; i < treelist.size(); i ++ )
      {
        morton2node[ treelist[ i ]->morton ] = treelist[ i ];
      }

      IndexPermuteTask<NODE> indexpermutetask;
      TraverseUp( indexpermutetask );
      ExecuteAllTasks();

    }; 
    vector<size_t> GetPermutation()
    {
      int n_nodes = 1 << this->depth;
      auto level_beg = this->treelist.begin() + n_nodes - 1;

      vector<size_t> perm;

      for ( int node_ind = 0; node_ind < n_nodes; node_ind ++ )
      {
        auto *node = *(level_beg + node_ind);
        auto gids = node->gids;
        perm.insert( perm.end(), gids.begin(), gids.end() );
      }

      return perm;
    }; 
    template<typename KNNTASK>
    Data<pair<T, size_t>> AllNearestNeighbor( size_t n_tree, size_t k, 
      size_t max_depth, pair<T, size_t> initNN,
      KNNTASK &dummy )
    {
      Data<pair<T, size_t>> NN( k, setup.ProblemSize(), initNN );

      setup.m = 4 * k;
      if ( setup.m < 32 ) setup.m = 32;
      setup.NN = &NN;

      if ( REPORT_ANN_STATUS )
      {
        printf( "========================================================\n");
      }

      for ( int t = 0; t < n_tree; t ++ )      
      {
        double knn_acc = 0.0;
        size_t num_acc = 0;
        TreePartition();
        TraverseLeafs( dummy );
        ExecuteAllTasks();

        size_t n_nodes = 1 << depth;
        auto level_beg = treelist.begin() + n_nodes - 1;
        for ( size_t node_ind = 0; node_ind < n_nodes; node_ind ++ )
        {
          auto *node = *(level_beg + node_ind);
          knn_acc += node->data.knn_acc;
          num_acc += node->data.num_acc;
        }
        if ( REPORT_ANN_STATUS )
        {
          printf( "ANN iter %2d, average accuracy %.2lf%% (over %4lu samples)\n", 
              t, knn_acc / num_acc, num_acc );
        }

        if ( knn_acc / num_acc < 0.8 )
        { 
          if ( 2.0 * setup.m < 2048 ) setup.m = 2.0 * setup.m;
        }
        else break;


#ifdef DEBUG_TREE
        printf( "Iter %2d NN 0 ", t );
        for ( size_t i = 0; i < NN.row(); i ++ )
        {
          printf( "%E(%lu) ", NN[ i ].first, NN[ i ].second );
        }
        printf( "\n" );
#endif
      } 
      if ( REPORT_ANN_STATUS )
      {
        printf( "========================================================\n\n");
      }

      #pragma omp parallel for
      for ( size_t j = 0; j < NN.col(); j ++ )
        sort( NN.data() + j * NN.row(), NN.data() + ( j + 1 ) * NN.row() );

      for ( auto &neig : NN )
      {
        if ( neig.second < 0 || neig.second >= NN.col() )
        {
          printf( "Illegle neighbor gid %lu\n", neig.second );
          break;
        }
      }

      return NN;

    }; 
    Data<int> CheckAllInteractions()
    {
      int total_depth = treelist.back()->l;
      int num_leafs = 1 << total_depth;
      Data<int> A( num_leafs, num_leafs, 0 );
      for ( int t = 1; t < treelist.size(); t ++ )
      {
        auto *node = treelist[ t ];
        for ( auto *it : node->NNNearNodes )
        {
          auto I = MortonHelper::Morton2Offsets( node->morton, total_depth );
          auto J = MortonHelper::Morton2Offsets(   it->morton, total_depth );
          for ( auto i : I ) for ( auto j : J ) A( i, j ) += 1; 
        }
        for ( auto *it : node->NNFarNodes )
        {
          auto I = MortonHelper::Morton2Offsets( node->morton, total_depth );
          auto J = MortonHelper::Morton2Offsets(   it->morton, total_depth );
          for ( auto i : I ) for ( auto j : J ) A( i, j ) += 1; 
        }
      }

      for ( size_t i = 0; i < num_leafs; i ++ )
      {
        for ( size_t j = 0; j < num_leafs; j ++ ) printf( "%d", A( i, j ) );
        printf( "\n" );
      }


      return A;
    }; 
    template<typename TASK, typename... Args>
    void TraverseLeafs( TASK &dummy, Args&... args )
    {
      assert( this->treelist.size() );

      int n_nodes = 1 << this->depth;
      auto level_beg = this->treelist.begin() + n_nodes - 1;

      if ( out_of_order_traversal )
      {
        for ( int node_ind = 0; node_ind < n_nodes; node_ind ++ )
        {
          auto *node = *(level_beg + node_ind);
          RecuTaskSubmit( node, dummy, args... );
        }
      }
      else
      {
        int nthd_glb = omp_get_max_threads();
        #pragma omp parallel for if ( n_nodes > nthd_glb / 2 ) schedule( dynamic )
        for ( int node_ind = 0; node_ind < n_nodes; node_ind ++ )
        {
          auto *node = *(level_beg + node_ind);
          RecuTaskExecute( node, dummy, args... );
        }
      }
    }; 
    template<typename TASK, typename... Args>
    void TraverseUp( TASK &dummy, Args&... args )
    {
      assert( this->treelist.size() );

      int local_begin_level = ( treelist[ 0 ]->l ) ? 1 : 0;

      for ( int l = this->depth; l >= local_begin_level; l -- )
      {
        size_t n_nodes = 1 << l;
        auto level_beg = this->treelist.begin() + n_nodes - 1;


        if ( out_of_order_traversal )
        {
          for ( size_t node_ind = 0; node_ind < n_nodes; node_ind ++ )
          {
            auto *node = *(level_beg + node_ind);
            RecuTaskSubmit( node, dummy, args... );
          }
        }
        else
        {
          int nthd_glb = omp_get_max_threads();
          #pragma omp parallel for if ( n_nodes > nthd_glb / 2 ) schedule( dynamic )
          for ( size_t node_ind = 0; node_ind < n_nodes; node_ind ++ )
          {
            auto *node = *(level_beg + node_ind);
            RecuTaskExecute( node, dummy, args... );
          }
        }
      }
    }; 
    template<typename TASK, typename... Args>
    void TraverseDown( TASK &dummy, Args&... args )
    {
      assert( this->treelist.size() );

      int local_begin_level = ( treelist[ 0 ]->l ) ? 1 : 0;

      for ( int l = local_begin_level; l <= this->depth; l ++ )
      {
        size_t n_nodes = 1 << l;
        auto level_beg = this->treelist.begin() + n_nodes - 1;

        if ( out_of_order_traversal )
        {
          for ( size_t node_ind = 0; node_ind < n_nodes; node_ind ++ )
          {
            auto *node = *(level_beg + node_ind);
            RecuTaskSubmit( node, dummy, args... );
          }
        }
        else
        {
          int nthd_glb = omp_get_max_threads();
          #pragma omp parallel for if ( n_nodes > nthd_glb / 2 ) schedule( dynamic )
          for ( size_t node_ind = 0; node_ind < n_nodes; node_ind ++ )
          {
            auto *node = *(level_beg + node_ind);
            RecuTaskExecute( node, dummy, args... );
          }
        }
      }
    }; 
    template<typename TASK, typename... Args>
    void TraverseUnOrdered( TASK &dummy, Args&... args )
    {
      TraverseDown( dummy, args... );
    }; 
    void DependencyCleanUp()
    {
      //for ( size_t i = 0; i < treelist.size(); i ++ )
      //{
      //  treelist[ i ]->DependencyCleanUp();
      //}
      for ( auto node : treelist ) node->DependencyCleanUp();

      for ( auto it : morton2node )
      {
        auto *node = it.second;
        if ( node ) node->DependencyCleanUp();
      }
    }; 
    void ExecuteAllTasks()
    {
      hmlp_run();
      DependencyCleanUp();
    }; 
    bool DoOutOfOrder() { return out_of_order_traversal; };


    template<typename SUMMARY>
    void Summary( SUMMARY &summary )
    {
      assert( N_CHILDREN == 2 );

      for ( std::size_t l = 0; l <= depth; l ++ )
      {
        size_t n_nodes = 1 << l;
        auto level_beg = treelist.begin() + n_nodes - 1;
        for ( size_t node_ind = 0; node_ind < n_nodes; node_ind ++ )
        {
          auto *node = *(level_beg + node_ind);
          summary( node );
        }
      }
    }; 
  private:

    bool out_of_order_traversal = true;

  protected:

    vector<size_t> global_indices;

}; 
}; 
}; 
#endif