Reverting changes in STAT waiting for fixes
authorfca <fca@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 11 Dec 2008 15:27:10 +0000 (15:27 +0000)
committerfca <fca@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 11 Dec 2008 15:27:10 +0000 (15:27 +0000)
STAT/CMake_libSTAT.txt
STAT/STATLinkDef.h
STAT/TKDTree.cxx [new file with mode: 0644]
STAT/TKDTree.h [new file with mode: 0644]
STAT/libSTAT.pkg

index fa7eaad..91e7bc2 100644 (file)
@@ -1,6 +1,7 @@
 # -*- mode: cmake -*-
 
 set(SRCS
+  TKDTree.cxx 
   TKDInterpolatorBase.cxx 
   TKDNodeInfo.cxx 
   TKDPDF.cxx 
index 5b09bac..3fa5253 100644 (file)
@@ -7,6 +7,12 @@
 #pragma link off all classes;
 #pragma link off all functions;
 
+#pragma link C++ class TKDTree<Int_t, Float_t>+;
+#pragma link C++ class TKDTree<Int_t, Double_t>+;
+#pragma link C++ typedef TKDTreeID;
+#pragma link C++ typedef TKDTreeIF;
+
+
 #pragma link C++ class TKDInterpolatorBase+;
 #pragma link C++ class TKDNodeInfo+;
 #pragma link C++ class TKDPDF+;
diff --git a/STAT/TKDTree.cxx b/STAT/TKDTree.cxx
new file mode 100644 (file)
index 0000000..907435a
--- /dev/null
@@ -0,0 +1,975 @@
+#include "TKDTree.h"
+#include "TRandom.h"
+
+#include "TString.h"
+#include <string.h>
+
+#ifndef R__ALPHA
+templateClassImp(TKDTree)
+#endif
+
+/*
+
+
+
+Content:
+1. What is kd-tree
+2. How to cosntruct kdtree - Pseudo code
+3. TKDTree implementation
+4. User interface
+
+
+
+1. What is kdtree ? ( http://en.wikipedia.org/wiki/Kd-tree )
+
+In computer science, a kd-tree (short for k-dimensional tree) is a space-partitioning data structure for organizing points in a k-dimensional space. kd-trees are a useful data structure for several applications, such as searches involving a multidimensional search key (e.g. range searches and nearest neighbour searches). kd-trees are a special case of BSP trees.
+
+A kd-tree uses only splitting planes that are perpendicular to one of the coordinate system axes. This differs from BSP trees, in which arbitrary splitting planes can be used. In addition, in the typical definition every node of a kd-tree, from the root to the leaves, stores a point.[1] This differs from BSP trees, in which leaves are typically the only nodes that contain points (or other geometric primitives). As a consequence, each splitting plane must go through one of the points in the kd-tree. kd-tries are a variant that store data only in leaf nodes. It is worth noting that in an alternative definition of kd-tree the points are stored in its leaf nodes only, although each splitting plane still goes through one of the points.[2]
+<img src=".....">
+
+2. Constructing a kd-tree ( Pseudo code)
+
+Since there are many possible ways to choose axis-aligned splitting planes, there are many different ways to construct kd-trees. The canonical method of kd-tree construction has the following constraints:
+
+    * As one moves down the tree, one cycles through the axes used to select the splitting planes. (For example, the root would have an x-aligned plane, the root's children would both have y-aligned planes, the root's grandchildren would all have z-aligned planes, and so on.)
+    * At each step, the point selected to create the splitting plane is the median of the points being put into the kd-tree, with respect to their coordinates in the axis being used. (Note the assumption that we feed the entire set of points into the algorithm up-front.)
+
+This method leads to a balanced kd-tree, in which each leaf node is about the same distance from the root. However, balanced trees are not necessarily optimal for all applications.
+
+Note also that it is not required to select the median point. In that case, the result is simply that there is no guarantee that the tree will be balanced. A simple heuristic to avoid coding a complex linear-time median-finding algorithm nor using an O(n log n) sort is to use sort to find the median of a fixed number of randomly selected points to serve as the cut line. Practically this technique results in nicely balanced trees.
+
+Given a list of n points, the following algorithm will construct a balanced kd-tree containing those points.
+<img src=".....">
+
+function kdtree (list of points pointList, int depth)
+{
+    if pointList is empty
+        return nil;
+    else
+    {
+        // Select axis based on depth so that axis cycles through all valid values
+        var int axis := depth mod k;
+
+        // Sort point list and choose median as pivot element
+        select median from pointList;
+
+        // Create node and construct subtrees
+        var tree_node node;
+        node.location := median;
+        node.leftChild := kdtree(points in pointList before median, depth+1);
+        node.rightChild := kdtree(points in pointList after median, depth+1);
+        return node;
+    }
+}
+
+
+3. TKDtree implementation. 
+
+// TKDtree is optimized to minimize memory consumption.
+// TKDTree is by default balanced. Nodes are mapped to the 1 dimensional arrays of size
+// fNnodes.
+// 
+// fAxix[fNnodes]  - Division axis (0,1,2,3 ...)
+// fValue[fNnodes] - Division value 
+//
+// Navigation to the Left and right doughter node is expressed by  analytical formula
+// Let's parent node id  is inode
+// Then: 
+// Left   = inode*2+1
+// Right  =  (inode+1)*2  
+// Let's  daughter node id is inode
+// Then:
+// Parent = inode/2
+// 
+// 
+// The mapping of the nodes to the 1D array enables us to use the indexing mechanism.
+// This is important if some additional kind of information 
+// (not directly connected to kd-tree, e.g function fit parameters in the node)
+// follow binary tree structure. The mapping can be reused.
+// 
+// Drawback:   Insertion to the TKDtree is not supported.  
+// Advantage:  Random access is supported -  simple formulas  
+// 
+// Number of division nodes and number of terminals :
+// fNnodes = (fNPoints/fBucketSize)
+// 
+// TKDTree is mapped to one dimensional array.
+// Mapping: (see for example the TKDTree::FindNode) 
+// 
+// The nodes are filled always from left side to the right side:
+//
+//
+// TERMINOLOGY: 
+// - inode  -  index of node 
+// - irow   -  index of row
+
+// Ideal case:
+// Number of terminal nodes = 2^N - N=3
+//
+//            INode
+// irow 0     0                                                                   -  1 inode 
+// irow 1     1                              2                                    -  2 inodes
+// irow 2     3              4               5               6                    -  4 inodes
+// irow 3     7       8      9      10       11     12       13      14           -  8 inodes
+//
+//
+// Non ideal case:
+// Number of terminal nodes = 2^N+k  - N=3  k=1
+
+//           INode
+// irow 0     0                                                                   - 1 inode
+// irow 1     1                              2                                    - 2 inodes
+// irow 2     3              4               5               6                    - 3 inodes
+// irow 3     7       8      9      10       11     12       13      14           - 8 inodes
+// irow 4     15  16                                                              - 2 inodes 
+//
+//
+// The division algorithm:
+// The  n points are divided to tho groups.
+//
+// n=2^k+rest
+//
+// Left  = 2^k-1 +(rest>2^k-2) ?  2^k-2      : rest
+// Right = 2^k-1 +(rest>2^k-2) ?  rest-2^k-2 : 0
+//
+//
+
+
+
+*/
+
+
+
+//////////////////////////////////////////////////////////////////////
+// Memory setup of protected data members:
+// 
+//     fDataOwner;  // Toggle ownership of the data
+//     fNnodes:
+//  size of branch node array, and index of first terminal node
+//     fNDim;       // number of variables
+//     fNpoints;    // number of multidimensional points
+//     fBucketSize; // number of data points per terminal node
+// 
+//     fIndPoints : array containing rearanged indexes according to nodes:
+//     | point indexes from 1st terminal node (fBucketSize elements) | ... | point indexes from the last terminal node (<= fBucketSize elements)
+// 
+//     Value   **fData;
+//     fRange : array containing the boundaries of the domain:
+//     | 1st dimension (min + max) | 2nd dimension (min + max) | ...
+//     fBoundaries : nodes boundaries
+//     | 1st node {1st dim * 2 elements | 2nd dim * 2 elements | ...} | 2nd node {...} | ...
+//     the nodes are arranged in the following order :
+//      - first fNnodes are primary nodes
+//      - after are the terminal nodes
+//     fNodes : array of primary nodes
+///////////////////////////////////////////////////////////////////////
+
+
+//_________________________________________________________________
+template <typename  Index, typename Value>
+TKDTree<Index, Value>::TKDTree() :
+  TObject()
+  ,fDataOwner(kFALSE)
+  ,fNnodes(0)
+  ,fNDim(0)
+  ,fNDimm(0)
+  ,fNpoints(0)
+  ,fBucketSize(0)
+  ,fAxis(0x0)
+  ,fValue(0x0)
+  ,fRange(0x0)
+  ,fData(0x0)
+  ,fBoundaries(0x0)
+  ,fkNNdim(0)
+  ,fkNN(0x0)
+  ,fkNNdist(0x0)
+  ,fDistBuffer(0x0)
+  ,fIndBuffer(0x0)
+  ,fIndPoints(0x0)
+  ,fRowT0(0)
+  ,fCrossNode(0)
+  ,fOffset(0)
+{
+// Default constructor. Nothing is built
+}
+
+
+//_________________________________________________________________
+template <typename  Index, typename Value>
+TKDTree<Index, Value>::TKDTree(Index npoints, Index ndim, UInt_t bsize, Value **data) :
+  TObject()
+  ,fDataOwner(kFALSE)
+  ,fNnodes(0)
+  ,fNDim(ndim)
+  ,fNDimm(2*ndim)
+  ,fNpoints(npoints)
+  ,fBucketSize(bsize)
+  ,fAxis(0x0)
+  ,fValue(0x0)
+  ,fRange(0x0)
+  ,fData(data) //Columnwise!!!!!
+  ,fBoundaries(0x0)
+  ,fkNNdim(0)
+  ,fkNN(0x0)
+  ,fkNNdist(0x0)
+  ,fDistBuffer(0x0)
+  ,fIndBuffer(0x0)
+  ,fIndPoints(0x0)
+  ,fRowT0(0)
+  ,fCrossNode(0)
+  ,fOffset(0)
+{
+  // Allocate data by hand. See TKDTree(TTree*, const Char_t*) for an automatic method.
+  
+  Build();
+}
+
+//_________________________________________________________________
+template <typename  Index, typename Value>
+TKDTree<Index, Value>::~TKDTree()
+{
+  // Destructor
+  if (fIndBuffer) delete [] fIndBuffer;
+  if (fDistBuffer) delete [] fDistBuffer;
+  if (fkNNdist) delete [] fkNNdist;
+  if (fkNN) delete [] fkNN;
+  if (fAxis) delete [] fAxis;
+  if (fValue) delete [] fValue;
+  if (fIndPoints) delete [] fIndPoints;
+  if (fRange) delete [] fRange;
+  if (fBoundaries) delete [] fBoundaries;
+  if (fDataOwner && fData){
+    for(int idim=0; idim<fNDim; idim++) delete [] fData[idim];
+    delete [] fData;
+  }
+}
+
+
+//_________________________________________________________________
+template <typename  Index, typename Value>
+void TKDTree<Index, Value>::Build(){
+  //
+  // Build binning
+  //
+  // 1. calculate number of nodes
+  // 2. calculate first terminal row
+  // 3. initialize index array
+  // 4. non recursive building of the binary tree
+
+  //
+  // Teh tree is recursivaly divided. The division point is choosen in such a way, that the left side 
+  // alway has 2^k points, while at the same time trying 
+  //
+  //1.
+  fNnodes = fNpoints/fBucketSize-1;
+  if (fNpoints%fBucketSize) fNnodes++;
+  
+  //2.
+  fRowT0=0;
+  for (;(fNnodes+1)>(1<<fRowT0);fRowT0++);
+  fRowT0-=1;
+  //         2 = 2**0 + 1
+  //         3 = 2**1 + 1
+  //         4 = 2**1 + 2
+  //         5 = 2**2 + 1
+  //         6 = 2**2 + 2
+  //         7 = 2**2 + 3
+  //         8 = 2**2 + 4
+  
+  //3.
+  // allocate space for boundaries
+  fRange = new Value[2*fNDim];
+  fIndPoints= new Index[fNpoints];
+  for (Index i=0; i<fNpoints; i++) fIndPoints[i] = i;
+  fAxis  = new UChar_t[fNnodes];
+  fValue = new Value[fNnodes];
+  //
+  fCrossNode = (1<<(fRowT0+1))-1;
+  if (fCrossNode<fNnodes) fCrossNode = 2*fCrossNode+1;
+  //
+  //  fOffset = (((fNnodes+1)-(1<<fRowT0)))*2;
+  Int_t   over   = (fNnodes+1)-(1<<fRowT0);
+  Int_t   filled = ((1<<fRowT0)-over)*fBucketSize;
+  fOffset = fNpoints-filled;
+  //
+  //   printf("Row0      %d\n", fRowT0);
+  //   printf("CrossNode %d\n", fCrossNode);
+  //   printf("Offset    %d\n", fOffset);
+  //
+  //
+  //4.
+  //    stack for non recursive build - size 128 bytes enough
+  Int_t rowStack[128];
+  Int_t nodeStack[128];
+  Int_t npointStack[128];
+  Int_t posStack[128];
+  Int_t currentIndex = 0;
+  Int_t iter =0;
+  rowStack[0]    = 0;
+  nodeStack[0]   = 0;
+  npointStack[0] = fNpoints;
+  posStack[0]   = 0;
+  //
+  Int_t nbucketsall =0;
+  while (currentIndex>=0){
+    iter++;
+    //
+    Int_t npoints  = npointStack[currentIndex];
+    if (npoints<=fBucketSize) {
+      //printf("terminal node : index %d iter %d\n", currentIndex, iter);
+      currentIndex--;
+      nbucketsall++;
+      continue; // terminal node
+    }
+    Int_t crow     = rowStack[currentIndex];
+    Int_t cpos     = posStack[currentIndex];
+    Int_t cnode    = nodeStack[currentIndex];          
+    //printf("currentIndex %d npoints %d node %d\n", currentIndex, npoints, cnode);
+    //
+    // divide points
+    Int_t nbuckets0 = npoints/fBucketSize;           //current number of  buckets
+    if (npoints%fBucketSize) nbuckets0++;            //
+    Int_t restRows = fRowT0-rowStack[currentIndex];  // rest of fully occupied node row
+    if (restRows<0) restRows =0;
+    for (;nbuckets0>(2<<restRows); restRows++);
+    Int_t nfull = 1<<restRows;
+    Int_t nrest = nbuckets0-nfull;
+    Int_t nleft =0, nright =0;
+    //
+    if (nrest>(nfull/2)){
+      nleft  = nfull*fBucketSize;
+      nright = npoints-nleft;
+    }else{
+      nright = nfull*fBucketSize/2;
+      nleft  = npoints-nright;
+    }
+    
+    //
+    //find the axis with biggest spread
+    Value maxspread=0;
+    Value tempspread, min, max;
+    Index axspread=0;
+    Value *array;
+    for (Int_t idim=0; idim<fNDim; idim++){
+      array = fData[idim];
+      Spread(npoints, array, fIndPoints+cpos, min, max);
+      tempspread = max - min;
+      if (maxspread < tempspread) {
+       maxspread=tempspread;
+       axspread = idim;
+      }
+      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);
+    fAxis[cnode]  = axspread;
+    fValue[cnode] = array[fIndPoints[cpos+nleft]];
+    //printf("Set node %d : ax %d val %f\n", cnode, node->fAxis, node->fValue);
+    //
+    //
+    npointStack[currentIndex] = nleft;
+    rowStack[currentIndex]    = crow+1;
+    posStack[currentIndex]    = cpos;
+    nodeStack[currentIndex]   = cnode*2+1;
+    currentIndex++;
+    npointStack[currentIndex] = nright;
+    rowStack[currentIndex]    = crow+1;
+    posStack[currentIndex]    = cpos+nleft;
+    nodeStack[currentIndex]   = (cnode*2)+2;
+    //
+    if (0){
+      // consistency check
+      Info("Build()", Form("points %d left %d right %d", npoints, nleft, nright));
+      if (nleft<nright) Warning("Build", "Problem Left-Right");
+      if (nleft<0 || nright<0) Warning("Build()", "Problem Negative number");
+    }
+  }    
+}
+
+
+//_________________________________________________________________
+template <typename  Index, typename Value>
+Bool_t TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int_t k, Index *&in, Value *&d)
+{
+  //
+  // NOT MI - to be checked
+  //
+  //
+  // Find "k" nearest neighbors to "point".
+  //
+  // Return true in case of success and false in case of failure.
+  // The indexes of the nearest k points are stored in the array "in" in
+  // increasing order of the distance to "point" and the maxim distance
+  // in "d".
+  //
+  // The arrays "in" and "d" are managed internally by the TKDTree.
+  
+  Index inode = FindNode(point);
+  if(inode < fNnodes){
+    Error("FindNearestNeighbors()", "Point outside data range.");
+    return kFALSE;
+  }
+  
+  UInt_t debug = 0;
+  Int_t nCalculations = 0;
+  
+  // allocate working memory
+  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.      
+  if(!fBoundaries) MakeBoundaries();
+  Value *bounds = &fBoundaries[inode*2*fNDim];
+  Index nBounds = 0;
+  for(int idim=0; idim<fNDim; idim++){
+    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
+  UChar_t ax;
+  Value   val, dif;
+  Int_t nAllNodes = fNnodes + fNpoints/fBucketSize + ((fNpoints%fBucketSize)?1:0);
+  Int_t nodeStack[128], nodeIn[128];
+  Index currentIndex = 0;
+  nodeStack[0]   = inode;
+  nodeIn[0] = inode;
+  while(currentIndex>=0){
+    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
+    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++;
+      // end all recursions
+      if(nBounds==2 * fNDim) break;
+      continue;
+    }
+    
+    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];
+      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
+       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++;
+       }
+      }
+      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
+    Int_t pn = (tnode-1)>>1;//tnode/2 + (tnode%2) - 1;
+    if(pn >= 0 && entry != pn){
+      // check if parent node is eligible at all
+      if(TMath::Abs(point[fAxis[pn]] - fValue[pn]) > fkNNdist[k-1]){
+       // mark bound
+       nBounds++;
+       // end all recursions
+       if(nBounds==2 * fNDim) break;
+      }
+      currentIndex++;
+      nodeStack[currentIndex]=pn;
+      nodeIn[currentIndex]=tnode;
+      if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]);
+    }
+    if(IsTerminal(tnode)) continue;
+    
+    // register children nodes
+    Int_t cn;
+    Bool_t kAllow[] = {kTRUE, kTRUE};
+    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;
+      cn = (tnode<<1)+ic;
+      if(cn < nAllNodes && entry != cn){
+       currentIndex++;
+       nodeStack[currentIndex] = cn;
+       nodeIn[currentIndex]=tnode;
+       if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]);
+      }
+    }          
+  }
+  // save results
+  in = fkNN;
+  d  = fkNNdist;  
+  return kTRUE;
+}
+
+
+
+//_________________________________________________________________
+template <typename  Index, typename Value>
+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;
+  while (currentIndex>=0){
+    inode    = stackNode[currentIndex];
+    if (IsTerminal(inode)) return inode;
+    
+    currentIndex--;
+    if (point[fAxis[inode]]<=fValue[inode]){
+      currentIndex++;
+      stackNode[currentIndex]=(inode<<1)+1;
+    }
+    if (point[fAxis[inode]]>=fValue[inode]){
+      currentIndex++;
+      stackNode[currentIndex]=(inode+1)<<1;
+    }
+  }
+  
+  return -1;
+}
+
+
+
+//_________________________________________________________________
+template <typename  Index, typename Value>
+void TKDTree<Index, Value>::FindPoint(Value * point, Index &index, Int_t &iter){
+  //
+  // find the index of point
+  // works only if we keep fData pointers
+  
+  Int_t stackNode[128];
+  Int_t currentIndex =0;
+  stackNode[0] = 0;
+  iter =0;
+  //
+  while (currentIndex>=0){
+    iter++;
+    Int_t inode    = stackNode[currentIndex];
+    currentIndex--;
+    if (IsTerminal(inode)){
+      // investigate terminal node
+      Int_t indexIP  = (inode >= fCrossNode) ? (inode-fCrossNode)*fBucketSize : (inode-fNnodes)*fBucketSize+fOffset;
+      printf("terminal %d indexP %d\n", inode, indexIP);
+      for (Int_t ibucket=0;ibucket<fBucketSize;ibucket++){
+       Bool_t isOK    = kTRUE;
+       indexIP+=ibucket;
+       printf("ibucket %d index %d\n", ibucket, indexIP);
+       if (indexIP>=fNpoints) continue;
+       Int_t index0   = fIndPoints[indexIP];
+       for (Int_t idim=0;idim<fNDim;idim++) if (fData[idim][index0]!=point[idim]) isOK = kFALSE;
+       if (isOK) index = index0;
+      }
+      continue;
+    }
+    
+    if (point[fAxis[inode]]<=fValue[inode]){
+      currentIndex++;
+      stackNode[currentIndex]=(inode*2)+1;
+    }
+    if (point[fAxis[inode]]>=fValue[inode]){
+      currentIndex++;
+      stackNode[currentIndex]=(inode*2)+2;
+    }
+  }
+  //
+  //  printf("Iter\t%d\n",iter);
+}
+
+//_________________________________________________________________
+template <typename  Index, typename Value>
+void TKDTree<Index, Value>::FindInRangeA(Value * point, Value * delta, Index *res , Index &npoints, Index & iter, Int_t bnode)
+{
+ //
+  // Find all points in the range specified by    (point +- range)
+  // res     - Resulting indexes are stored in res array 
+  // npoints - Number of selected indexes in range
+  // NOTICE:
+  // For some cases it is better to don't keep data - because of memory consumption
+  // If the data are not kept - only boundary conditions are investigated
+  // some of the data can be outside of the range
+  // What is guranteed in this mode: All of the points in the range are selected + some fraction of others (but close)
+
+       Index stackNode[128];
+  iter=0;
+  Index currentIndex = 0;
+  stackNode[currentIndex] = bnode; 
+  while (currentIndex>=0){
+    iter++;
+    Int_t inode    = stackNode[currentIndex];
+    //
+    currentIndex--;
+    if (!IsTerminal(inode)){
+      // not terminal
+      //TKDNode<Index, Value> * node = &(fNodes[inode]);
+      if (point[fAxis[inode]] - delta[fAxis[inode]] <  fValue[inode]) {
+       currentIndex++; 
+       stackNode[currentIndex]= (inode*2)+1;
+      }
+      if (point[fAxis[inode]] + delta[fAxis[inode]] >= fValue[inode]){
+       currentIndex++; 
+       stackNode[currentIndex]= (inode*2)+2;
+      }
+    }else{
+      Int_t indexIP  = 
+       (inode >= fCrossNode) ? (inode-fCrossNode)*fBucketSize : (inode-fNnodes)*fBucketSize+fOffset;
+      for (Int_t ibucket=0;ibucket<fBucketSize;ibucket++){
+       if (indexIP+ibucket>=fNpoints) break;
+       res[npoints]   = fIndPoints[indexIP+ibucket];
+       npoints++;
+      }
+    }
+  }
+  if (fData){
+    //
+    //  compress rest if data still accesible
+    //
+    Index npoints2 = npoints;
+    npoints=0;
+    for (Index i=0; i<npoints2;i++){
+      Bool_t isOK = kTRUE;
+      for (Index idim = 0; idim< fNDim; idim++){
+       if (TMath::Abs(fData[idim][res[i]]- point[idim])>delta[idim]) 
+         isOK = kFALSE;        
+      }
+      if (isOK){
+       res[npoints] = res[i];
+       npoints++;
+      }
+    }    
+  }
+}
+
+
+//_________________________________________________________________
+template <typename  Index, typename Value>
+void TKDTree<Index, Value>::FindInRangeB(Value * point, Value * delta, Index *res , Index &npoints,Index & iter, Int_t bnode)
+{
+  //
+  Long64_t  goldStatus = (1<<(2*fNDim))-1;  // gold status
+  Index stackNode[128];
+  Long64_t   stackStatus[128]; 
+  iter=0;
+  Index currentIndex   = 0;
+  stackNode[currentIndex]   = bnode; 
+  stackStatus[currentIndex] = 0;
+  while (currentIndex>=0){
+    Int_t inode     = stackNode[currentIndex];
+    Long64_t status = stackStatus[currentIndex];
+    currentIndex--;
+    iter++;
+    if (IsTerminal(inode)){
+      Int_t indexIP  = 
+       (inode >= fCrossNode) ? (inode-fCrossNode)*fBucketSize : (inode-fNnodes)*fBucketSize+fOffset;
+      for (Int_t ibucket=0;ibucket<fBucketSize;ibucket++){
+       if (indexIP+ibucket>=fNpoints) break;
+       res[npoints]   = fIndPoints[indexIP+ibucket];
+       npoints++;
+      }
+      continue;
+    }      
+    // not terminal    
+    if (status == goldStatus){
+      Int_t ileft  = inode;
+      Int_t iright = inode;
+      for (;ileft<fNnodes; ileft   = (ileft<<1)+1);
+      for (;iright<fNnodes; iright = (iright<<1)+2);
+      Int_t indexL  = 
+       (ileft >= fCrossNode) ? (ileft-fCrossNode)*fBucketSize : (ileft-fNnodes)*fBucketSize+fOffset;
+      Int_t indexR  = 
+       (iright >= fCrossNode) ? (iright-fCrossNode)*fBucketSize : (iright-fNnodes)*fBucketSize+fOffset;
+      if (indexL<=indexR){
+       Int_t endpoint = indexR+fBucketSize;
+       if (endpoint>fNpoints) endpoint=fNpoints;
+       for (Int_t ipoint=indexL;ipoint<endpoint;ipoint++){
+         res[npoints]   = fIndPoints[ipoint];
+         npoints++;
+       }         
+       continue;
+      }
+    }
+    //
+    // TKDNode<Index, Value> * node = &(fNodes[inode]);
+    if (point[fAxis[inode]] - delta[fAxis[inode]] <  fValue[inode]) {
+      currentIndex++; 
+      stackNode[currentIndex]= (inode<<1)+1;
+      if (point[fAxis[inode]] + delta[fAxis[inode]] > fValue[inode])
+       stackStatus[currentIndex]= status | (1<<(2*fAxis[inode]));
+    }
+    if (point[fAxis[inode]] + delta[fAxis[inode]] >= fValue[inode]){
+      currentIndex++; 
+      stackNode[currentIndex]= (inode<<1)+2;
+      if (point[fAxis[inode]] - delta[fAxis[inode]]<fValue[inode])
+       stackStatus[currentIndex]= status | (1<<(2*fAxis[inode]+1));
+    }
+  }
+  if (fData){
+    // compress rest
+    Index npoints2 = npoints;
+    npoints=0;
+    for (Index i=0; i<npoints2;i++){
+      Bool_t isOK = kTRUE;
+      for (Index idim = 0; idim< fNDim; idim++){
+       if (TMath::Abs(fData[idim][res[i]]- point[idim])>delta[idim]) 
+         isOK = kFALSE;        
+      }
+      if (isOK){
+       res[npoints] = res[i];
+       npoints++;
+      }
+    }    
+  }
+}
+
+
+
+//_________________________________________________________________
+template <typename  Index, typename Value>
+void TKDTree<Index, Value>::SetData(Index npoints, Index ndim, UInt_t bsize, Value **data)
+{
+// TO DO
+// 
+// Check reconstruction/reallocation of memory of data. Maybe it is not
+// necessary to delete and realocate space but only to use the same space
+
+       Clear();
+  
+       //Columnwise!!!!
+   fData = data;
+   fNpoints = npoints;
+   fNDim = ndim;
+   fBucketSize = bsize;
+
+       Build();
+}
+
+
+
+//_________________________________________________________________
+template <typename  Index, typename Value>
+void TKDTree<Index, Value>::Spread(Index ntotal, Value *a, Index *index, Value &min, Value &max) const
+{
+  //Value min, max;
+  Index i;
+  min = a[index[0]];
+  max = a[index[0]];
+  for (i=0; i<ntotal; i++){
+    if (a[index[i]]<min) min = a[index[i]];
+    if (a[index[i]]>max) max = a[index[i]];
+  }
+}
+
+
+//_________________________________________________________________
+template <typename  Index, typename Value>
+Value TKDTree<Index, Value>::KOrdStat(Index ntotal, Value *a, Index k, Index *index) const
+{
+  //
+   //copy of the TMath::KOrdStat because I need an Index work array
+
+   Index i, ir, j, l, mid;
+   Index arr;
+   Index temp;
+
+   Index rk = k;
+   l=0;
+   ir = ntotal-1;
+  for(;;) {
+      if (ir<=l+1) { //active partition contains 1 or 2 elements
+         if (ir == l+1 && a[index[ir]]<a[index[l]])
+           {temp = index[l]; index[l]=index[ir]; index[ir]=temp;}
+         Value tmp = a[index[rk]];
+         return tmp;
+      } else {
+         mid = (l+ir) >> 1; //choose median of left, center and right
+         {temp = index[mid]; index[mid]=index[l+1]; index[l+1]=temp;}//elements as partitioning element arr.
+         if (a[index[l]]>a[index[ir]])  //also rearrange so that a[l]<=a[l+1]
+           {temp = index[l]; index[l]=index[ir]; index[ir]=temp;}
+
+         if (a[index[l+1]]>a[index[ir]])
+           {temp=index[l+1]; index[l+1]=index[ir]; index[ir]=temp;}
+
+         if (a[index[l]]>a[index[l+1]])
+           {temp = index[l]; index[l]=index[l+1]; index[l+1]=temp;}
+
+         i=l+1;        //initialize pointers for partitioning
+         j=ir;
+         arr = index[l+1];
+         for (;;) {
+           do i++; while (a[index[i]]<a[arr]);
+           do j--; while (a[index[j]]>a[arr]);
+           if (j<i) break;  //pointers crossed, partitioning complete
+              {temp=index[i]; index[i]=index[j]; index[j]=temp;}
+         }
+         index[l+1]=index[j];
+         index[j]=arr;
+         if (j>=rk) ir = j-1; //keep active the partition that
+         if (j<=rk) l=i;      //contains the k_th element
+      }
+   }
+}
+
+//_________________________________________________________________
+template <typename Index, typename Value>
+void TKDTree<Index, Value>::MakeBoundaries(Value *range)
+{
+// Build boundaries for each node.
+
+
+       if(range) memcpy(fRange, range, fNDimm*sizeof(Value));
+       // total number of nodes including terminal nodes
+       Int_t totNodes = fNnodes + fNpoints/fBucketSize + ((fNpoints%fBucketSize)?1:0);
+       fBoundaries = new Value[totNodes*fNDimm];
+       //Info("MakeBoundaries(Value*)", Form("Allocate boundaries for %d nodes", totNodes));
+
+       
+       // loop
+       Value *tbounds = 0x0, *cbounds = 0x0;
+       Int_t cn;
+       for(int inode=fNnodes-1; inode>=0; inode--){
+               tbounds = &fBoundaries[inode*fNDimm];
+               memcpy(tbounds, fRange, fNDimm*sizeof(Value));
+                               
+               // check left child node
+               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
+               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>
+void TKDTree<Index, Value>::CookBoundaries(const Int_t node, Bool_t LEFT)
+{
+       // define index of this terminal node
+       Int_t index = (node<<1) + (LEFT ? 1 : 2);
+       //Info("CookBoundaries()", Form("Node %d", index));
+       
+       // build and initialize boundaries for this node
+       Value *tbounds = &fBoundaries[fNDimm*index];
+       memcpy(tbounds, fRange, fNDimm*sizeof(Value));
+       Bool_t flag[256];  // cope with up to 128 dimensions 
+       memset(flag, kFALSE, fNDimm);
+       Int_t nvals = 0;
+               
+       // recurse parent nodes
+       Int_t pn = node;
+       while(pn >= 0 && nvals < fNDimm){
+               if(LEFT){
+                       index = (fAxis[pn]<<1)+1;
+                       if(!flag[index]) {
+                               tbounds[index] = fValue[pn];
+                               flag[index] = kTRUE;
+                               nvals++;
+                       }
+               } else {
+                       index = fAxis[pn]<<1;
+                       if(!flag[index]) {
+                               tbounds[index] = fValue[pn];
+                               flag[index] = kTRUE;
+                               nvals++;
+                       }
+               }
+               LEFT = pn&1;
+               pn =  (pn - 1)>>1;
+       }
+}
+
+
+
+//______________________________________________________________________
+TKDTreeIF *TKDTreeTestBuild(const Int_t npoints, const Int_t bsize){  
+  //
+  // Example function to  
+  //
+  Float_t *data0 =  new Float_t[npoints*2];
+  Float_t *data[2];
+  data[0] = &data0[0];
+  data[1] = &data0[npoints];
+  for (Int_t i=0;i<npoints;i++) {
+    data[1][i]= gRandom->Rndm();
+    data[0][i]= gRandom->Rndm();
+  }
+  TKDTree<Int_t, Float_t> *kdtree = new TKDTreeIF(npoints, 2, bsize, data);
+  return kdtree;       
+}
+
+
+
+template class TKDTree<Int_t, Float_t>;
+template class TKDTree<Int_t, Double_t>;
+
diff --git a/STAT/TKDTree.h b/STAT/TKDTree.h
new file mode 100644 (file)
index 0000000..4b65c44
--- /dev/null
@@ -0,0 +1,130 @@
+#ifndef ROOT_TKDTree
+#define ROOT_TKDTree
+
+#ifndef ROOT_TObject
+#include "TObject.h"
+#endif
+
+#include "TMath.h"
+
+
+template <typename Index, typename Value> class TKDTree : public TObject
+{
+public:
+       
+  TKDTree();
+  TKDTree(Index npoints, Index ndim, UInt_t bsize, Value **data);
+  ~TKDTree();
+  // Get indexes of left and right daughter nodes
+  Int_t GetLeft(Int_t inode)  const    {return inode*2+1;}
+  Int_t GetRight(Int_t inode) const    {return (inode+1)*2;}
+  Int_t GetParent(Int_t inode) const  {return inode/2;}
+  //  
+  // getters
+  Index*  GetPointsIndexes(Int_t node) const {
+    if(node < fNnodes) return 0x0;
+    Int_t offset = (node >= fCrossNode) ? (node-fCrossNode)*fBucketSize : fOffset+(node-fNnodes)*fBucketSize;
+    return &fIndPoints[offset];
+  }
+  UChar_t GetNodeAxis(Int_t id) const {return (id < 0 || id >= fNnodes) ? 0 : fAxis[id];}
+  Value   GetNodeValue(Int_t id) const {return (id < 0 || id >= fNnodes) ? 0 : fValue[id];}
+  Int_t   GetNNodes() const {return fNnodes;}
+  Value*  GetBoundaries();
+  Value*  GetBoundary(const Int_t node);
+  static  Int_t   GetIndex(Int_t row, Int_t collumn){return collumn+(1<<row);}
+  static void GetCoord(Int_t index, Int_t &row, Int_t &collumn){for (row=0; index>=(16<<row);row+=4); for (; index>=(2<<row);row++);collumn= index-(1<<row);};
+  Bool_t  FindNearestNeighbors(const Value *point, const Int_t kNN, Index *&i, Value *&d);
+  Index   FindNode(const Value * point);
+  void    FindPoint(Value * point, Index &index, Int_t &iter);
+  void    FindInRangeA(Value * point, Value * delta, Index *res , Index &npoints,Index & iter, Int_t bnode);
+  void    FindInRangeB(Value * point, Value * delta, Index *res , Index &npoints,Index & iter, Int_t bnode);
+  void    FindBNodeA(Value * point, Value * delta, Int_t &inode);
+  Bool_t  IsTerminal(Index inode) const {return (inode>=fNnodes);}
+  Value   KOrdStat(Index ntotal, Value *a, Index k, Index *index) const;
+  void    MakeBoundaries(Value *range = 0x0);
+  void    SetData(Index npoints, Index ndim, UInt_t bsize, Value **data);
+  void    Spread(Index ntotal, Value *a, Index *index, Value &min, Value &max) const;
+  
+
+ protected:
+  void            Build();  // build tree
+  
+ private:
+  TKDTree(const TKDTree &); // not implemented
+  TKDTree<Index, Value>& operator=(const TKDTree<Index, Value>&); // not implemented
+  void            CookBoundaries(const Int_t node, Bool_t left);
+  
+  
+ protected:
+  Bool_t  fDataOwner;  //! Toggle ownership of the data
+  Int_t   fNnodes;     // size of node array
+  Index   fNDim;       // number of dimensions
+  Index   fNDimm;      // dummy 2*fNDim
+  Index   fNpoints;    // number of multidimensional points
+  Index   fBucketSize; // limit statistic for nodes 
+  UChar_t *fAxis;      //[fNnodes] nodes cutting axis
+  Value   *fValue;     //[fNnodes] nodes cutting value
+  //
+  Value          *fRange;     //[fNDimm] range of data for each dimension
+  Value   **fData;     //! data points
+  Value          *fBoundaries;//! nodes boundaries - check class doc
+  
+  // kNN related data
+ protected:
+  Int_t   fkNNdim;     //! current kNN arrays allocated dimension
+  Index   *fkNN;       //! k nearest neighbors indexes
+  Value   *fkNNdist;   //! k nearest neighbors distances
+  Value   *fDistBuffer;//! working space for kNN
+  Index   *fIndBuffer; //! working space for kNN
+  
+ private:
+  Index   *fIndPoints; //! array of points indexes
+  Int_t   fRowT0;      //! smallest terminal row
+  Int_t   fCrossNode;  //! cross node
+  Int_t   fOffset;     //! offset in fIndPoints
+  
+  ClassDef(TKDTree, 1)  // KD tree
+};
+
+
+typedef TKDTree<Int_t, Double_t> TKDTreeID;
+typedef TKDTree<Int_t, Float_t> TKDTreeIF;
+
+// Test functions:
+TKDTreeIF *  TKDTreeTestBuild();
+
+
+
+//_________________________________________________________________
+template <typename  Index, typename Value>
+void TKDTree<Index, Value>::FindBNodeA(Value *point, Value *delta, Int_t &inode){
+  //
+  // find the smallest node covering the full range - start
+  //
+  inode =0; 
+  for (;inode<fNnodes;){
+    if (TMath::Abs(point[fAxis[inode]] - fValue[inode])<delta[fAxis[inode]]) break;
+    inode = (point[fAxis[inode]] < fValue[inode]) ? (inode*2)+1: (inode*2)+2;
+  }
+}
+
+//_________________________________________________________________
+template <typename  Index, typename Value>
+Value* TKDTree<Index, Value>::GetBoundaries()
+{
+  // Get the boundaries
+       if(!fBoundaries) MakeBoundaries();
+       return fBoundaries;
+}
+
+//_________________________________________________________________
+template <typename  Index, typename Value>
+Value* TKDTree<Index, Value>::GetBoundary(const Int_t node)
+{
+  // Get a boundary
+       if(!fBoundaries) MakeBoundaries();
+       return &fBoundaries[node*2*fNDim];
+}
+
+#endif
+
index 859ee71..dbfc720 100644 (file)
@@ -1,6 +1,6 @@
 #-*- Mode: Makefile -*-
 
-SRCS=  \
+SRCS=  TKDTree.cxx \
   TKDInterpolatorBase.cxx \
   TKDNodeInfo.cxx \
   TKDPDF.cxx \