]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STAT/TKDTree.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
CPU and Memory tests, coding violations fixed, library
[u/mrichter/AliRoot.git] / STAT / TKDTree.cxx
diff --git a/STAT/TKDTree.cxx b/STAT/TKDTree.cxx
index e679ed3cf67eda7eccc16df044e05620466f2c9a..aba2ba7ba8c39e42eccdd4b77e31273bd96a5c29 100644 (file)
--- a/STAT/TKDTree.cxx
+++ b/STAT/TKDTree.cxx
@@ -29,6 +29,16 @@ templateClassImp(TKDTree)
 //      - after are the terminal nodes
 //     fNodes : array of primary nodes
 ///////////////////////////////////////////////////////////////////////
 //      - after are the terminal nodes
 //     fNodes : array of primary nodes
 ///////////////////////////////////////////////////////////////////////
+
+#include "malloc.h"
+int memory()
+{
+       static struct mallinfo debug;
+       debug = mallinfo();
+       //printf("arena[%d] usmblks[%d] uordblks[%d] fordblks[%d]\n", debug.arena, debug.usmblks, debug.uordblks, debug.fordblks);
+       return debug.uordblks;
+}
+
 //_________________________________________________________________
 template <typename  Index, typename Value>
 TKDTree<Index, Value>::TKDTree() :
 //_________________________________________________________________
 template <typename  Index, typename Value>
 TKDTree<Index, Value>::TKDTree() :
@@ -36,13 +46,19 @@ TKDTree<Index, Value>::TKDTree() :
        ,kDataOwner(kFALSE)
        ,fNnodes(0)
        ,fNDim(0)
        ,kDataOwner(kFALSE)
        ,fNnodes(0)
        ,fNDim(0)
+       ,fNDimm(0)
        ,fNpoints(0)
        ,fBucketSize(0)
        ,fNpoints(0)
        ,fBucketSize(0)
-       ,fData(0x0)
+       ,fAxis(0x0)
+       ,fValue(0x0)
        ,fRange(0x0)
        ,fRange(0x0)
+       ,fData(0x0)
        ,fBoundaries(0x0)
        ,fBoundaries(0x0)
-       ,fNodes(0x0)
+       ,fkNNdim(0)
        ,fkNN(0x0)
        ,fkNN(0x0)
+       ,fkNNdist(0x0)
+       ,fDistBuffer(0x0)
+       ,fIndBuffer(0x0)
        ,fIndPoints(0x0)
        ,fRowT0(0)
        ,fCrossNode(0)
        ,fIndPoints(0x0)
        ,fRowT0(0)
        ,fCrossNode(0)
@@ -59,29 +75,41 @@ TKDTree<Index, Value>::TKDTree(Index npoints, Index ndim, UInt_t bsize, Value **
        ,kDataOwner(kFALSE)
        ,fNnodes(0)
        ,fNDim(ndim)
        ,kDataOwner(kFALSE)
        ,fNnodes(0)
        ,fNDim(ndim)
+       ,fNDimm(2*ndim)
        ,fNpoints(npoints)
        ,fBucketSize(bsize)
        ,fNpoints(npoints)
        ,fBucketSize(bsize)
-       ,fData(data) //Columnwise!!!!!
+       ,fAxis(0x0)
+       ,fValue(0x0)
        ,fRange(0x0)
        ,fRange(0x0)
+       ,fData(data) //Columnwise!!!!!
        ,fBoundaries(0x0)
        ,fBoundaries(0x0)
-       ,fNodes(0x0)
+       ,fkNNdim(0)
        ,fkNN(0x0)
        ,fkNN(0x0)
+       ,fkNNdist(0x0)
+       ,fDistBuffer(0x0)
+       ,fIndBuffer(0x0)
        ,fIndPoints(0x0)
        ,fRowT0(0)
        ,fCrossNode(0)
        ,fOffset(0)
 {
        ,fIndPoints(0x0)
        ,fRowT0(0)
        ,fCrossNode(0)
        ,fOffset(0)
 {
-       // Allocate data by hand. See TKDTree(TTree*, const Char_t*) for an automatic method.
-       
+// Allocate data by hand. See TKDTree(TTree*, const Char_t*) for an automatic method.
+
+       //printf("start building kdTree %d\n", memory());
        Build();
        Build();
+       //printf("finish building kdTree %d\n", memory());
 }
 
 //_________________________________________________________________
 template <typename  Index, typename Value>
 TKDTree<Index, Value>::~TKDTree()
 {
 }
 
 //_________________________________________________________________
 template <typename  Index, typename Value>
 TKDTree<Index, Value>::~TKDTree()
 {
+       if (fIndBuffer) delete [] fIndBuffer;
+       if (fDistBuffer) delete [] fDistBuffer;
+       if (fkNNdist) delete [] fkNNdist;
        if (fkNN) delete [] fkNN;
        if (fkNN) delete [] fkNN;
-       if (fNodes) delete [] fNodes;
+       if (fAxis) delete [] fAxis;
+       if (fValue) delete [] fValue;
        if (fIndPoints) delete [] fIndPoints;
        if (fRange) delete [] fRange;
        if (fBoundaries) delete [] fBoundaries;
        if (fIndPoints) delete [] fIndPoints;
        if (fRange) delete [] fRange;
        if (fBoundaries) delete [] fBoundaries;
@@ -124,7 +152,8 @@ void TKDTree<Index, Value>::Build(){
        fRange = new Value[2*fNDim];
        fIndPoints= new Index[fNpoints];
        for (Index i=0; i<fNpoints; i++) fIndPoints[i] = i;
        fRange = new Value[2*fNDim];
        fIndPoints= new Index[fNpoints];
        for (Index i=0; i<fNpoints; i++) fIndPoints[i] = i;
-       fNodes  = new TKDNode[fNnodes];
+       fAxis  = new UChar_t[fNnodes];
+       fValue = new Value[fNnodes];
        //
        fCrossNode = (1<<(fRowT0+1))-1;
        if (fCrossNode<fNnodes) fCrossNode = 2*fCrossNode+1;
        //
        fCrossNode = (1<<(fRowT0+1))-1;
        if (fCrossNode<fNnodes) fCrossNode = 2*fCrossNode+1;
@@ -166,7 +195,6 @@ void TKDTree<Index, Value>::Build(){
                Int_t crow     = rowStack[currentIndex];
                Int_t cpos     = posStack[currentIndex];
                Int_t cnode    = nodeStack[currentIndex];               
                Int_t crow     = rowStack[currentIndex];
                Int_t cpos     = posStack[currentIndex];
                Int_t cnode    = nodeStack[currentIndex];               
-               TKDNode * node = &(fNodes[cnode]);
                //printf("currentIndex %d npoints %d node %d\n", currentIndex, npoints, cnode);
                //
                // divide points
                //printf("currentIndex %d npoints %d node %d\n", currentIndex, npoints, cnode);
                //
                // divide points
@@ -201,12 +229,14 @@ void TKDTree<Index, Value>::Build(){
                                maxspread=tempspread;
                                axspread = idim;
                        }
                                maxspread=tempspread;
                                axspread = idim;
                        }
-                       if(cnode==0) {fRange[2*idim] = min; fRange[2*idim+1] = max;}
+                       if(cnode) continue;
+                       //printf("set %d %6.3f %6.3f\n", idim, min, max);
+                       fRange[2*idim] = min; fRange[2*idim+1] = max;
                }
                array = fData[axspread];
                KOrdStat(npoints, array, nleft, fIndPoints+cpos);
                }
                array = fData[axspread];
                KOrdStat(npoints, array, nleft, fIndPoints+cpos);
-               node->fAxis  = axspread;
-               node->fValue = array[fIndPoints[cpos+nleft]];
+               fAxis[cnode]  = axspread;
+               fValue[cnode] = array[fIndPoints[cpos+nleft]];
                //printf("Set node %d : ax %d val %f\n", cnode, node->fAxis, node->fValue);
                //
                //
                //printf("Set node %d : ax %d val %f\n", cnode, node->fAxis, node->fValue);
                //
                //
@@ -235,7 +265,7 @@ void TKDTree<Index, Value>::Build(){
 
 //_________________________________________________________________
 template <typename  Index, typename Value>
 
 //_________________________________________________________________
 template <typename  Index, typename Value>
-Bool_t TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int_t k, Index *&in, Value &d)
+Bool_t TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int_t k, Index *&in, Value *&d)
 {
 // Find "k" nearest neighbors to "point".
 //
 {
 // Find "k" nearest neighbors to "point".
 //
@@ -244,7 +274,7 @@ Bool_t      TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int
 // increasing order of the distance to "point" and the maxim distance
 // in "d".
 //
 // increasing order of the distance to "point" and the maxim distance
 // in "d".
 //
-// The array "in" is managed internally by the TKDTree.
+// The arrays "in" and "d" are managed internally by the TKDTree.
 
        Index inode = FindNode(point);
        if(inode < fNnodes){
 
        Index inode = FindNode(point);
        if(inode < fNnodes){
@@ -252,28 +282,46 @@ Bool_t    TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int
                return kFALSE;
        }
 
                return kFALSE;
        }
 
+       UInt_t debug = 0;
+       Int_t nCalculations = 0;
+       
        // allocate working memory
        // allocate working memory
-       if(!fkNN) fkNN = new Index[k];
-       if(sizeof(fkNN) < k*sizeof(Index)) fkNN = new(fkNN) Index[k];
-       for(int i=0; i<k; i++) fkNN[i] = -1;    
-       Index *fkNN_tmp = new Index[k];
-       Value *fkNN_dist = new Value[k];
-       for(int i=0; i<k; i++) fkNN_dist[i] = 9999.;
-       Value *fkNN_dist_tmp = new Value[k];
-       Value *dist = new Value[fBucketSize];
-       Index *index = new Index[fBucketSize];
-
+       if(!fDistBuffer){
+               fDistBuffer = new Value[fBucketSize];
+               fIndBuffer  = new Index[fBucketSize];
+       }
+       if(fkNNdim < k){
+               if(debug>=1){
+                       if(fkNN) printf("Reallocate memory %d -> %d \n", fkNNdim, 2*k);
+                       else printf("Allocate %d memory\n", 2*k);
+               }
+               fkNNdim  = 2*k;
+               if(fkNN){
+                       delete [] fkNN; 
+                       delete [] fkNNdist;
+               }
+               fkNN     = new Index[fkNNdim];
+               fkNNdist = new Value[fkNNdim];
+       }
+       memset(fkNN, -1, k*sizeof(Index));
+       for(int i=0; i<k; i++) fkNNdist[i] = 9999.;
+       Index itmp, jtmp; Value ftmp, gtmp;
+       
        // calculate number of boundaries with the data domain. 
        // calculate number of boundaries with the data domain. 
-       Index nBounds = 0;
        if(!fBoundaries) MakeBoundaries();
        Value *bounds = &fBoundaries[inode*2*fNDim];
        if(!fBoundaries) MakeBoundaries();
        Value *bounds = &fBoundaries[inode*2*fNDim];
+       Index nBounds = 0;
        for(int idim=0; idim<fNDim; idim++){
        for(int idim=0; idim<fNDim; idim++){
-               if(bounds[2*idim] == fRange[2*idim]) nBounds++;
-               if(bounds[2*idim+1] == fRange[2*idim+1]) nBounds++;
+               if((bounds[2*idim] - fRange[2*idim]) < 1.E-10) nBounds++;
+               if((bounds[2*idim+1] - fRange[2*idim+1]) < 1.E-10) nBounds++;
        }
        }
-       
+       if(debug>=1) printf("Calculate boundaries [nBounds = %d]\n", nBounds);
+
+               
        // traverse tree
        // traverse tree
-       TKDNode *node;
+       UChar_t ax;
+       Value   val, dif;
+       Int_t nAllNodes = fNnodes + GetNTNodes();
        Int_t nodeStack[128], nodeIn[128];
        Index currentIndex = 0;
        nodeStack[0]   = inode;
        Int_t nodeStack[128], nodeIn[128];
        Index currentIndex = 0;
        nodeStack[0]   = inode;
@@ -282,14 +330,16 @@ Bool_t    TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int
                Int_t tnode = nodeStack[currentIndex];
                Int_t entry = nodeIn[currentIndex];
                currentIndex--;
                Int_t tnode = nodeStack[currentIndex];
                Int_t entry = nodeIn[currentIndex];
                currentIndex--;
-
+               if(debug>=1) printf("Analyse tnode %d entry %d\n", tnode, entry);
+               
                // check if node is still eligible
                // check if node is still eligible
-               node = &fNodes[tnode/2 + (tnode%2) - 1];
-               if((TMath::Abs(point[node->fAxis] - node->fValue) > fkNN_dist[k-1])
-                                                                                       &&
-                       ((point[node->fAxis] > node->fValue && tnode%2) ||
-                        (point[node->fAxis] < node->fValue && !tnode%2))) {
-                       //printf("\tREMOVE NODE\n");
+               Int_t inode = (tnode-1)>>1; //tnode/2 + (tnode%2) - 1;
+               ax = fAxis[inode];
+               val = fValue[inode];
+               dif = point[ax] - val;
+               if((TMath::Abs(dif) > fkNNdist[k-1]) &&
+                       ((dif > 0. && (tnode&1)) || (dif < 0. && !(tnode&1)))) {
+                       if(debug>=1) printf("\tremove %d\n", (tnode-1)>>1);
 
                        // mark bound
                        nBounds++;
 
                        // mark bound
                        nBounds++;
@@ -299,85 +349,91 @@ Bool_t    TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int
                }
                
                if(IsTerminal(tnode)){
                }
                
                if(IsTerminal(tnode)){
+                       if(debug>=2) printf("\tProcess terminal node %d\n", tnode);
                        // Link data to terminal node
                        Int_t offset = (tnode >= fCrossNode) ? (tnode-fCrossNode)*fBucketSize : fOffset+(tnode-fNnodes)*fBucketSize;
                        Index *indexPoints = &fIndPoints[offset];
                        // Link data to terminal node
                        Int_t offset = (tnode >= fCrossNode) ? (tnode-fCrossNode)*fBucketSize : fOffset+(tnode-fNnodes)*fBucketSize;
                        Index *indexPoints = &fIndPoints[offset];
-                       Int_t nPoints = (tnode == 2*fNnodes) ? fNpoints%fBucketSize : fBucketSize;
-                       for(int idx=0; idx<nPoints; idx++){
+                       Int_t nbs = (tnode == 2*fNnodes) ? fNpoints%fBucketSize : fBucketSize;
+                       nbs = nbs ? nbs : fBucketSize;
+                       nCalculations += nbs;
+                       
+                       Int_t npoints = 0;
+                       for(int idx=0; idx<nbs; idx++){
                                // calculate distance in the L1 metric
                                // calculate distance in the L1 metric
-                               dist[idx] = 0.;
-                               for(int idim=0; idim<fNDim; idim++) dist[idx] += TMath::Abs(point[idim] - fData[idim][indexPoints[idx]]);
+                               fDistBuffer[npoints] = 0.;
+                               for(int idim=0; idim<fNDim; idim++) fDistBuffer[npoints] += TMath::Abs(point[idim] - fData[idim][indexPoints[idx]]);
+                               // register candidate neighbor
+                               if(fDistBuffer[npoints] < fkNNdist[k-1]){
+                                       fIndBuffer[npoints] = indexPoints[idx];
+                                       npoints++;
+                               }
                        }
                        }
-                       // arrange increasingly distances array and update kNN indexes
-                       TMath::Sort(nPoints, dist, index, kFALSE);
-                       for(int ibrowse=0; ibrowse<nPoints; ibrowse++){
-                               // exit if the current distance calculated in this node is
-                               // larger than the largest distance stored in the kNN array
-                               if(dist[index[ibrowse]] >= fkNN_dist[k-1]) break;
-                               for(int iNN=0; iNN<k; iNN++){
-                                       if(dist[index[ibrowse]] < fkNN_dist[iNN]){
-                                               //insert neighbor. In principle this is only one call to STL vector. Maybe we can change to it ?!
-                                               //printf("\t\tinsert data %d @ %d distance %f\n", indexPoints[index[ibrowse]], iNN, dist[index[ibrowse]]);
-                       
-                                               memcpy(fkNN_tmp, &fkNN[iNN], (k-iNN-1)*sizeof(Index));
-                                               fkNN[iNN] = indexPoints[index[ibrowse]];
-                                               memcpy(&fkNN[iNN+1], fkNN_tmp, (k-iNN-1)*sizeof(Index));
-
-                                               memcpy(fkNN_dist_tmp, &fkNN_dist[iNN], (k-iNN-1)*sizeof(Value));
-                                               fkNN_dist[iNN] = dist[index[ibrowse]];
-                                               memcpy(&fkNN_dist[iNN+1], fkNN_dist_tmp, (k-iNN-1)*sizeof(Value));
-                                               break;
-                                       }
+                       for(int ibrowse=0; ibrowse<npoints; ibrowse++){
+                               if(fDistBuffer[ibrowse] >= fkNNdist[k-1]) continue;
+                               //insert neighbor
+                               int iNN=0;
+                               while(fDistBuffer[ibrowse] >= fkNNdist[iNN]) iNN++;
+                               if(debug>=2) printf("\t\tinsert data %d @ %d distance %f\n", fIndBuffer[ibrowse], iNN, fDistBuffer[ibrowse]);
+                               
+                               itmp = fkNN[iNN]; ftmp = fkNNdist[iNN];
+                               fkNN[iNN]     = fIndBuffer[ibrowse];
+                               fkNNdist[iNN] = fDistBuffer[ibrowse];
+                               for(int ireplace=iNN+1; ireplace<k; ireplace++){
+                                       jtmp = fkNN[ireplace]; gtmp = fkNNdist[ireplace];
+                                       fkNN[ireplace]     = itmp; fkNNdist[ireplace] = ftmp;
+                                       itmp = jtmp; ftmp = gtmp;
+                                       if(ftmp == 9999.) break;
                                }
                                }
+                               if(debug>=3){
+                                       for(int i=0; i<k; i++){
+                                               if(fkNNdist[i] == 9999.) break;
+                                               printf("\t\tfkNNdist[%d] = %f\n", i, fkNNdist[i]);
+                                       }
+                               }                               
                        }
                }
                
                // register parent
                        }
                }
                
                // register parent
-               Int_t parent_node = tnode/2 + (tnode%2) - 1;
-               if(parent_node >= 0 && entry != parent_node){
+               Int_t pn = (tnode-1)>>1;//tnode/2 + (tnode%2) - 1;
+               if(pn >= 0 && entry != pn){
                        // check if parent node is eligible at all
                        // check if parent node is eligible at all
-                       node = &fNodes[parent_node];
-                       if(TMath::Abs(point[node->fAxis] - node->fValue) > fkNN_dist[k-1]){
+                       if(TMath::Abs(point[fAxis[pn]] - fValue[pn]) > fkNNdist[k-1]){
                                // mark bound
                                nBounds++;
                                // end all recursions
                                if(nBounds==2 * fNDim) break;
                        }
                        currentIndex++;
                                // mark bound
                                nBounds++;
                                // end all recursions
                                if(nBounds==2 * fNDim) break;
                        }
                        currentIndex++;
-                       nodeStack[currentIndex]=tnode/2 + (tnode%2) - 1;
+                       nodeStack[currentIndex]=pn;
                        nodeIn[currentIndex]=tnode;
                        nodeIn[currentIndex]=tnode;
-                       //printf("\tregister %d\n", nodeStack[currentIndex]);
+                       if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]);
                }
                }
-
+               if(IsTerminal(tnode)) continue;
+               
                // register children nodes
                // register children nodes
-               Int_t child_node;
+               Int_t cn;
                Bool_t kAllow[] = {kTRUE, kTRUE};
                Bool_t kAllow[] = {kTRUE, kTRUE};
-               node = &fNodes[tnode];
-               if(TMath::Abs(point[node->fAxis] - node->fValue) > fkNN_dist[k-1]){
-                       if(point[node->fAxis] > node->fValue) kAllow[0] = kFALSE;
+               ax = fAxis[tnode];
+               val = fValue[tnode];
+               if(TMath::Abs(point[ax] - val) > fkNNdist[k-1]){
+                       if(point[ax] > val) kAllow[0] = kFALSE;
                        else kAllow[1] = kFALSE;
                }
                for(int ic=1; ic<=2; ic++){
                        if(!kAllow[ic-1]) continue;
                        else kAllow[1] = kFALSE;
                }
                for(int ic=1; ic<=2; ic++){
                        if(!kAllow[ic-1]) continue;
-                       child_node = (tnode*2)+ic;
-                       if(child_node < fNnodes + GetNTerminalNodes() && entry != child_node){
+                       cn = (tnode<<1)+ic;
+                       if(cn < nAllNodes && entry != cn){
                                currentIndex++;
                                currentIndex++;
-                               nodeStack[currentIndex] = child_node;
+                               nodeStack[currentIndex] = cn;
                                nodeIn[currentIndex]=tnode;
                                nodeIn[currentIndex]=tnode;
-                               //printf("\tregister %d\n", nodeStack[currentIndex]);
+                               if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]);
                        }
                        }
-               }
+               }               
        }
        // save results
        in = fkNN;
        }
        // save results
        in = fkNN;
-       d  = fkNN_dist[k-1];
+       d  = fkNNdist;
        
        
-       delete [] index;
-       delete [] dist;
-       delete [] fkNN_tmp;
-       delete [] fkNN_dist;
-       delete [] fkNN_dist_tmp;
-
        return kTRUE;
 }
 
        return kTRUE;
 }
 
@@ -388,24 +444,22 @@ template <typename  Index, typename Value>
 Index TKDTree<Index, Value>::FindNode(const Value * point){
   //
   // find the terminal node to which point belongs
 Index TKDTree<Index, Value>::FindNode(const Value * point){
   //
   // find the terminal node to which point belongs
-  
+
        Index stackNode[128], inode;
        Int_t currentIndex =0;
        stackNode[0] = 0;
        Index stackNode[128], inode;
        Int_t currentIndex =0;
        stackNode[0] = 0;
-       //
        while (currentIndex>=0){
                inode    = stackNode[currentIndex];
        while (currentIndex>=0){
                inode    = stackNode[currentIndex];
-               currentIndex--;
                if (IsTerminal(inode)) return inode;
                
                if (IsTerminal(inode)) return inode;
                
-               TKDNode & node = fNodes[inode];
-               if (point[node.fAxis]<=node.fValue){
+               currentIndex--;
+               if (point[fAxis[inode]]<=fValue[inode]){
                        currentIndex++;
                        currentIndex++;
-                       stackNode[currentIndex]=(inode*2)+1;
+                       stackNode[currentIndex]=(inode<<1)+1;
                }
                }
-               if (point[node.fAxis]>=node.fValue){
+               if (point[fAxis[inode]]>=fValue[inode]){
                        currentIndex++;
                        currentIndex++;
-                       stackNode[currentIndex]=(inode*2)+2;
+                       stackNode[currentIndex]=(inode+1)<<1;
                }
        }
        
                }
        }
        
@@ -446,12 +500,11 @@ void TKDTree<Index, Value>::FindPoint(Value * point, Index &index, Int_t &iter){
                        continue;
                }
                
                        continue;
                }
                
-               TKDNode & node = fNodes[inode];
-               if (point[node.fAxis]<=node.fValue){
+               if (point[fAxis[inode]]<=fValue[inode]){
                        currentIndex++;
                        stackNode[currentIndex]=(inode*2)+1;
                }
                        currentIndex++;
                        stackNode[currentIndex]=(inode*2)+1;
                }
-               if (point[node.fAxis]>=node.fValue){
+               if (point[fAxis[inode]]>=fValue[inode]){
                        currentIndex++;
                        stackNode[currentIndex]=(inode*2)+2;
                }
                        currentIndex++;
                        stackNode[currentIndex]=(inode*2)+2;
                }
@@ -485,12 +538,12 @@ void TKDTree<Index, Value>::FindInRangeA(Value * point, Value * delta, Index *re
     currentIndex--;
     if (!IsTerminal(inode)){
       // not terminal
     currentIndex--;
     if (!IsTerminal(inode)){
       // not terminal
-      TKDNode * node = &(fNodes[inode]);
-      if (point[node->fAxis] - delta[node->fAxis] <  node->fValue) {
+      //TKDNode<Index, Value> * node = &(fNodes[inode]);
+      if (point[fAxis[inode]] - delta[fAxis[inode]] <  fValue[inode]) {
        currentIndex++; 
        stackNode[currentIndex]= (inode*2)+1;
       }
        currentIndex++; 
        stackNode[currentIndex]= (inode*2)+1;
       }
-      if (point[node->fAxis] + delta[node->fAxis] >= node->fValue){
+      if (point[fAxis[inode]] + delta[fAxis[inode]] >= fValue[inode]){
        currentIndex++; 
        stackNode[currentIndex]= (inode*2)+2;
       }
        currentIndex++; 
        stackNode[currentIndex]= (inode*2)+2;
       }
@@ -573,18 +626,18 @@ void TKDTree<Index, Value>::FindInRangeB(Value * point, Value * delta, Index *re
       }
     }
     //
       }
     }
     //
-    TKDNode * node = &(fNodes[inode]);
-    if (point[node->fAxis] - delta[node->fAxis] <  node->fValue) {
+    // TKDNode<Index, Value> * node = &(fNodes[inode]);
+    if (point[fAxis[inode]] - delta[fAxis[inode]] <  fValue[inode]) {
       currentIndex++; 
       stackNode[currentIndex]= (inode<<1)+1;
       currentIndex++; 
       stackNode[currentIndex]= (inode<<1)+1;
-      if (point[node->fAxis] + delta[node->fAxis] > node->fValue) 
-       stackStatus[currentIndex]= status | (1<<(2*node->fAxis));
+      if (point[fAxis[inode]] + delta[fAxis[inode]] > fValue[inode])
+       stackStatus[currentIndex]= status | (1<<(2*fAxis[inode]));
     }
     }
-    if (point[node->fAxis] + delta[node->fAxis] >= node->fValue){
+    if (point[fAxis[inode]] + delta[fAxis[inode]] >= fValue[inode]){
       currentIndex++; 
       stackNode[currentIndex]= (inode<<1)+2;
       currentIndex++; 
       stackNode[currentIndex]= (inode<<1)+2;
-      if (point[node->fAxis] - delta[node->fAxis]<node->fValue) 
-       stackStatus[currentIndex]= status | (1<<(2*node->fAxis+1));
+      if (point[fAxis[inode]] - delta[fAxis[inode]]<fValue[inode])
+       stackStatus[currentIndex]= status | (1<<(2*fAxis[inode]+1));
     }
   }
   if (fData){
     }
   }
   if (fData){
@@ -700,62 +753,69 @@ void TKDTree<Index, Value>::MakeBoundaries(Value *range)
 // Build boundaries for each node.
 
 
 // Build boundaries for each node.
 
 
-       if(range) memcpy(fRange, range, 2*fNDim*sizeof(Value));
+       if(range) memcpy(fRange, range, fNDimm*sizeof(Value));
        // total number of nodes including terminal nodes
        // total number of nodes including terminal nodes
-       Int_t totNodes = fNnodes + GetNTerminalNodes();
-       fBoundaries = new Value[totNodes*2*fNDim];
-       printf("Allocate %d nodes\n", totNodes);
+       Int_t totNodes = fNnodes + GetNTNodes();
+       fBoundaries = new Value[totNodes*fNDimm];
+       //Info("MakeBoundaries(Value*)", Form("Allocate boundaries for %d nodes", totNodes));
+
+       
        // loop
        // loop
-       Int_t child_index;
+       Value *tbounds = 0x0, *cbounds = 0x0;
+       Int_t cn;
        for(int inode=fNnodes-1; inode>=0; inode--){
        for(int inode=fNnodes-1; inode>=0; inode--){
-               //printf("\tAllocate node %d\n", inode);
-               memcpy(&fBoundaries[inode*2*fNDim], fRange, 2*fNDim*sizeof(Value));
+               tbounds = &fBoundaries[inode*fNDimm];
+               memcpy(tbounds, fRange, fNDimm*sizeof(Value));
                                
                // check left child node
                                
                // check left child node
-               child_index = 2*inode+1;
-               if(IsTerminal(child_index)) CookBoundariesTerminal(inode, kTRUE);
-               for(int idim=0; idim<fNDim; idim++) fBoundaries[2*fNDim*inode+2*idim] = fBoundaries[2*fNDim*child_index+2*idim];
+               cn = (inode<<1)+1;
+               if(IsTerminal(cn)) CookBoundaries(inode, kTRUE);
+               cbounds = &fBoundaries[fNDimm*cn];
+               for(int idim=0; idim<fNDim; idim++) tbounds[idim<<1] = cbounds[idim<<1];
                
                // check right child node
                
                // check right child node
-               child_index = 2*inode+2;
-               if(IsTerminal(child_index)) CookBoundariesTerminal(inode, kFALSE);
-               for(int idim=0; idim<fNDim; idim++) fBoundaries[2*fNDim*inode+2*idim+1] = fBoundaries[2*fNDim*child_index+2*idim+1];
+               cn = (inode+1)<<1;
+               if(IsTerminal(cn)) CookBoundaries(inode, kFALSE);
+               cbounds = &fBoundaries[fNDimm*cn];
+               for(int idim=0; idim<fNDim; idim++) tbounds[(idim<<1)+1] = cbounds[(idim<<1)+1];
        }
 }
 
 //_________________________________________________________________
 template <typename Index, typename Value>
        }
 }
 
 //_________________________________________________________________
 template <typename Index, typename Value>
-void TKDTree<Index, Value>::CookBoundariesTerminal(Int_t parent_node, Bool_t LEFT)
+void TKDTree<Index, Value>::CookBoundaries(const Int_t node, Bool_t LEFT)
 {
        // define index of this terminal node
 {
        // define index of this terminal node
-       Int_t index = 2 * parent_node + (LEFT ? 1 : 2);
+       Int_t index = (node<<1) + (LEFT ? 1 : 2);
+       //Info("CookBoundaries()", Form("Node %d", index));
        
        // build and initialize boundaries for this node
        
        // build and initialize boundaries for this node
-       //printf("\tAllocate terminal node %d\n", index);
-       memcpy(&fBoundaries[2*fNDim*index], fRange, 2*fNDim*sizeof(Value));
+       Value *tbounds = &fBoundaries[fNDimm*index];
+       memcpy(tbounds, fRange, fNDimm*sizeof(Value));
        Bool_t flag[256];  // cope with up to 128 dimensions 
        Bool_t flag[256];  // cope with up to 128 dimensions 
-       memset(flag, kFALSE, 2*fNDim);
+       memset(flag, kFALSE, fNDimm);
        Int_t nvals = 0;
                
        // recurse parent nodes
        Int_t nvals = 0;
                
        // recurse parent nodes
-       TKDNode *node = 0x0;
-       while(parent_node >= 0 && nvals < 2 * fNDim){
-               node = &(fNodes[parent_node]);
+       Int_t pn = node;
+       while(pn >= 0 && nvals < fNDimm){
                if(LEFT){
                if(LEFT){
-                       if(!flag[2*node->fAxis+1]) {
-                               fBoundaries[2*fNDim*index+2*node->fAxis+1] = node->fValue;
-                               flag[2*node->fAxis+1] = kTRUE;
+                       index = (fAxis[pn]<<1)+1;
+                       if(!flag[index]) {
+                               tbounds[index] = fValue[pn];
+                               flag[index] = kTRUE;
                                nvals++;
                        }
                } else {
                                nvals++;
                        }
                } else {
-                       if(!flag[2*node->fAxis]) {
-                               fBoundaries[2*fNDim*index+2*node->fAxis] = node->fValue;
-                               flag[2*node->fAxis] = kTRUE;
+                       index = fAxis[pn]<<1;
+                       if(!flag[index]) {
+                               tbounds[index] = fValue[pn];
+                               flag[index] = kTRUE;
                                nvals++;
                        }
                }
                                nvals++;
                        }
                }
-               LEFT = parent_node%2;
-               parent_node =  parent_node/2 + (parent_node%2) - 1;
+               LEFT = pn&1;
+               pn =  (pn - 1)>>1;
        }
 }
 
        }
 }